GOLF rheology in ParStokes Solver?

Discussion about coding and new developments
Post Reply
kfourteau
Posts: 15
Joined: 25 Jun 2020, 16:45
Antispam: Yes

GOLF rheology in ParStokes Solver?

Post by kfourteau »

Hello Everyone,

In Brief: I would like to use the anisotropic GOLF material law with the ParStokes solver for large simulations. It requires some modification of the ParStokes.F90 file, in order to use non scalar viscosities, but I do not know how to concretly do it.

My Problem:
I have tried to perform some simulations of the deformation of viscoplastic porous medium with different non-linear rheologies. My meshes have complicated geometries and typically have about 10 million elements. I am performing steady state simulations, with the Picard method for the non-linear iterations, and a mix of Dirichlet and free surfaces boundary conditions.

For isotropic rheologies I am using ParStokes with an iterative method and Glen's law as as the material model. I am really satisfied with this solver, the convergence is relatively fast and the results quite stable. As a comparison, I have also tried the standard Stokes solver. With the standard Stokes solver I have some convergence issues, and the results are not so satisfactory. Thus, for non-linear isotropic materials, the ParStokes solver, and its block preconditioning, perfectly suits my needs.
For my study I'd relly like to also use some anisotropic material. For that I started using the GOLF material law with the AIFlow solver. Unfortunately, I run into the same problems as with the standard Stokes solver: convergence is difficult and the results are not so robust.

I think what could fix my problem would be to be able to use the GOLF rheology with the ParStokes Solver (the best of both worlds).

What I Have Done So Far:
I have created a modified version of the ParStokes.F90 which could hopefully use the GOLF rheology. Unfortunately, it is not straight forward to add GOLF into ParStokes, notably because by default ParStokes requires a scalar viscosity (mu in the code), which has no direct equivalent with the anisotropic GOLF law.

Concretly, I have:
- Copied the Subroutine GetMaterialDefs() from AIFlowSolver_nlD2.F90 to the modified ParStokes.F90. This subroutine reads the Material Parameters required for GOLF.
- Copied the block from AIFlowSolver_nlD2.F90 that reads the Fabric variables (I also read the FlowWidth variable, but I do not intent to use it)
- Modified the LocalMatrix() function to accept more variables (Fabric, etc)
- In LocalMatrix(), copied the block from AIFlowSolver_nlD2.F90 that computes the 'non relative viscosity matrix', called C in the code. It is not clear to me what it is exactly, but I guess this is what should replace the scalar viscosity.
- In LocalMatrix(), copied the block from AIFlowSolver_nlD2.F90 that computes a 3x3 block of the Stiffness matrix. I think this 3x3 block corresponds to the viscosity term in the Galerkin formulation, but as I am not sure I have not try to use it in the ParStokes assembly.

I have checked that the names of the added variables to ParStokes are compatible with those already present. Up to this point I have thus loaded and computed all necessary new variables for GOLF in the modified ParStokes solver, but I have not modified the assembly process. When I run a simulation, the results are the same as with the normal ParStokes, which suggests me that I did not break anything with these new computed variables.

What Still Needs To Be Done:
My problem is now to modify the assembly step, in order to replace the scalar viscosity mu with the appropriate GOLF rheology. I have identified three parts to be changed in the original ParStokes.F90 file:
- L1327/1328 for what I guess is the viscosity term in the sitffness matrix
- L1350 for the velocity block preconditioning
- L1391 for the mass matrix used in the preconditioner P

For L1327/1328 I suspect that injecting the 3x3 block computed earlier with the matrix C of the GOLF rheology might work, but I am really not sure of it. For the two other parts I have no idea.

Do you have an idea on how to properly modify this assembly step to use the GOLF rheology? It is even possible?
I have attached a zip file with a commented copy of my modified ParStokes solver (to be compiled with elmerf90 ParStokes_Modofied.F90 Golf_law.F90 -o ParStokes_Modified.so). There is also a small simulation set-up where the Solver is called.

Thanks a lot !
Regards,
Kévin
Attachments
Example_ModifiedParStokes.zip
(146.54 KiB) Downloaded 92 times
mika
Posts: 207
Joined: 15 Sep 2009, 07:44

Re: GOLF rheology in ParStokes Solver?

Post by mika »

Here an open question is to identify an effective approximation to the Schur complement (the coefficient matrix for the pressure updates in the preconditioning). The current choice is designed for a scalar viscosity and the generalization to anisotropy seems to be a research topic. Other modifications should be more straightforward to implement.

-- Mika
kfourteau
Posts: 15
Joined: 25 Jun 2020, 16:45
Antispam: Yes

Re: GOLF rheology in ParStokes Solver?

Post by kfourteau »

Hi Mika,

Thanks a lot for your answer. I have done a bit a bibliography, and indeed there appears to be no clear and easy way to define a good approximation of the Schur complement with anisotropic viscosities. From what I understand (and this is not my original field of competence so I might be wrong here), the Schur complement can be approximated with:
- A simple mass matrix. This is super simple, but gives a poor approximation.
- A mass matrix weighted by the inverse of the scalar viscosity. This is what is currently done in Elmer, but requires a scalar viscosity.
- A mass matrix weighted by the "isotropic part of the viscosity". I'm not entirely sure what it means, but I guess it might be something like the trace of the fourth order viscosity tensor.
- For anisotropic viscosity, there are some more rigorous but more complicated methods, such as: https://arxiv.org/pdf/2107.00820.pdf . It is certainly quite difficult to integrate in the current version of ParStokes and is probably not worth the effort anyway.

I've modify a bit my ParStokes code to try include these modifications, but without any clear success.

Do you think it would be possible to implement something along the lines of:
- Having the viscosity expressed as a non-isotropic fourth order tensor (something similar to Eq 7 of https://agupubs.onlinelibrary.wiley.com ... 99JB900146 rotated in the general reference frame based on the orientation of the crystal axis).
- Assemble the Stiffness matrix with this fourth order viscosity.
- Skip the diagonalization of the velocity preconditionner for simplicity (hence setting Block Diagonal A to False in the sif).
- Take the mass matrix divided by the trace of the fourth order viscosity as the approximated Schur complement.

I think I can do most of it myself, except for the assembly of the stiffness matrix. I have tried to do the maths to derive the standard Galerkin formulation of the Stokes equation with a tensorial viscosity, but I systematically got lost. Thus I do not know the expression of the stiffness matrix coefficients.

Thanks again,
Kévin
mika
Posts: 207
Joined: 15 Sep 2009, 07:44

Re: GOLF rheology in ParStokes Solver?

Post by mika »

Yes, I also think that using the average of the diagonal entries of the constitutive operator may be the simplest strategy so that the scalar viscosity parameter can be used in the pressure preconditioner. Did you find any article where this strategy is applied, so that one would know beforehand how strong anisotropy could be treated in this way and whether this might be a realistic approach here?

To handle anisotropy, I'd switch to a matrix representation (of the material law) similar to the version used often in linearized elasticity. This usually simplifies the implementation.

-- Mika
kfourteau
Posts: 15
Joined: 25 Jun 2020, 16:45
Antispam: Yes

Re: GOLF rheology in ParStokes Solver?

Post by kfourteau »

I have found one paper that uses the isotropic part of the tensorial viscosity to compute the Schur complement: https://www.sciencedirect.com/science/a ... 1514#b0085 . In the article they say that they use "the lumped mass matrix weighted by the inverse of the isotropic part of the viscosity". They also cite a book to justify this choice, but I do not have access to it so I do not know exactly what is said in it.

I would be quite eager to try this method, to see if it solves my problem. So what I have in mind right now would be to modify ParStokes for (i) non-linear anisotropic viscosity, (ii) Schur complement approximation with isotropic part of viscosity, (iii) Picard method only, and (iv) no diagonalization of velocity preconditionner.

To implement it I would:

- First, add in the Material section the crystalographic orientation of each body with the colatitude and longitude angles.
- Then, modify the ParStokes solver to something like this pseudo-code:

Code: Select all

 !Modify Assembly Step
 DO t=1,Active !Loop over active Elements

  ! Get Material Properties to compute viscosity matrix
  theta_angle(1:n)  = GetReal( Material, 'Colatitude' )
  phi_angle(1:n)  = GetReal( Material, 'Longitude' )
  n_exp(1:n)  = GetReal( Material, 'Non Linear Exponent' )
  Visco_fact(1:n)  = GetReal( Material, 'Viscosity Prefactor' )
  beta(1:n)  = GetReal( Material, 'Anisotropy Parameter' )

  !Enter LocalMatrix() Subroutine    
    DO t=1,IP % n
       !mu_tensor is the fourth order linearized viscosity tensor in the general reference frame
       !mu_tensor0 is an anisotropic viscosity tensor in the grain frame
       !R is the rotation matrix based on the crystalographic angles
       IF (SkipPowerLaw) THEN
         mu_tensor = Visco_fact * Isotropic_tensor             !Skip non-linear law and just put a constant isotropic viscosity (identity in Voigt notation)
         mu_trace = Visco_fact
       ELSE
         mu_tensor = TENSOR_ROTATION(mu_tensor0, R)            ! Rotate Viscosity Matrix in general reference frame
         mu_tensor = mu_tensor * SR_inv**((1/nexp -1)/2)       ! Multiply by Strain Rate second invariant to obtain the effective viscosity of the non-linear rheology
         mu_trace = (mu_tensor(1,1,1,1) + mu_tensor(2,2,2,2) + & 
                     mu_tensor(3,3,3,3) + mu_tensor(1,2,1,2) + & 
                     mu_tensor(1,3,1,3) + mu_tensor(2,3,2,3)) / 6      ! Compute the isotropic part of viscosity tensor (trace of tensor in Voigt notation)
       END IF
     
       DO p=1,nd   !Loop over basis
         DO q=1,nd !Loop over basis

         DO i=1,dim
           DO j = 1,dim
            ! WHAT TO PUT HERE FOR THE VISCOSITY TERM OF THE STIFFNESS MATRIX?
           END DO !End Loop over j
       
           IF (q <= n) & 
             A(i,dim+1) = A(i,dim+1) - s * Basis(q) * dBasisdx(p,i)   !Discrete Divergence Operator
           IF (p <= n) &   
             A(dim+1,i) = A(dim+1,i) - s * dBasisdx(q,i) * Basis(p)   !Discrete Gradient Operator 
         END DO !End Loop over i
         
       IF (DiagonalA) CALL Fatal( 'Modified ParStokes', 'Diagonal A not Implemented' )
       
       END DO !End Loop over q
      !Leave Force vector as is
      END DO !End Loop over p
     
      DO p=1,n
       DO q=1,n       
         Mass(p,q) = Mass(p,q) - s * 1.0d0/mu_trace * Basis(p) * Basis(q)  !New Schur Complement      
       END DO
      END DO       
             
    END DO !End Loop over Integration Points
  !Leave LocalMatrix() Subroutine
 END DO !End Loop over Elements
 
 !Keep the rest of ParStokes as it is
I still do not know how to express properly the viscosity term of the sitffness matrix as a function of the local tensorial viscosity and of the local basis functions. All articles I managed to find do not give the complete expression of the matrix term, and just say it is the discretisation of the viscous operator.
Also in my pseudo-code I used the full form of the viscosity tensor (i.e. a 3x3x3x3 tensor), but perhaps the Voigt notation is more approriate.

Thanks,
Kévin
mika
Posts: 207
Joined: 15 Sep 2009, 07:44

Re: GOLF rheology in ParStokes Solver?

Post by mika »

I think the book mentioned in the paper doesn't really treat anisotropy.

If you can create a code which enables to represent the material law as T + p I = C^{ijkl}D_{kl}, with the righ-hand side expressed by using the Voigt notation (that is, the tensor C is then expected to be a (6x6)-matrix), I could try to help in finalizing the code such that the dot product of the stress and the gradient of a test function is computed correctly in the integration of the stiffnesss matrix. I suppose the code could potentially be an update of the solver in the Git repository.

Best regards,
Mika
kfourteau
Posts: 15
Joined: 25 Jun 2020, 16:45
Antispam: Yes

Re: GOLF rheology in ParStokes Solver?

Post by kfourteau »

I have attached a modified version of the ParStokes Solver that computes an anisotropic viscosity under a (6x6) matrix form.

I tried to follow the original structure of the ParStokes solver. Notably, the computation of the effective viscosity matrix is done in a separate module file, similar to original MaterialModels module.

From what I understood of the Elmer code, the Voigt notation used in Elmer is one where the deviatoric stress is represented as s = [s11, s22, s33, s12, s32, s13]. That is the convention I followed for the computation of the viscosity matrix.

For the computation of the isotropic part of the viscosity, I took the trace as in Eq. 3.10 of Betten (1982): https://www.emis.de/journals/HOA/IJMMS/5/87.pdf , and normalized to obtained the standard viscosity in the isotropic case. I have check that it is indeed invariant under rotation.

I have attached my codes with a small test-case, as well as pdf detailing all the modifications done to the code.
Let me know if you need me to do some more modifications.

Thanks a lot!
Kévin
Attachments
Modifications_ParStokes.pdf
(120.49 KiB) Downloaded 86 times
ParStokes_Aniso.zip
(114.15 KiB) Downloaded 83 times
mika
Posts: 207
Joined: 15 Sep 2009, 07:44

Re: GOLF rheology in ParStokes Solver?

Post by mika »

I committed changes which should be quite close to what are needed here, see

https://github.com/ElmerCSC/elmerfem/co ... b298d39eba

and a small fix

https://github.com/ElmerCSC/elmerfem/co ... 678d1e4c05

I followed the convention that the stress tensor is vectorized by listing its components as [T11,T22,T33,T12,T13,T23], while the symmetric part of the velocity gradient is vectorized as [D11,D22,D33,2*D12,2*D13,2*D23]. This way is convenient as one may then evaluate the energy density as the dot product of these 6-tuples. Now everything is contained within ParStokes.F90, so it also includes the module AnisotropicMaterialModels which should still be finalized so that the necessary material models are available.

If the case of a scalar viscosity is handled in the tensor form, I checked that the computed solution essentially agrees with the solution obtained by handling it without activating the tensor formulation, so I think the tensor formulation satisfies a minimal consistency test.

A good practice in coding is that one avoids having unnecessary copies of the same code in different places if this can be avoided by using USE constructs. In this way I was able to shorten the content of the module AnisotropicMaterialModels to one subroutine.

The remaining task should be to revise the module AnisotropicMaterialModels and further testing. So the variable mu_tensor should be defined such that it acts on the vectorized version of the symmetric part of the velocity gradient to give a stress vector. The stress caused by the pressure is handled separately.

Best regards,
Mika
kfourteau
Posts: 15
Joined: 25 Jun 2020, 16:45
Antispam: Yes

Re: GOLF rheology in ParStokes Solver?

Post by kfourteau »

Thanks a lot for the implementation of anisotropic viscosities! :D

I haven't had the time to test it for now, but I'll do it as soon as possible. I will let you know the tests that I have done and how they went.

Thanks again,
Kévin
Post Reply