Hi everyone,
I’m working on a FEniCSx simulation involving a 2-equation system with two variables: temperature (tm) and hydration degree (hyddeg). Here’s a quick overview of my model:
Initial conditions:
-
tm(temperature) starts at 293 K and rises over time. -
hyddeg(hydration degree) starts at 0.1 and progresses to a maximum of 1.0. -
Problem: My simulation produces unexpected results for the temperature (
tm):- On even timesteps,
tmis above 293 K, as expected. - On odd timesteps,
tminexplicably drops to 0 K, which shouldn’t happen as the temperature should remain above the initial value.
- On even timesteps,
-
Suspected Issue: I believe the problem is related to how I’m saving and updating my history variable (
tm_n). Specifically:- If I replace
tmwith a fixed value liketm_0(293.0) in my equations (see the line 184 in the code), the simulation runs without errors but still shows the errors at even/uneven timesteps. - However, when using the primary variable
tm, the solver crashes with ananerror caused by a division by zero (e.g., an exponent function with1/0).
- If I replace
This makes me think I may have made a mistake in setting or updating the initial conditions, particularly for the history variable tm_n. I suspect there might be an issue with how tm_n is being updated during the timestep loop.
Blockquote
from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import VTXWriter, XDMFFile
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 dolfinx.mesh import create_unit_square, locate_entities, meshtags
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,
#------------------------------meshio-------------------------------#
# folder = "model_2eq_trialruns"
#------------------------------quad---------------------------------#
# Geometry
Length = 0.6 # [m]
Height = 0.04 # [m]
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([Length, Height])],
[20, 20], 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)
# checking if it is properly imlemented
tm_n.name = "tm"
hyddeg_n.name = "hyddeg"
with dolfinx.io.XDMFFile(domain.comm, "out_elasticity/displacements.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(tm_n)
file.write_function(hyddeg_n)
#----------------------Boundary Conditions---------------------------#
# We start by identifying the facets contained in each boundary and create a
# custom integration measure ds
boundaries = [(1, lambda x: np.isclose(x[1], 0)), # bottom
(2, lambda x: np.isclose(x[0], Length)), # right
(3, lambda x: np.isclose(x[1], Height)), # top
(4, lambda x: np.isclose(x[0], 0))] # left
# We now loop through all the boundary conditions and create MeshTags
# identifying the facets for each boundary condition.
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(domain, fdim, 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_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
# Debugging boundary condition
# To debug boundary conditions, the easiest thing to do is to visualize the boundary in Paraview
# by writing the MeshTags to file. We can then inspect individual boundaries using the Threshold-filter.
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with XDMFFile(domain.comm, "facet_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(facet_tag, domain.geometry)
#-------------------------------------------------------------------#
list_bc = []
# Setup
ns_0 = 0.75
pc_0 = 8.352 * 1e6
tm_0 = 293.15
# Material Parameters
mat_rho_SR0 = 2.65e3
mat_rho_LR0 = 1e3
mat_rho_GR0 = 1.0
# Gawin saturation
mat_a = 18.6237 * 1e6 # in [Pa]
mat_b = 2.2748 # [-] parameter for saturation function
# thermal capacity of mixture
mat_Cp_s0 = 800.0 # specific heat of solid at t=0
mat_Cp_l0 = 4184.0 # specific heat of water at t=0
mat_Cp_g0 = 1618.0 # specific heat of air at t=0
# Hydration evolution 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]
mat_rhohat_hydr_inf = 69.0 # [kg/(m^3 s)] final stage
#-------------------------------------------------------------------#
# Time Derivatives
dtm_dt = (tm - tm_n) / dt # tm's
dhyddeg_dt = (hyddeg - hyddeg_n) / dt # hyddeg's
# Saturation
sL = (1.0+(pc_0/mat_a)**(mat_b/(mat_b-1.0)))**(-1.0/mat_b)
nF = 1.0 - ns_0
nL = sL * nF
nG = nF - nL
# Mass growth
rho_hat_S = dhyddeg_dt * mat_rhohat_hydr_inf
# thermal capacity of mixture
mat_As = 0.000195
rho_SR = mat_rho_SR0 + mat_As * tm
rhoCp_s = ns_0 * rho_SR * mat_Cp_s0
rhoCp_eff = rhoCp_s + nL * mat_rho_LR0 * mat_Cp_l0 + nG * mat_rho_GR0 * mat_Cp_g0
# Thermal conductivity
lambda_eff = 1.67 * ( 1.0 + 0.0005 * (tm - tm_0) )
# Hydration evolution equation -> ne radi razliku
A_Gamma = mat_A_1 * ((mat_A_2 / mat_kappa_inf) + (mat_kappa_inf * hyddeg)) * (1.0 - hyddeg) * ufl.exp(-mat_eta_bar * hyddeg)
# -> initial condition for the first time step is off (!)
C_todo = ufl.exp(-mat_Ea_R/(tm)) # with "tm" as primary does NOT work; switch to tm_0 = 293.15 [K] and then it converges but still gives the error at uneven timesteps for tm
#-------------------------------------------------------------------#
# Integration measures
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, 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 * C_todo)) * del_hyddeg) * dx
weak_form = Ener_Mix + Hyd_M
#-------------------------------------------------------------------#
# Solver
# Initialise non-linear problem
problem = NonlinearProblem(weak_form, up, list_bc)
# Initialise newton solver
solver = NewtonSolver(domain.comm, problem)
# Set Newton solver options
# solver.damping = 0.8 # Apply damping to Newton updates
solver.atol = 1e-3
solver.rtol = 1e-4
solver.convergence_criterion = "residual" # "incremental" / "residual"
solver.max_it = 50
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)
tm_sol.name = "temperature"
hyddeg_sol.name = "hyddeg"
file1 = dolfinx.io.VTXWriter(domain.comm, "output1.bp", tm_sol, "bp4")
file2 = dolfinx.io.VTXWriter(domain.comm, "output2.bp", hyddeg_sol, "bp4")
#---------------------------Solve problem-------------------------------------#
duration_solve = 0.0
# print
num_steps = int(num_steps)
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)
file1.write(t)
file2.write(t)
# Update time
t += dt
file1.close()
file2.close()
# outfile.close()
If anyone has encountered a similar issue or has insights on how to debug or fix this, I’d greatly appreciate your help! I’ve attached my code for reference. Let me know if additional details are needed.
Thanks in advance!