**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.**

*In Brief:***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