Hi Jean-Paul,
> The indexing is correct - I've taken that part of the code from some
other codes that I've fully verified. But I thought it would still be a
useful exercise to re-implement the assembly routine for the other two
parameterisations. Maybe this would further convince you of its
correctness.
Convinced ;)
> I've attached a quickly modified version of step-44 that does this (see
the functions assemble_system_one_cell_tau(), assemble_system_one_cell_S()
and assemble_system_one_cell_P()), and uses the aforementioned push/pull
operations to transform the stresses and material tangents. (I didn't feel
like reimplementing the constitutive laws for the different
parameterisations -- I leave that as an exercise to you, if you feel so
inclined.) I've also attached some result logs for the default
configuration for step-44 plus one additional global refinement step.
ah very nice, thanks again.
> That's great! To me, this implies that you fixed a bug. Ultimately the
RHS that is computed via any of these parameterisations should, in
principle, lead to exactly the same result. They are equivalent, and it is
sometimes simply convenience (in terms of material laws, numerical
efficiency,...) that dictates which you might choose. They are all
expressing the weak form of the balance of linear momentum (i.e. the
different power conjugate pairs collectively quantify the same mechanical
power), and can happily be transformed from one to the other by exploiting
the relationships between the various stresses and strains. This implies
that the exact linearisation, even if its expressed in terms of some other
stress-strain measures, will in fact be the linearisation of any one of
these three expressions of the residual associated with displacement DoFs.
The linearisation can be similarly computed for one parameterisation and
transformed into the others. Its a nice exercise to go though, if you have
the time and patience to do so.
Ah ok, I was not aware of that! In your first answer, you wrote "So you may
as well start off by parameterising the problem in terms of F and P, and
naturally you have to adjust the linearisation as well." I thought I need
to adjust the linearization in any case, when I start to rephrase the
residual assembly, but I guess you talked about an entirely rephrased
problem (i.e., residual and linearization), which requires an adaption of
the linearization. I checked now the convergence behavior of my two-point
residual assembly in combination with the spatial linearization and I
observe indeed quadratic convergence (though there are some deviations
depending on the particular iteration). Rephrasing the linearization is
therefore as an exercise left for the next pandemic.
One last point I don't understand right now is the magnitude of the machine
precision you talk about:
> You'll see that they all exhibit quadratic convergence, with only a
slight difference in the convergence history that can probably be
attributed to accumulated round-off errors. Nevertheless, for all
implementations the absolute residual at the end of each timestep can be
pushed down to near machine precision with only a few Newton iterations.
The absolute errors at the end of each iteration are of the order \mathcal
{O}(10^{-10}). Also, if I compare my two-point residual assembly using
matrix-free with the 'usual' spatial residual assembly and compute the l2
norm of the difference, i.e. l2_norm(r_mf - r-_tau), I obtain a difference
of the order \mathcal{O}(10^{-10}). But from what I learnt, double
precision should be at least accurate up to ~10^{-15}. There are five
orders of magnitude in between. I probably miss something here in between.
Any idea?
Regards,
David
Jean-Paul Pelteret schrieb am Mittwoch, 30. Dezember 2020 um 21:27:30 UTC+1:
> HI David,
>
> You're welcome! I'm glad that you found the explanation useful :-)
>
>
> the required changes for the tangent operator seem acceptable (assuming
> the indexing in the get_HH() function does the right thing, otherwise it's
> probably not acceptable)
>
>
> The indexing is correct - I've taken that part of the code from some other
> codes that I've fully verified. But I thought it would still be a useful
> exercise to re-implement the assembly routine for the other two
> parameterisations. Maybe this would further convince you of its
> correctness.
>
> I've attached a quickly modified version of step-44 that does this (see
> the functions assemble_system_one_cell_tau(), assemble_system_one_cell_S()
> and assemble_system_one_cell_P()), and uses the aforementioned push/pull
> operations to transform the stresses and material tangents. (I didn't feel
> like reimplementing the constitutive laws for the different
> parameterisations -- I leave that as an exercise to you, if you feel so
> inclined.) I've also attached some result logs for the default
> configuration for step-44 plus one additional global refinement step.
> You'll see that they all exhibit quadratic convergence, with only a slight
> difference in the convergence history that can probably be attributed to
> accumulated round-off errors. Nevertheless, for all implementations the
> absolute residual at the end of each timestep can be pushed down to near
> machine precision with only a few Newton iterations.
>
> Comparing the P-formulation to the rest, you can see that its actually
> quite straight forward to implement, so hopefully converting it to matrix
> free wouldn't pose too much of a challenge. Note that in both the P- and
> S-formulations, all of the linearisation terms coupled to the displacement
> also change slightly. The coupling terms derive from the definition of the
> hydrostatic stress, which is expressed differently for each
> parameterisation.
>
>
> As an interesting side note: If I remove my wrong 'symm_grad_Nx'
> contribution from my shared code snippet, the remaining code corresponds
> already to a rephrased formulation using the Kirchoff stress and
> deformation gradient (as you have already noticed). Using this residual in
> combination with the "tau-Jc" tangent from step-44 (one-field adopted)
> leads to a convergence of the nonlinear scheme. So, it seems the
> linearization is still "correct enough" for the modified (two-point) frame
> of the residual. I will -of course- not rely on this half baked solution.
>
>
> That's great! To me, this implies that you fixed a bug. Ultimately the RHS
> that is computed via any of these parameterisations should, in principle,
> lead to exactly the same result. They are equivalent, and it is sometimes
> simply convenience (in terms of material laws, numerical efficiency,...)
> that dictates which you might choose. They are all expressing the weak form
> of the balance of linear momentum (i.e. the different power conjugate pairs
> collectively quantify the same mechanical power), and can happily be
> transformed from one to the other by exploiting the relationships between
> the various stresses and strains. This implies that the exact
> linearisation, even if its expressed in terms of some other stress-strain
> measures, will in fact be the linearisation of any one of these three
> expressions of the residual associated with displacement DoFs. The
> linearisation can be similarly computed for one parameterisation and
> transformed into the others. Its a nice exercise to go though, if you have
> the time and patience to do so.
>
> Best,
> Jean-Paul
>
>
> On 30.12.20 15:07, 'David' via deal.II User Group wrote:
>
> Hi Jean-Paul,
>
> first of all, thank you very much for your comprehensive answer. It's a
> little bit unfortunate that there is no way for the deformed test gradient
> evaluation using matrix-free. However, your code snippets and the
> explanation are traceable and the required changes for the tangent operator
> seem acceptable (assuming the indexing in the get_HH() function does the
> right thing, otherwise it's probably not acceptable). As always, I was a
> little bit too ambitious with my time schedule for this task. I will
> document your suggestions in an issue and try to implement the rephrased
> problem in the medium run. In the end, the change is entirely performance
> motivated, since the residual assembly (without matrx-free) is
> comparatively expensive considering that it is only performed once per
> nonlinear iteration.
>
> As an interesting side note: If I remove my wrong 'symm_grad_Nx'
> contribution from my shared code snippet, the remaining code corresponds
> already to a rephrased formulation using the Kirchoff stress and
> deformation gradient (as you have already noticed). Using this residual in
> combination with the "tau-Jc" tangent from step-44 (one-field adopted)
> leads to a convergence of the nonlinear scheme. So, it seems the
> linearization is still "correct enough" for the modified (two-point) frame
> of the residual. I will -of course- not rely on this half baked solution.
>
> Best regards,
> David
>
> Jean-Paul Pelteret schrieb am Dienstag, 29. Dezember 2020 um 20:06:17
> UTC+1:
>
>> Hi David,
>>
>> So as I'm sure that you know, if you want to assemble the linear system
>> for quasi-static non-linear (hyper)elasticity with a referential DoFHandler
>> (using no Eulerian mapping to transform the shape functions to the spatial
>> setting) AND without doing a similar transformation by hand (as step-44
>> does) then you are pretty much limited to using one of two
>> parameterisations. Either you have to choose the (two-point) Kirchoff
>> stress (P) and deformation gradient (F) as the work conjugate pairs, or the
>> (fully-referential) Piola-Kirchhoff stress (S) and the Green-Lagrange
>> strain tensor (E = 0.5(C-I) with C being the right Cauchy-Green tensor).
>>
>> Thinking about just the residual term: to the best of my knowledge (and
>> from what I recall from the work that we did here
>> <https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6336>), the
>> variation dC = symm(F^{T}.dF) = symm(F^{T} d(Grad(u))) -- with u
>> representing the displacement -- cannot be framed a way that the
>> matrix-free framework supports: the test function cannot be pre-multiplied
>> by a field variable. One would have to somehow rephrase the whole problem
>> such that you end up with just the test function d(Grad(u)) as a
>> pre-multiplication to your stress variable, and in the end you come up with
>> the exactly the other parameterisation because P = F.S and dF = d(I +
>> Grad(u)) == d(Grad(u)). So you may as well start off by parameterising the
>> problem in terms of F and P, and naturally you have to adjust the
>> linearisation as well. The linearisation of the two-point tensor P contains
>> both the material and geometric tangent (you can identify each by
>> linearising the decomposition P = F.S to get these two quantities). I get
>> the sense from what I see in the code that you've shared you're trying
>> to accomplish exactly this (although I'm not sure its quite right, since
>> tau = P.F^{T} --> P = tau.F^{-T}, and you seem to have some extra
>> manipulation involving what you call symm_grad_Nx -- this looks a bit
>> suspect to me).
>>
>> I have some code to share that will do the transformation from the stress
>> and tangent quantities already used in step-44 to those that you need for
>> this F-P parameterisation. You can put these in the PointHistory class
>>
>> template <int dim>
>> class PointHistory
>> { ...
>> // Fully referential configuration
>> SymmetricTensor<2, dim>
>> get_S() const
>> {
>> return Physics::Transformations::Contravariant::pull_back(get_tau(), F);
>> }
>> SymmetricTensor<4, dim>
>> get_H() const
>> {
>> return Physics::Transformations::Contravariant::pull_back(get_Jc(), F);
>> }
>> // Mixed / two-point configuration
>> Tensor<2, dim>
>> get_P() const
>> {
>> return get_F() * Tensor<2, dim>(get_S());
>> }
>> Tensor<4, dim>
>> get_HH() const
>> {
>> // Reference: Simo1984a eq 4.4 // Simo, J. C. and Marsden, J. // On the
>> rotated stress tensor and the material version of the {Doyle-Ericksen}
>> formula // DOI: 10.1007/BF00281556
>> const SymmetricTensor<2, dim> S = get_S();
>> const SymmetricTensor<4, dim> H = get_H();
>> // Push forward index 2 of material stress contribution
>> // This index operation relies on the symmetry of the
>> // last two indices, so therefore (in case it makes a
>> // difference) we transform it first.
>> Tensor<4, dim> tmp1;
>> for (unsigned int I = 0; I < dim; ++I)
>> for (unsigned int J = 0; J < dim; ++J)
>> for (unsigned int k = 0; k < dim; ++k)
>> for (unsigned int L = 0; L < dim; ++L)
>> for (unsigned int K = 0; K < dim; ++K)
>> tmp1[I][J][k][L] += get_F()[k][K] * H[I][J][K][L];
>> Tensor<4, dim> HH_mixed;
>> const Tensor<2, dim> I_ns = Physics::Elasticity::StandardTensors<dim>::I;
>> for (unsigned int i = 0; i < dim; ++i)
>> for (unsigned int J = 0; J < dim; ++J)
>> for (unsigned int k = 0; k < dim; ++k)
>> for (unsigned int L = 0; L < dim; ++L)
>> {
>> // Push forward index 0 of material stress contribution
>> for (unsigned int I = 0; I < dim; ++I)
>> HH_mixed[i][J][k][L] += get_F()[i][I] * tmp1[I][J][k][L];
>> // Add geometric stress contribution
>> HH_mixed[i][J][k][L] += I_ns[i][k] * S[J][L];
>> }
>> return HH_mixed;
>> }
>> }
>> Naturally, the above could be simplified a bit for this fixed
>> parameterisation.
>>
>> With this, the (matrix-based) assembly would looks something like this
>> for the RHS
>> // Excuse the messiness -- coding in an email client is not a good idea!
>> for (unsigned int q_point = 0; q_point < this->n_q_points; ++q_point) {
>> const
>> Tensor<2,dim> P = lqph[q_point]->get_P();
>> for (unsigned int i = 0; i < this->dofs_per_cell; ++i) {
>> const unsigned int i_group = this->fe.system_to_base_index(i).first.first;
>> const Tensor<2,dim> dF_I = fe_values_ref[this->u_fe].gradient(i, q_point);
>>
>> if (i_group == this->u_dof)
>> data.cell_rhs(i) -= double_contract(dF_I, P) * JxW;
>> ... // See the rest of step-44
>>
>> and the salient part of the tangent matrix assembly would be something
>> like
>> for (unsigned int q_point = 0; q_point < this->n_q_points; ++q_point) { //
>> Linearisation of P with respect to F; // contains both material and
>> geometric tangent contributions const Tensor<4,dim> HH =
>> lqph[q_point]->get_HH();
>>
>> for (unsigned int i = 0; i < this->dofs_per_cell; ++i) {
>> const unsigned int i_group = this->fe.system_to_base_index(i).first.first
>> ; const Tensor<2,dim> dF_I = fe_values_ref[this->u_fe].gradient(i,
>> q_point);
>> for (unsigned int j = 0; j <= i; ++j)
>> {
>> const unsigned int j_group = this->fe.system_to_base_index(j).first.first
>> ;
>> const Tensor<2,dim> dF_J = fe_values_ref[this->u_fe].gradient(j, q_point)
>> ;
>>
>> if ((i_group == j_group) && (i_group == this->u_dof))
>> {
>> cell_matrix(i, j) += scalar_product(dF_I, HH, dF_J) * JxW; } ... // See
>> the rest of step-44
>>
>> So IIRC the call to phi.submit_gradient() that is bound to the
>> undeformed/non-mapped DoFHandler is effectively the same as testing against
>> what I've called dF_I in the above code. I think that the transformations
>> from tau->P and Jc->HH should make it simple enough to adjust what remains.
>>
>> I hope that this helps you implement what you're trying to do.
>> Jean-Paul
>>
>>
>>
>> On 29.12.20 15:59, 'David' via deal.II User Group wrote:
>>
>> Maybe as an edit: what I currently do looks the following way:
>>
>> ```
>> // Get gradient in reference frame
>> const Tensor<2, dim, VectorizedArrayType> grad_u = phi.get_gradient(
>> q_point);
>> // Compute deformation gradient
>> const Tensor<2, dim, VectorizedArrayType> F = Physics::Elasticity::
>> Kinematics::F(grad_u);
>> const SymmetricTensor<2, dim, VectorizedArrayType> b = Physics::
>> Elasticity::Kinematics::b(F);
>>
>> const VectorizedArrayType det_F = determinant(F);
>> // Invert F
>> const Tensor<2, dim, VectorizedArrayType> F_inv = invert(F);
>> const SymmetricTensor<2, dim, VectorizedArrayType> b_bar = Physics::
>> Elasticity::Kinematics::b( Physics::Elasticity::Kinematics::F_iso(F));
>> // Get tau from material description
>> SymmetricTensor<2, dim, VectorizedArrayType> tau; material->get_tau(tau,
>> det_F, b_bar, b);
>>
>> // Gradient itslef is included in the integrate call, apply push forward
>> with F^{-1}
>> const Tensor<2, dim, VectorizedArrayType> symm_grad_Nx = symmetrize(F_inv
>> );
>> // Compute the result
>> const Tensor<2, dim, VectorizedArrayType> result = symm_grad_Nx * Tensor<
>> 2, dim, VectorizedArrayType>(tau);
>> // Apply an additional push forward with F^{-T}
>> phi.submit_gradient(-result * transpose(F_inv), q_point); // end loop
>> over quadrature points
>>
>> // For each element
>> phi.integrate(false, true);
>> ```
>>
>> It works, but I don't get the quadratic convergence order of the Newton
>> scheme anymore. Hence, I assume this is not sufficient.
>> David schrieb am Dienstag, 29. Dezember 2020 um 12:41:35 UTC+1:
>>
>>>
>>> Hi all,
>>>
>>> I'm currently trying to implement a vectorized variant of the residual
>>> assembly as it is done in step-44/one of the corresponding code-gallery
>>> examples using FEEvaluation objects in combination with matrix-free.
>>>
>>> I was not able to find an appropriate solution for the given code line
>>> <https://github.com/dealii/dealii/blob/57ca3f894d1bb31f6a6dd448864f4386dc6692ad/examples/step-44/step-44.cc#L2040>and
>>>
>>> the subsequent application for the assembly process: the major problem is
>>> that the shape function gradients are defined with regard to the current
>>> configuration, whereas my matrix-free object holds a mapping, which refers
>>> to the initial configuration. The reference configuration mapping of my
>>> matrix-free object is desired as the final integration is actually
>>> performed in the reference frame.
>>> A valid option is probably (I have not tried it) to use a matrix-free
>>> object with a different mapping (MappingQEulerian) which directly evaluates
>>> the shape gradients in the deformed configuration and use another
>>> referential matrix-free object for the integration part.
>>> The main reason to avoid this approach is that it would require a
>>> reinitialization of the 'deformed' matrix-free object in each nonlinear
>>> step and the mapping needs to be recomputed in each reinizialization. On
>>> the other hand, the expensive reinitialization is not really required as
>>> the mapping between the referential and the current configuration is known
>>> due to the (inverse) deformation gradient.
>>> Therefore, I'm looking for an opportunity to access the shape gradients
>>> and perform a push forward similarly to the approach in step-44 in order to
>>> evaluate the desired quantities.
>>>
>>> Until now I was not able to find an obvious solution for this
>>> computations using the FEEvaluation class. Maybe anyone else has a better
>>> idea for the described problem.
>>>
>>> Thanks in advance,
>>> David
>>>
>> --
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> ---
>> You received this message because you are subscribed to the Google Groups
>> "deal.II User Group" group.
>> To unsubscribe from this group and stop receiving emails from it, send an
>> email to dealii+un...@googlegroups.com.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/7b3c1488-cf75-4526-8442-d1045dc73398n%40googlegroups.com
>>
>> <https://groups.google.com/d/msgid/dealii/7b3c1488-cf75-4526-8442-d1045dc73398n%40googlegroups.com?utm_medium=email&utm_source=footer>
>> .
>>
>>
>> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com.
>
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/1062eaa4-0e32-43f8-8fc4-70dfe2638a81n%40googlegroups.com
>
> <https://groups.google.com/d/msgid/dealii/1062eaa4-0e32-43f8-8fc4-70dfe2638a81n%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/e0ccf5a8-636d-4f13-95aa-16810da97f0dn%40googlegroups.com.
