Hello,
I’m a new user of FEniCSx and I’m trying to implement a poroelastic model written by mfront 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.05m) because my sample is also only several centimeters in size.
I found the simulation became very slow (dt<0.01s), 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).
I have attached my MFront and FEniCSx code for review. I would appreciate any advice on how to improve the file to accelerate the calculation.
- mfront code :
@DSL ImplicitGenericBehaviour;
@Author junfeng;
@Behaviour PoroElasticity_Junfeng_Implicit;
@Algorithm NewtonRaphson_NumericalJacobian ;
@Theta 1. ;
@Epsilon 1e-14 ;
@Description{
An isotropic poroelasticity model
}
@Flux StressStensor sig;
sig.setGlossaryName("Stress");
@Gradient StrainStensor eto;
eto.setGlossaryName("Strain");
@StateVariable StrainStensor eel;
eel.setGlossaryName("ElasticStrain");
@ExternalStateVariable stress pf;
pf.setEntryName("PorePressure");
@Gradient TVector Dpf;
Dpf.setEntryName("PorePressureGradient");
@Flux TVector w;
w.setEntryName("FluidFlux");
// Elasticity
@MaterialProperty stress young;
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu;
nu.setGlossaryName("PoissonRatio");
@LocalVariable stress lambda;
@LocalVariable stress mu;
// Hydraulic
@Parameter real Ks = 21.7e9; // Solid modulus
Ks.setEntryName("UnjacketedModulus"); // example : Pa
/* Belmokhtar M, Delage P, Ghabezloo S, Tang A-M, Menaceur H, Conil N. Poroelasticity of the Callovo–Oxfordian Claystone. Rock Mech Rock Eng 2017;50:871–89. https://doi.org/10.1007/s00603-016-1137-3.*/
@MaterialProperty real kF ;
kF.setEntryName("FluidPermeability");
@MaterialProperty real biot ;
biot.setEntryName("biot") ;
@MaterialProperty real Phi0 ;
Phi0.setEntryName("InitialPorosity") ;
@AuxiliaryStateVariable real Phi ;
Phi.setGlossaryName("Porosity") ;
@AuxiliaryStateVariable real mf;
mf.setEntryName("FluidMass");
@LocalVariable real rhof;
rhof.setEntryName("FluidDensity");
@LocalVariable real Kf ;
Kf.setEntryName("FluidBulkModulus");
@LocalVariable real muf ;
muf.setEntryName("FluidDynamicViscosity");
@LocalVariable real NN ;
NN.setEntryName("BiotN");
@LocalVariable real MM ;
MM.setEntryName("BiotM");
@TangentOperatorBlocks{dsig_ddeto, dsig_ddpf, dmf_ddeto, dmf_ddpf, dw_ddpf , dw_ddDpf};
@InitLocalVariables{
if(Phi == 0 ){
Phi = Phi0;
}
lambda = computeLambda(young, nu);
mu = computeMu(young, nu);
}
@Integrator{
sig = lambda * trace(eel + theta * deel) * I₂ + 2* mu * (eel + theta * deel) - biot * (pf + theta * dpf) * I₂ ;
feel = deel - deto ;
}
@UpdateAuxiliaryStateVariables{
Kf = 2.3e+9 ;
muf = 1.e-3 ;
rhof = 1000.*(1+(pf+dpf)/Kf) ;
MM = 1./ ( 2*(1-biot)* ((1-nu)*biot/young - nu*biot/(young)) + (1-biot)/young*(biot-2*nu*biot) + Phi*(1./Kf-1./Ks) );
NN = 1./( 1/MM - Phi/Kf );
// NN = 1./((biot - Phi)/Ks) ;
Phi += biot * trace(deto) + 1/NN*dpf ;
mf = rhof * Phi ;
w = - kF * rhof / muf * (Dpf + dDpf) ;
}
@TangentOperator{
dsig_ddeto = lambda * (I₂ ⊗ I₂) + 2 * mu * I₄;
dsig_ddpf = - biot * I₂ ;
dmf_ddeto = biot * rhof * I₂ ;
dmf_ddpf = rhof*Phi/Kf + rhof/NN ;
dw_ddDpf = - kF * tmatrix<N, N, real>::Id() * rhof / muf ;
dw_ddpf = - kF * rhof / muf / Kf * (Dpf+dDpf);
}
- FEniCSx code
from dolfinx import mesh, fem, default_scalar_type, io
import ufl
import basix
from dolfinx.cpp.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pandas as pd
import time
from dolfinx_materials.quadrature_map import QuadratureMap
from dolfinx_materials.solvers import NonlinearMaterialProblem
from dolfinx_materials.material.mfront import MFrontMaterial
from dolfinx_materials.utils import symmetric_tensor_to_vector
import os
from tqdm import tqdm
# import sys
current_path = os.getcwd()
start_time =time.time()
Lx = 1
Ly = 1
Nx = 20
Ny = 20
domain = mesh.create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (Lx, Ly)], [Nx, Ny], mesh.CellType.quadrilateral)
x = ufl.SpatialCoordinate(domain)
def axi_grad(u):
return ufl.as_vector([u[0].dx(0), u[1].dx(1), u[0]/x[0], (u[1].dx(0)+u[0].dx(1))/2])
mat_prop = {
"YoungModulus": EE,
"PoissonRatio": nu,
"FluidPermeability": kf,
"biot": Biot,
"InitialPorosity": Phi0,
}
material = MFrontMaterial(
os.path.join(current_path, "src/libBehaviour.dylib"),
"PoroElasticity_Junfeng_Implicit",
hypothesis="axisymmetric",
material_properties=mat_prop
)
V_u = basix.ufl.element("Lagrange", domain.basix_cell(), 1, shape=(2,))
V_pf = basix.ufl.element("Lagrange", domain.basix_cell(), 1)
V = fem.functionspace(domain, basix.ufl.mixed_element([V_u, V_pf]))
v = fem.Function(V, name="Unknown")
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], Lx)
def top(x):
return np.isclose(x[1], Ly)
def bottom(x):
return np.isclose(x[1], 0)
# Initial condition
pf_ini = 1e6
v.sub(1).interpolate(lambda x: np.full((x.shape[1],), pf_ini))
# Boundary condition
top_facets = mesh.locate_entities_boundary(domain, 1, top)
bottom_facets = mesh.locate_entities_boundary(domain, 1, bottom)
left_facets = mesh.locate_entities_boundary(domain, 1, left)
right_facets = mesh.locate_entities_boundary(domain, 1, right)
# Displacement
bottom_dofs_y = fem.locate_dofs_topological(V.sub(0).sub(1), 1, bottom_facets)
bc_u_bottom = fem.dirichletbc(default_scalar_type(0), bottom_dofs_y, V.sub(0).sub(1))
# Pore pressure
bottom_dofs_pf = fem.locate_dofs_topological(V.sub(1), 1, bottom_facets)
bc_pf_bottom = fem.dirichletbc(pf_ini, bottom_dofs_pf, V.sub(1))
bcs = [bc_u_bottom, bc_pf_bottom]
# Loadings
boundaries = [(1, top), (2, right)]
facet_indices, facet_markers = [], []
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, 1, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags = mesh.meshtags(domain, 1, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
f1 = fem.Constant(domain, np.array([0.0, 0.0], dtype=np.float64)) # Force on the top
f2 = fem.Constant(domain, np.array([0.0, 0.0], dtype=np.float64)) # Force on the bottom
# Monotonic loading
F_ini_value = np.array([-1.e6], dtype=np.float64)
F_fin_value = np.array([-5.e6], dtype=np.float64)
F_load_time_value = np.array([1000], dtype=np.float64)
F_ini = fem.Constant(domain, F_ini_value)
F_fin = fem.Constant(domain, F_fin_value)
F_load_time = fem.Constant(domain, F_load_time_value)
v.sub(1).interpolate(lambda x: np.full((x.shape[1],), F_ini.value[0]))
F_rate_value = (F_fin_value - F_ini_value) / F_load_time_value
F_rate = fem.Constant(domain, F_rate_value)
deg_quad = 2
qmap = QuadratureMap(domain, deg_quad, material)
(u, pf) = ufl.split(v)
qmap.register_gradient("Strain", axi_grad(u))
sig = qmap.fluxes["Stress"]
qmap.register_external_state_variable("PorePressure", pf)
qmap.register_gradient("PorePressureGradient", ufl.grad(pf))
w = qmap.fluxes["FluidFlux"]
mf = qmap.internal_state_variables["FluidMass"]
mf_n = mf.copy()
dt = fem.Constant(domain, 0.0)
v_ = ufl.TestFunction(V)
(u_, pf_) = ufl.split(v_)
dv = ufl.TrialFunction(V)
F = (
( (mf-mf_n)/dt*pf_ - ufl.dot(w, ufl.grad(pf_)) ) *2*np.pi*x[0]*qmap.dx
+ ufl.dot(sig, axi_grad(u_))* 2*np.pi*x[0]*qmap.dx
- ufl.dot(f1, u_)*2*np.pi*x[0]*ds(1)
- ufl.dot(f2, u_)*2*np.pi*x[0]*ds(2) )
Jac = qmap.derivative(F, v, dv)
problem = NonlinearMaterialProblem(qmap, F, Jac, v, bcs)
qmap.update()
mf.vector.copy(mf_n.vector)
newton = NewtonSolver(MPI.COMM_WORLD)
newton.rtol = 1e-4
newton.atol = 1e-4
newton.convergence_criterion = "incremental"
newton.report = True
newton.max_it = 50
max_time = 6000
initial_dt = 10
multiple_scale = 2
devisor = 2
min_dt = 20
max_dt = 100
current_time = 0
dt_value = initial_dt
while current_time < max_time:
# back up of current results
backup_fluxes = {name: flux.vector.copy() for (name, flux) in qmap.fluxes.items() }
backup_isv = {name: flux.vector.copy() for (name, flux) in qmap.internal_state_variables.items() }
backup_v = v.vector.copy()
# print("current time = ", current_time, ", dt = ",dt_value)
trial_dt = min(dt_value * multiple_scale, max_time - current_time, max_dt)
dt.value = trial_dt
if(current_time + trial_dt < F_load_time.value[0]):
f1.value[1] = F_ini.value[0] + F_rate.value[0] * (current_time + trial_dt)
f2.value[0] = F_ini.value[0] + F_rate.value[0] * (current_time + trial_dt)
else:
f1.value[1] = F_fin.value[0]
f2.value[0] = F_fin.value[0]
try:
converged, it = problem.solve(newton, print_solution=False)
except RuntimeError as e:
while trial_dt > min_dt:
print("current time = ", current_time, f"Step failed with trial dt: {trial_dt}, reducing dt and retrying...")
# reassign current values to flux isv and v(the problem name)
backup_v.copy(v.vector)
for (name, flux) in qmap.fluxes.items():
backup_fluxes[name].copy(flux.vector)
for (name, flux) in qmap.internal_state_variables.items():
backup_isv[name].copy(flux.vector)
trial_dt = max(trial_dt / devisor, min_dt)
dt.value = np.array([trial_dt], dtype=np.float64)
if current_time + trial_dt < F_load_time.value[0]:
f1.value[1] = F_ini.value[0] + F_rate.value[0] * (current_time + trial_dt)
f2.value[0] = F_ini.value[0] + F_rate.value[0] * (current_time + trial_dt)
else:
f1.value[1] = F_fin.value[0]
f2.value[0] = F_fin.value[0]
try:
converged, it = problem.solve(newton, print_solution=False)
if converged:
break
except RuntimeError as e:
if trial_dt == min_dt:
print("Step failed at minimum dt. Exiting...")
raise RuntimeError("Solver failed even with reduced time step")
continue
mf.vector.copy(mf_n.vector)
current_time += trial_dt
dt_value = trial_dt
end_time = time.time()
print(f"Elipsed time = {(end_time-start_time):.1f} s")
Thank you in advance !
Junfeng REN