Fix possible problem with finite strain #1249
Conversation
|
The two new tests seem to run in 3d but the stored output files were generated in 2d, or the other way around. Can you update them? |
|
In your description above, you state the ODE |
| @@ -31,8 +31,9 @@ namespace aspect | |||
| IntegratedStrain<dim>::initialize_one_particle_property(const Point<dim> &, | |||
| std::vector<double> &data) const | |||
| { | |||
| for (unsigned int i = 0; i < SymmetricTensor<2,dim>::n_independent_components ; ++i) | |||
| data.push_back(0.0); | |||
| const Tensor<2,dim> identity = Tensor<2,dim>(unit_symmetric_tensor<dim>()); | |||
There was a problem hiding this comment.
I don't think you need the extra cast on the rhs, or do you? In either case, you can probably make the variable static.
| // (= time step * (velocity gradient tensor - symmetric_part)) | ||
| // + unit tensor | ||
| const Tensor<2,dim> rotation (dt * (grad_u - deformation_rate) + unit_symmetric_tensor<dim>()); | ||
| const Tensor<2,dim> k1 = grad_u * old_strain * dt; |
There was a problem hiding this comment.
can you add a comment before this block that (i) states the ODE you are integrating, and (ii) just says that you are using an RK4 integrator?
|
@gassmoeller - Thanks for catching this! At this point you know as much or more than I do. I am going through other codes to exactly how others have done it and so far (one code) it appears slightly different than the new approach you've outlined. It would not surprise me if different groups are not calculating exactly the same quantity. I'll come back to this first thing next week after I have a chance to write things out. |
gassmoeller
changed the title
|
Any news? |
|
@bangerth - yes and thanks for the reminder. I'll post comments this weekend. |
|
@gassmoeller - sorry of the long delay in replying. I took a look at a different code I use and their approach involves calculating the new strain by adding the (current_strain_rate_tensor * current_time_ step) to the old strain. With this approach they can match an analytical solution for rotation of a circle undergoing pure or simple shear. So, very similar to the new approach but it is not clear to me why they do not need to solve the ODE? Perhaps I'm just overlooking something very basic here. In your method above, I'm also not clear as to where you are separating |
|
Do not worry, I will take a look later this week. Concerning the analytical solution: Does the code reproduce each tensor component or just the second invariant? |
|
I don't think you can do that, can you? I mean, think about a flow that shears a volume, but then you reverse the flow and it undoes all of the shearing. At the end, the accumulated strain is zero, as is its second invariant. Can you express this kind of thing as an evolution equation on only the second invariant? |
|
@gassmoeller - I almost positive that the test is for the strain axes and not just the invariant. As @bangerth mentioned, I'm not sure if one could test rotation with just the invariant. I just went through a second code and their method of calculating finite strain is slightly different as well. In this case, the authors simply calculate the invariant directly (second_strain_rate_invariant * time_step) rather than calculate the individual strain tensor components and then calculate the invariant. All of this is somewhat disconcerting given that the accumulated strain can play a large role in deformation patterns through strain softening. I've forwarded this discussion to Cedric Thieulot ( @cedrict ) and Susanne Buiter, who each have written their own codes. Susanne Buiter and Susan Ellis' code (SULEC) is the first code I have been referring to. More to follow. |
|
I have confirmed that the new way of calculating finite strain gives correct analytical results for pure shear. Will check tomorrow for simple shear and our old formulation, and if necessary update the cookbook. |
gassmoeller
force-pushed
the
fix_finite_strain
branch
from
0bfd0f1
to
9d16165
Compare
|
I can confirm that the new formulation matches analytical expressions for finite strain in simple and pure shear as described in McKenzie & Jackson 1983, and I have added the according benchmark (including analytical expressions) to this PR. On Monday I will modify and add the compositional field calculation method and the cookbook. |
gassmoeller
force-pushed
the
fix_finite_strain
branch
3 times, most recently
from
6928a2c
to
63c4fe5
Compare
gassmoeller
changed the title
|
I think this PR is ready for a review. The new finite strain implementation passes the benchmarks, the particle and compositional field implementations produce identical results, and I updated the cookbook with a more detailed description about the math, and a new figure. |
| # | ||
| # where u is the velocity, t is time, and a comma represents a derivative | ||
| # in that particular direction. u_x,x in this example is equivalent to | ||
| # -u_y,y so that also: F_11 = e^(u_y,y * t) which is intuitive. |
| # -u_y,y so that also: F_11 = e^(u_y,y * t) which is intuitive. | ||
| # For this benchmark it is important to recognize that the deformation in | ||
| # this model is actually only "pure" pure shear at the origin, since the | ||
| # the deformation field around is a mixture of pure shear and a rotational |
| @@ -0,0 +1,5 @@ | |||
| # This is testdata for the Ascii file tracer generator | |||
| # The last line is outside of the test domain and tests | ||
| # that plugins simply ignore points outside of the domain | ||
| # Columns: x_coordinate y_coordinate | ||
| 0.0 0.0 |
There was a problem hiding this comment.
The comment is puzzling since there is only one line, so to talk about the "last line" requires deeper thought than I can offer right now :-)
| # F_xx = 1.0 | ||
| # F_xy = u_x,y * t | ||
| # F_yx = 0 | ||
| # F_yy = 1.0, |
There was a problem hiding this comment.
Is it correct that the finite strain becomes non-symmetric?
| # The last line is outside of the test domain and tests | ||
| # that plugins simply ignore points outside of the domain | ||
| # Columns: x_coordinate y_coordinate | ||
| 0.5 0.5 |
| @@ -1,70 +1,61 @@ | |||
| # Listing of Parameters | |||
| # --------------------- | |||
| # Model for melting at a mid-ocean ridge | |||
| # Model for global mantle convection with tracking of the finite strain tensor. | |||
There was a problem hiding this comment.
thanks for cleaning this file up a bit!
| The first one, $R e^{n-1} R^T$, rotates $e^{n-1}$, the accumulated strain from all the previous time steps; | ||
| and the second part adds the strain of the current time step. | ||
| Let us denote the accumulated strain at time step $n$ as $\mathbf E^n$. We can calculate its time derivative | ||
| as the tensor product of the current velocity gradient $\mathbf G = \frac{\partial u_i}{\partial x_j}$ and the accumulated strain $\mathbf E^{n-1}$ at the previous time step. |
There was a problem hiding this comment.
Shouldn't this be
There was a problem hiding this comment.
"and the strain
There was a problem hiding this comment.
the term "tensor product" typically means something different; I think you want to say something like "the product of two tensors, namely the current velocity gradient ... and the ...."
| Let us denote the accumulated strain at time step $n$ as $\mathbf E^n$. We can calculate its time derivative | ||
| as the tensor product of the current velocity gradient $\mathbf G = \frac{\partial u_i}{\partial x_j}$ and the accumulated strain $\mathbf E^{n-1}$ at the previous time step. | ||
| While we refer to other studies \cite{McKenzie1983, dahlen1998theoretical, Becker2003} for a derivation of | ||
| this relationship, we can give an intuitive example for the necessity to apply the velocity gradient to the already accumulated strain, instead of simply integrating the velocity gradient over time. Consider a simple one-dimensional ``grain'', where the strain tensor only has one component, the compression in x-direction. If one applies a convergent velocity gradient of 0.5 to it (e.g. a velocity of zero at its left end at $x=0.0$, and a velocity of $-0.5$ at $x=1.0$), integrating the velocity gradient would suggest that the grain reaches a length of zero after two units of time, and would then ``flip'' its orientation, which is clearly non-physical. Instead what happens is its length is cut in half in one unit of time, and continues two be halved per unit of time, which is identical to $\mathbf E^n = 0.5 \cdot \mathbf E^{n-1}$, or $\mathbf E^n = \mathbf G \cdot \mathbf E^{n-1}$. |
There was a problem hiding this comment.
There was a problem hiding this comment.
x-direction ->
There was a problem hiding this comment.
continues two -> continues to
There was a problem hiding this comment.
I stumbled about the fact that this velocity field is not incompressible and tried to come up with a complicated scheme of this being the axis of a 2d flow. But it doesn't have to be incompressible. Maybe just say something like "If one embed this grain into a convergent flow field for a compressible medium where the velocity gradient is..."
| @@ -5983,19 +5979,25 @@ \subsubsection{Tracking finite strain} | |||
|
|
|||
| \begin{figure} | |||
| \centering | |||
| \includegraphics[width=0.5\textwidth]{cookbooks/finite_strain/finite_strain.pdf} | |||
| \includegraphics[width=0.75\textwidth]{cookbooks/finite_strain/finite_strain.pdf} | |||
There was a problem hiding this comment.
Yes, really nicely done picture. How did you do this?
| It would have been feasible to run the same model with tracer particles that track the finite strain, as | ||
| additionally implemented and tested in the simple shear and pure shear benchmarks mentioned in this section. | ||
| Both approaches have specific advantages, and for scientific computations one needs to evaluate the more suitable | ||
| strategy. Compositional fields cover the whole domain, but contain numerical diffusion, effectively reducing the |
There was a problem hiding this comment.
"are affected by numerical diffusion"?
| additionally implemented and tested in the simple shear and pure shear benchmarks mentioned in this section. | ||
| Both approaches have specific advantages, and for scientific computations one needs to evaluate the more suitable | ||
| strategy. Compositional fields cover the whole domain, but contain numerical diffusion, effectively reducing the | ||
| maximum accumulated strain. Tracers only present finite strain values at discrete position, but can be used in |
There was a problem hiding this comment.
provide or represent instead of "present"?
| Both approaches have specific advantages, and for scientific computations one needs to evaluate the more suitable | ||
| strategy. Compositional fields cover the whole domain, but contain numerical diffusion, effectively reducing the | ||
| maximum accumulated strain. Tracers only present finite strain values at discrete position, but can be used in | ||
| restricted domains (and are much faster in this case). If however |
There was a problem hiding this comment.
What do you mean by that? Something like: "but can, if this is desired, used in only a part of the domain"?
| $e_{ij} = \frac{1}{2}\left(\frac{\partial d_i}{\partial x_j} + \frac{\partial d_j}{\partial x_i}\right)$ | ||
| with the displacement $\mathbf d$ the difference between a point's | ||
| position at time $t$ and at time zero.} | ||
| \note{In this section, following \cite{Becker2003, dahlen1998theoretical} we denote the velocity gradient tensor as $\mathbf G$, where $\mathbf G = \nabla \mathbf u^T$, |
There was a problem hiding this comment.
put a comma after the reference
| The first one, $R e^{n-1} R^T$, rotates $e^{n-1}$, the accumulated strain from all the previous time steps; | ||
| and the second part adds the strain of the current time step. | ||
| Let us denote the accumulated strain at time step $n$ as $\mathbf E^n$. We can calculate its time derivative | ||
| as the tensor product of the current velocity gradient $\mathbf G = \frac{\partial u_i}{\partial x_j}$ and the accumulated strain $\mathbf E^{n-1}$ at the previous time step. |
There was a problem hiding this comment.
the term "tensor product" typically means something different; I think you want to say something like "the product of two tensors, namely the current velocity gradient ... and the ...."
| Let us denote the accumulated strain at time step $n$ as $\mathbf E^n$. We can calculate its time derivative | ||
| as the tensor product of the current velocity gradient $\mathbf G = \frac{\partial u_i}{\partial x_j}$ and the accumulated strain $\mathbf E^{n-1}$ at the previous time step. | ||
| While we refer to other studies \cite{McKenzie1983, dahlen1998theoretical, Becker2003} for a derivation of | ||
| this relationship, we can give an intuitive example for the necessity to apply the velocity gradient to the already accumulated strain, instead of simply integrating the velocity gradient over time. Consider a simple one-dimensional ``grain'', where the strain tensor only has one component, the compression in x-direction. If one applies a convergent velocity gradient of 0.5 to it (e.g. a velocity of zero at its left end at $x=0.0$, and a velocity of $-0.5$ at $x=1.0$), integrating the velocity gradient would suggest that the grain reaches a length of zero after two units of time, and would then ``flip'' its orientation, which is clearly non-physical. Instead what happens is its length is cut in half in one unit of time, and continues two be halved per unit of time, which is identical to $\mathbf E^n = 0.5 \cdot \mathbf E^{n-1}$, or $\mathbf E^n = \mathbf G \cdot \mathbf E^{n-1}$. |
There was a problem hiding this comment.
x-direction ->
| Let us denote the accumulated strain at time step $n$ as $\mathbf E^n$. We can calculate its time derivative | ||
| as the tensor product of the current velocity gradient $\mathbf G = \frac{\partial u_i}{\partial x_j}$ and the accumulated strain $\mathbf E^{n-1}$ at the previous time step. | ||
| While we refer to other studies \cite{McKenzie1983, dahlen1998theoretical, Becker2003} for a derivation of | ||
| this relationship, we can give an intuitive example for the necessity to apply the velocity gradient to the already accumulated strain, instead of simply integrating the velocity gradient over time. Consider a simple one-dimensional ``grain'', where the strain tensor only has one component, the compression in x-direction. If one applies a convergent velocity gradient of 0.5 to it (e.g. a velocity of zero at its left end at $x=0.0$, and a velocity of $-0.5$ at $x=1.0$), integrating the velocity gradient would suggest that the grain reaches a length of zero after two units of time, and would then ``flip'' its orientation, which is clearly non-physical. Instead what happens is its length is cut in half in one unit of time, and continues two be halved per unit of time, which is identical to $\mathbf E^n = 0.5 \cdot \mathbf E^{n-1}$, or $\mathbf E^n = \mathbf G \cdot \mathbf E^{n-1}$. |
There was a problem hiding this comment.
continues two -> continues to
| Let us denote the accumulated strain at time step $n$ as $\mathbf E^n$. We can calculate its time derivative | ||
| as the tensor product of the current velocity gradient $\mathbf G = \frac{\partial u_i}{\partial x_j}$ and the accumulated strain $\mathbf E^{n-1}$ at the previous time step. | ||
| While we refer to other studies \cite{McKenzie1983, dahlen1998theoretical, Becker2003} for a derivation of | ||
| this relationship, we can give an intuitive example for the necessity to apply the velocity gradient to the already accumulated strain, instead of simply integrating the velocity gradient over time. Consider a simple one-dimensional ``grain'', where the strain tensor only has one component, the compression in x-direction. If one applies a convergent velocity gradient of 0.5 to it (e.g. a velocity of zero at its left end at $x=0.0$, and a velocity of $-0.5$ at $x=1.0$), integrating the velocity gradient would suggest that the grain reaches a length of zero after two units of time, and would then ``flip'' its orientation, which is clearly non-physical. Instead what happens is its length is cut in half in one unit of time, and continues two be halved per unit of time, which is identical to $\mathbf E^n = 0.5 \cdot \mathbf E^{n-1}$, or $\mathbf E^n = \mathbf G \cdot \mathbf E^{n-1}$. |
There was a problem hiding this comment.
I stumbled about the fact that this velocity field is not incompressible and tried to come up with a complicated scheme of this being the axis of a 2d flow. But it doesn't have to be incompressible. Maybe just say something like "If one embed this grain into a convergent flow field for a compressible medium where the velocity gradient is..."
| @@ -5983,19 +5979,25 @@ \subsubsection{Tracking finite strain} | |||
|
|
|||
| \begin{figure} | |||
| \centering | |||
| \includegraphics[width=0.5\textwidth]{cookbooks/finite_strain/finite_strain.pdf} | |||
| \includegraphics[width=0.75\textwidth]{cookbooks/finite_strain/finite_strain.pdf} | |||
There was a problem hiding this comment.
Yes, really nicely done picture. How did you do this?
| \ref{sec:finite-strain}.} | ||
| \ref{sec:finite-strain} at a time of 67.6~Ma. Top panel: Temperature distribution. Bottom panel: | ||
| Natural strain distribution. Additional black crosses are the scaled eigenvectors of the | ||
| stretching tensor $\mathbf L$, showing the direction of stretching and compression.} |
There was a problem hiding this comment.
So the rotation matrix is effectively unimportant to find the major stress axes?
| \label{fig:finite_strain} | ||
| \end{figure} | ||
|
|
||
| The plugin was tested against analytical solutions for the finite strain tensor in simple and pure shear as described in \url{../benchmark/finite_strain/pure_shear.prm} and \url{../benchmark/finite_strain/simple_shear.prm}. |
There was a problem hiding this comment.
IIRC, you shouldn't need the ../ at the beginning of the paths, right?
| direction (short vertical lines) and stretched in horizontal direction (long horizontal lines). | ||
| The sides of the plume head show the opposite effect. Shear occurs mostly at the | ||
| edges of the plume head, in the plume tail, and in the bottom boundary layer | ||
| (black areas in the natural strain distribution). |
There was a problem hiding this comment.
So is it ok that unstrained material would have a cross with equal length arms?
| It would have been feasible to run the same model with tracer particles that track the finite strain, as | ||
| additionally implemented and tested in the simple shear and pure shear benchmarks mentioned in this section. | ||
| Both approaches have specific advantages, and for scientific computations one needs to evaluate the more suitable | ||
| strategy. Compositional fields cover the whole domain, but contain numerical diffusion, effectively reducing the |
|
I have addressed all comments, and reformulated things a bit (also the example calculation in the manual was intuitive, but wrong, I fixed that). I now consistently use F instead of E (to stay consistent with the papers from where I got the derivation). This also makes more clear that F is the (non-symmetric) deformation gradient tensor, while L is the symmetric stretching tensor, one would associate with the finite strain tensor. To your questions:
It would be nice if someone could read over the new text again, but apart from this this should be ready from my side (I can squash the commits if requested, just leaving them there to make it easier to review). |
|
I'll read through the text again as requested if you take care of #1254 :-) |
|
deal 😄 |
|
@gassmoeller - sorry, missed your note until now. I never heard back regarding the details of the test the other code did, but I believe they tracked the rotation of the finite strain axes for a case of simple shear or pure shear. They certainly could have not gone to large enough strains to see the difference. Thank you for adding the benchmarks! |
| The first one, $R e^{n-1} R^T$, rotates $e^{n-1}$, the accumulated strain from all the previous time steps; | ||
| and the second part adds the strain of the current time step. | ||
| Let us denote the accumulated deformation at time step $n$ as $\mathbf F^n$. We can calculate its time derivative | ||
| as the product of two tensors, namely the current velocity gradient $\mathbf G_{ij} = \frac{\partial u_i}{\partial x_j}$ and the deformation gradient $\mathbf F^{n-1}$ accumulated up to the previous time step, in other words $\frac{\partial \mathbf F}{\partial t} = \mathbf G \mathbf F$, and $\mathbf F_0 = \mathbf I$, with $\mathbf I$ being the identity tensor. |
There was a problem hiding this comment.
I think you will want to call it F^0 since you've had the time step as an upper index above.
| Let us denote the accumulated deformation at time step $n$ as $\mathbf F^n$. We can calculate its time derivative | ||
| as the product of two tensors, namely the current velocity gradient $\mathbf G_{ij} = \frac{\partial u_i}{\partial x_j}$ and the deformation gradient $\mathbf F^{n-1}$ accumulated up to the previous time step, in other words $\frac{\partial \mathbf F}{\partial t} = \mathbf G \mathbf F$, and $\mathbf F_0 = \mathbf I$, with $\mathbf I$ being the identity tensor. | ||
| While we refer to other studies \cite{McKenzie1983, dahlen1998theoretical, Becker2003} for a derivation of | ||
| this relationship, we can give an intuitive example for the necessity to apply the velocity gradient to the already accumulated deformation, instead of simply integrating the velocity gradient over time. Consider a simple one-dimensional ``grain'', where the deformation tensor only has one component, the compression in $x$-direction. If one embeds this grain into a convergent flow field for a compressible medium where the dimensionless velocity gradient is $-0.5$ (e.g. a velocity of zero at its left end at $x=0.0$, and a velocity of $-0.5$ at its right end at $x=1.0$), simply integrating the velocity gradient would suggest that the grain reaches a length of zero after two units of time, and would then ``flip'' its orientation, which is clearly non-physical. |
There was a problem hiding this comment.
Do you want to imply that the grain has size one? Or is it infinitesimal?
gassmoeller
force-pushed
the
fix_finite_strain
branch
from
9a749f5
to
761ed95
Compare
After reading somewhat more into the topic, I suspect the finite strain computation we currently use is intuitive, but might be wrong. I base this on the study of the papers: Becker et. al 2003, GJI ("Comparison of azimuthal seismic anisotropy from surface waves and finite strain from global mantle-circulation models") and McKenzie & Jackson, 1983, EPSL ("The relationship between strain rates, crustal thickening, paleomagnetism, finite strain and fault movements within a deforming zone"). If somebody has suggestions for other helpful studies about finite strain, I would be happy to get suggestions (@naliboff, you probably know more about this than I do). Further down L means the vector gradient of the velocity (a full tensor) and I follow the notation of the McKenzie paper.
What we did so far was:
I + 0.5 (L - L^T), or the rotation an infinitely small hard particle would experience in this flow field).0.5(L+L^T)).The Tensor was initialized with zeroes, which was wrong in itself (it should be initialized with a identity tensor).
However, I think the whole approach above is different from the mentioned papers. They track the deformation tensor
Fby solving the differential equationd/dt (F) = L * F, i.e. they track a full (non-symmetric) tensor, and then afterwards split the full tensor into a symmetric (V) and antisymmetric (R) part asF = V*R.Vthen represents the finite deformation we are interested in. McKenzie mentions in the appendix explicitly that one has to be careful not mixing the non-deformed rotation W with the deformed rotation R, because they are essentially only the same, when nothing deforms (sounds simple, if you write it like this 😄). I have tested both approaches, and can only recover the analytical solutions of the McKenzie paper with their approach. However, if someone has suggestions for good test cases or other analytical solutions, I would be happy to take a look, to feel more confident about what I wrote above.About the code below: It is a bit harder to read, because I now use an RK4 scheme to solve the ODE, but essentially it just says the new strain is the old strain multiplied by the current velocity_gradient tensor times the time step length.