Hello everyone,
I’m new to FEniCSx, transitioning from PANDAS for finite element analysis. My material model, which works correctly in FreeFEM, involves two equations: an energy equation (temperature “tm”) and an evolution equation driving a temperature spike (hydration degree “hyddeg”). The setup uses a sealed rectangular specimen with zero environmental fluxes.
However, in FEniCSx, I’m seeing a strange issue where the error decreases as the timestep size increases, which doesn’t make sense. Has anyone else encountered this or have any advice?
Thanks!
Here you can see my whole simplified model:
from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import VTXWriter
from dolfinx import nls
import numpy as np
import ufl
import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc #solvers
from pathlib import Path
from dolfinx.io import gmshio
import meshio # gmsh,
#------------------------------mesh-------------------------------#
Length = 0.6
Width = 0.04
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([Length, Width])],
[30, 2], cell_type=mesh.CellType.quadrilateral)
# Dimensions
fdim = domain.topology.dim - 1
dim = domain.topology.dim
#-------------------------------------------------------------------#
# Function Spaces and Mixed Space
tm_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # tm: 1. primary variable, linear element
hyddeg_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # hyddeg_el: 2. primary variable, linear element
element = ufl.MixedElement([tm_el, hyddeg_el]) # Mixed Element
V = fem.FunctionSpace(domain, element) # Mixed Function space for all dofs
up = fem.Function(V) # total dofs in the problem/ func for all dofs
print("up dofs:", len(up.x.array))
tm_sp = fem.FunctionSpace(domain, tm_el) # tm Function space
hyddeg_sp = fem.FunctionSpace(domain, hyddeg_el) # hyddeg Function space
Vtm, up_to_tm = V.sub(0).collapse() # dofs in u func subspace
tm_n = fem.Function(Vtm) # dofs in u func
Vhyd, up_to_hyddeg = V.sub(1).collapse() # dofs in u func subspace
hyddeg_n = fem.Function(Vhyd) # dofs in u func
tm, hyddeg = ufl.split(up) # (trail function) give two subfunctions for two elements/dofs from mixed space. Otherwise decoupled. for nonlinear don't use trial fn
(del_tm, del_hyddeg) = ufl.TestFunctions(V)
#-------------------------------------------------------------------#
# Define temporal parameters
t = 0.0 # Start time
T = 172800.0 # Final time (2days)
num_steps = 1000
dt = T / num_steps # time step size
#---------------------------Initial Conditions---------------------------#
# Initialize history values
up.x.array[:] = 0.0
tm_n.x.array[:] = 0.0
hyddeg_n.x.array[:] = 0.0
def tm_init(x):
values = np.zeros((1, x.shape[1]))
values[0] = 293.15 # [K]
return values
def hyddeg_init(x):
values = np.zeros((1, x.shape[1]))
values[0] = 0.1 # [-]
return values
tm_n.interpolate(tm_init)
hyddeg_n.interpolate(hyddeg_init)
#-------------------------------------------------------------------#
# Material Parameters
# thermal capacity of mixture
rhoCp_eff = 2.35e6 # [J/m^3K]
# Hydration evolution equation parameters
mat_Ea_R = 5000.0 # [K]
mat_A_1 = 0.95e4 # [-]
mat_A_2 = 0.5e-5 # [-]
mat_kappa_inf = 0.72 # [-]
mat_eta_bar = 5.3 # [-]
# Mass growth parameters
mat_h_hydr = 2.5e6 # [J/kg] heat of hydration process
mat_rhohat_hydr_inf = 69.0 # [kg/(m^3 s)] final stage of hydration
#-------------------------------------------------------------------#
# Time Derivatives
dtm_dt = (tm - tm_n) / dt # tm's
dhyddeg_dt = (hyddeg - hyddeg_n) / dt # hyddeg's
# Mass growth
rho_hat_S = dhyddeg_dt * mat_rhohat_hydr_inf
# Thermal conductivity
lambda_eff = 1.67 * ( 1.0 + 0.0005 * (tm - 293.15) )
# Hydration evolution equation
A_Gamma = mat_A_1 * ((mat_A_2 / mat_kappa_inf) + (mat_kappa_inf * hyddeg)) * (1.0 - hyddeg) * ufl.exp(-mat_eta_bar * hyddeg)
#-------------------------------------------------------------------#
# Integration measures
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, metadata=metadata) #line/area element
dx = ufl.Measure("dx", domain=domain, metadata=metadata) #volume element
# Weak Form
# Energy balance of mixture
Ener_Mix = (rhoCp_eff * dtm_dt * del_tm )* dx
Ener_Mix += (lambda_eff * ufl.inner(ufl.grad(tm) , ufl.grad(del_tm))) * dx
Ener_Mix += - (rho_hat_S * mat_h_hydr) * del_tm * dx
# Hydration evolution equation
Hyd_M = ((dhyddeg_dt - (A_Gamma * ufl.exp(-mat_Ea_R/tm))) * del_hyddeg) * dx
weak_form = Ener_Mix + Hyd_M
#-------------------------------------------------------------------#
# Solver
# Initialise non-linear problem and newton solver
problem = NonlinearProblem(weak_form, up)
solver = NewtonSolver(domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-3
solver.rtol = 1e-4
solver.convergence_criterion = "residual" #"incremental"
solver.max_it = 10
solver.report = True
# Configure mumps
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
#-----------------------------Solution Settings-----------------------------#
Vtm_sol, up_to_tm_sol = V.sub(0).collapse() # dofs in u func subspace
tm_sol = fem.Function(Vtm_sol)
Vhyd_sol, up_to_hyddeg_sol = V.sub(1).collapse() # dofs in u func subspace
hyddeg_sol = fem.Function(Vhyd_sol)
#---------------------------Solve problem-------------------------------------#
duration_solve = 0.0
for n in range(num_steps):
# Solve newton-steps
log.set_log_level(log.LogLevel.INFO)
duration_solve -= MPI.Wtime()
niter, converged = solver.solve(up)
duration_solve += MPI.Wtime()
PETSc.Sys.Print(
"Phys. Time {:.4f}, Calc. Time {:.4f}, Num. Iter. {}".format(
t, duration_solve, niter
)
)
tm_n.x.array[:] = tm_sol.x.array
hyddeg_n.x.array[:] = hyddeg_sol.x.array
tm_sl = up.sub(0).collapse()
hyddeg_sl = up.sub(1).collapse()
tm_sl.x.scatter_forward()
hyddeg_sl.x.scatter_forward()
tm_sol.interpolate(tm_sl)
hyddeg_sol.interpolate(hyddeg_sl)
# Update time
t += dt