Optimizing MFront Poroelastic Model in FEniCSx for Low Permeability Simulations

[yy74774]

Hello,

I’m a new user of MFront and I’m trying to implement a poroelastic model into FEniCSx to simulate tests I have conducted in the lab. I simulated a mechanical loading process for the axisymmetric case while keeping the bottom edge drained.

Initially, I ran my code with a larger element size ($dz = 0.05m$) and high permeability ($k = 10^{-16}m^2$). I checked the results against analytical solutions, and it seemed to work well.

However, when I adjusted the parameters to the real ones:

• Extremely low permeability ($<10^{-20}m^2$) for the Callovo-Oxfordian claystone,
• Smaller element size ($dz = 0.002m$) because my sample is also only several centimeters in size.

I found the simulation became very slow ($dt < 0.01 s$), and it would take several days to complete a simulation of $t_{total}=6000s$. I’m not sure which part causes the problem, or if it’s inherently supposed to be this slow (which would be very unfortunate).

The parameters I used for validation of the model:
Lx = 1, Ly = 1, Nx = 20, Ny = 20, kf = 1e-16
And the for the real tests
Lx = 0.02, Ly = 0.04, Nx = 10, Ny = 20, kf = 1e-21

I have attached my MFront file and the FEniCSx code for review. I would appreciate any advice on how to improve the MFront file to accelerate the calculation.

Thank you in advance !
Junfeng REN

Poroelasticity.zip

[thelfer]

Hi @yy74774,
I'll be back from vacations next week, I'll try to have a look. Would bé far easier however if you could exhibit a simple test with MTest
Best,
Thomas

[yy74774]

Hi Thomas,

Thanks a lot for your rapid reponse.

I wirte a simple mtest file (in attachement) imposing pore pressure $p_f$ and the gradient of pore pressure $\nabla p_f$, it seems like it works well.
Poroelasticity_mtest.txt
Actually, I tried different permeabilities in FEniCSx. The simulation converged quickly when $k_f>8\times 10^{-18}m^2$, and the timestep decreased suddenly to $dt<0.1s$ when $k_f=6\times 10^{-18}m^2$. I'm not sure if this infomation is useful but I think maybe there are some errors in FEniCSx between these two permeabiliies.

Best,
Junfeng

[thelfer]

Ok, i'll review This next week. Don't hésitante to ping me if I am too slow to answer

[yy74774]

Hello,

I believe the issue lies within the FEniCSx part of the simulation. When I changed the unit from meters to millimeters (resulting in the permeability becoming $10^{-12}mm^2$), the calculation speed improved significantly.
For now, it seems the simulation is working fine with millimeters. Thanks a lot for your valuable time:)

Best,
Junfeng

[thelfer]

@yy74774 At first sight I don't understand your trouble.
Your behaviour is mechanically linear. The other parts seems weakly non linear if not linear, so no convergence issue are to be expected.
What is precisely the trouble ? Is it the linear solver that fails or the Newton method.

BTW, the implementation could be much simplier byt using the DefaultGenericBehaviour DSL, why are you doing a nonlinear resolution for a linear behaviour ?

Moreover, there are mix between gradient and auxiliary state variables. I think that it would be way clearer in this case to use only gradients and dual thermodynamic forces.

[yy74774]

Hello,

Thank you for your feedback on my code. I’d like to clarify the issue further.

I am working on modeling poroelastic behavior (a linear behavior) using MFront and implementing it into FEniCSx. To add plasticity in the future, I am using a non-linear solver in FEniCSx and the ImplicitGenericBehaviour in MFront.

When the permeability is larger than 1e-17 m², the timestep can be as large as 100 seconds, and the convergence requires only 2 iterations with the NewtonSolver, which I think is reasonable.

The issue arises when the permeability k is decreased from 8e-18 to 6e-18 m². In this range, the timestep significantly decreases to 0.04 seconds, and 3 iterations are needed for convergence. This might suggest a threshold between these permeability values that causes the NewtonSolver to struggle.

My solution is to change the unit system to $kg, mm, s$ to avoid potential numerical calculation issues. Now, I use 1e-12 mm² instead of 1e-18 m² (although they are equivalent), and this approach works well.

Lastly, I am uncertain about the mix between gradient and auxiliary state variables. Does this imply that it would be clearer to change the type of variable of mf (fluid mass) in my MFront code or something else ?

I sincerely appreciate your response. As a beginner with MFront and FEniCSx, I have much to learn and understand. If you have any further thoughts, please feel free to share them with me.

Best regards,
Junfeng REN

[thelfer]

Are you using an iterative solver ? Change in units may change conditionning and convergence criterion

[yy74774]

Rebonjour, I confirm that I am using an iterative solver.

[thelfer]

Hi @yy74774,
This is the nonlinear solver which is built on top of a linear solver. My question was specifically about the linear solver :)