Problem with setting initial conditions and history variable

Silvio · November 27, 2024, 3:58pm

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, tm is above 293 K, as expected.
- On odd timesteps, tm inexplicably drops to 0 K, which shouldn’t happen as the temperature should remain above the initial value.
Suspected Issue: I believe the problem is related to how I’m saving and updating my history variable (tm_n). Specifically:
- If I replace tm with a fixed value like tm_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 a nan error caused by a division by zero (e.g., an exponent function with 1/0).

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!

Topic		Replies	Views
Timestep and Convergence Error Behavior - Need Help General	1	42	August 25, 2024
Problem with convergence General	0	40	July 18, 2024
Issue for assigning initial condition in Fenicsx General	8	690	November 17, 2022
How to set the initial conditions in FEniCSx General	2	1285	December 7, 2022
Dolfinx: Time-dependent problem equilibrating instantly dolfinx	1	554	April 6, 2022

Problem with setting initial conditions and history variable

Related topics