Hello everyone,
I’m new to FEniCSx and have recently transitioned from using PANDAS for finite element analysis. My background has been primarily in implementing my material model in FreeFEM, but I’ve heard positive things about FEniCSx and decided to give it a try.
I have simplified my model to consist of two governing equations: an energy equation (with a temperature (“tm” in code) as first primary variable) and an evolution equation driving a temperature spike (with hydration degree (“hyddeg”) of the material for the second primary variable). The setup involves a rectangular geometry with a completely sealed specimen, thus all environmental fluxes are zero. This setup was functioning well in FreeFEM, however, I am facing challenges with my FEniCSx implementation.
One illogical issue I’ve encountered is regarding the error behavior in relation to the timestep size. Oddly, the error decreases as the timestep size increases, which is counterintuitive and suggests there might be a fundamental error in how I’ve implemented the timestep in FEniCSx.
Attached are my simplified material model and the “case1_coarsemesh.msh” file for the geometry (embedded at the end of my fenicsx code cause of the restrictions in this forum).
Could anyone provide insights or suggestions on what might be causing these issues? Any help would be greatly appreciated as I navigate this new tool.
Thank you for your help!
Silvio
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,
#------------------------------meshio-------------------------------#
#------------------------------quad-------------------------------#
proc = MPI.COMM_WORLD.rank
# Extract data and returns meshio mesh including physical markers
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data.astype(np.int32)]})
return out_mesh
if proc == 0:
# Read in mesh
msh = meshio.read("case1_coarsemesh.msh")
print(msh)
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "quad", prune_z=True)
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mt.xdmf", line_mesh)
MPI.COMM_WORLD.barrier()
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid") # cell tags # mesh info
ct = xdmf.read_meshtags(domain, name="Grid") # numberingsinfo: element info
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
facet_tag = xdmf.read_meshtags(domain, name="Grid") # faces tags
# 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, 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 * ufl.exp(-mat_Ea_R/tm))) * del_hyddeg) * dx
weak_form = Ener_Mix + Hyd_M
#-------------------------------------------------------------------#
# Solver
# Initialise non-linear problem
problem = NonlinearProblem(weak_form, up)
# Initialise newton solver
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
# --------------------------------------------------------------------------------------------------------------------
# ------------------ This next part is my "case1_coarsemesh.msh" with a rectangular shape
$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
1
2 5 "My surface"
$EndPhysicalNames
$Entities
4 4 1 0
1 0 0 0 0
2 0.6 0 0 0
3 0.6 0.04 0 0
4 0 0.04 0 0
1 0 0 0 0.6 0 0 1 1 2 1 -2
2 0.6 0 0 0.6 0.04 0 1 2 2 2 -3
3 0 0.04 0 0.6 0.04 0 1 3 2 3 -4
4 0 0 0 0 0.04 0 1 4 2 4 -1
1 0 0 0 0.6 0.04 0 1 5 4 1 2 3 4
$EndEntities
$Nodes
9 96 1 96
0 1 0 1
1
0 0 0
0 2 0 1
2
0.6 0 0
0 3 0 1
3
0.6 0.04 0
0 4 0 1
4
0 0.04 0
1 1 0 29
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
0.01999999999996189 0 0
0.03999999999992054 0 0
0.05999999999988979 0 0
0.07999999999985184 0 0
0.09999999999979356 0 0
0.1199999999997407 0 0
0.1399999999996824 0 0
0.1599999999996241 0 0
0.1799999999995658 0 0
0.1999999999995075 0 0
0.2199999999994492 0 0
0.2399999999994017 0 0
0.2599999999993434 0 0
0.2799999999992852 0 0
0.2999999999992377 0 0
0.3199999999992748 0 0
0.3399999999993278 0 0
0.3599999999993836 0 0
0.3799999999994396 0 0
0.3999999999994955 0 0
0.4199999999995471 0 0
0.4399999999995949 0 0
0.4599999999996537 0 0
0.4799999999997017 0 0
0.499999999999746 0 0
0.5199999999997968 0 0
0.5399999999998368 0 0
0.5599999999998884 0 0
0.5799999999999441 0 0
1 2 0 1
34
0.6 0.01999999999994783 0
1 3 0 29
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
0.5800000000000303 0.04 0
0.5600000000000506 0.04 0
0.540000000000065 0.04 0
0.5200000000001228 0.04 0
0.5000000000002239 0.04 0
0.4800000000002818 0.04 0
0.4600000000003396 0.04 0
0.4400000000003974 0.04 0
0.4200000000004552 0.04 0
0.4000000000005131 0.04 0
0.3800000000005709 0.04 0
0.360000000000607 0.04 0
0.3400000000006648 0.04 0
0.3200000000007226 0.04 0
0.3000000000007696 0.04 0
0.280000000000732 0.04 0
0.2600000000006786 0.04 0
0.2400000000006222 0.04 0
0.2200000000005658 0.04 0
0.2000000000005094 0.04 0
0.1800000000004574 0.04 0
0.1600000000004089 0.04 0
0.1400000000003496 0.04 0
0.1200000000003012 0.04 0
0.1000000000002564 0.04 0
0.08000000000020513 0.04 0
0.0600000000001647 0.04 0
0.04000000000011261 0.04 0
0.02000000000005631 0.04 0
1 4 0 1
64
0 0.02000000000005288 0
2 1 0 32
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
0.5379917528533162 0.01998513415406266 0
0.03330671636172638 0.01968565975286473 0
0.573941044755064 0.02015520091789318 0
0.07895351594187906 0.0199749808748457 0
0.499072067356839 0.0200121627013391 0
0.119080810547096 0.01998734684171777 0
0.1590820312500057 0.01998625355836172 0
0.4590820312500046 0.02001373008672747 0
0.1990820312499973 0.01998622712727763 0
0.4190820312500133 0.02001266683056389 0
0.2390821838378858 0.01998617861867655 0
0.3790808105468954 0.01999881287096215 0
0.2790927124023353 0.0199855743366759 0
0.3390173339843798 0.01999466841497333 0
0.3192407226562431 0.01998867619492974 0
0.5511484879897564 0.01989454532022246 0
0.0584459085708992 0.01992084936176544 0
0.09938075565530209 0.01999157978131541 0
0.5192419559899495 0.02000239237629373 0
0.1393992614746478 0.01999167757381342 0
0.4793981685758517 0.02000825907470487 0
0.1793994140625054 0.01999148612281571 0
0.4393994140625044 0.02000834086772885 0
0.3993992614746216 0.02000485428370751 0
0.219399433135987 0.01999147537988102 0
0.2994367599487262 0.0199913014316351 0
0.3593887329101658 0.01999357233589276 0
0.2594010734558107 0.01999132014690709 0
0.5836269411603322 0.01211951515143555 0
0.01800898215186078 0.02779700918260565 0
0.01808332327820281 0.01184373289902159 0
0.5835935479344708 0.02808100248725591 0
$EndNodes
$Elements
5 127 1 127
1 1 1 30
1 1 5
2 5 6
3 6 7
4 7 8
5 8 9
6 9 10
7 10 11
8 11 12
9 12 13
10 13 14
11 14 15
12 15 16
13 16 17
14 17 18
15 18 19
16 19 20
17 20 21
18 21 22
19 22 23
20 23 24
21 24 25
22 25 26
23 26 27
24 27 28
25 28 29
26 29 30
27 30 31
28 31 32
29 32 33
30 33 2
1 2 1 2
31 2 34
32 34 3
1 3 1 30
33 3 35
34 35 36
35 36 37
36 37 38
37 38 39
38 39 40
39 40 41
40 41 42
41 42 43
42 43 44
43 44 45
44 45 46
45 46 47
46 47 48
47 48 49
48 49 50
49 50 51
50 51 52
51 52 53
52 53 54
53 54 55
54 55 56
55 56 57
56 57 58
57 58 59
58 59 60
59 60 61
60 61 62
61 62 63
62 63 4
1 4 1 2
63 4 64
64 64 1
2 1 3 63
65 23 76 91 22
66 25 74 88 24
67 27 72 87 26
68 29 69 85 28
69 31 65 83 30
70 60 61 81 68
71 48 49 90 79
72 58 59 82 70
73 56 57 84 71
74 54 55 86 73
75 52 53 89 75
76 50 51 92 77
77 21 78 79 20
78 36 37 65 80
79 32 67 36 80
80 7 81 66 6
81 38 39 69 83
82 9 82 68 8
83 11 84 70 10
84 40 41 72 85
85 13 86 71 12
86 42 43 74 87
87 15 89 73 14
88 44 45 76 88
89 17 92 75 16
90 46 47 78 91
91 19 90 77 18
92 62 63 94 66
93 32 33 93 67
94 77 90 49 50
95 78 21 22 91
96 75 92 51 52
97 76 23 24 88
98 73 89 53 54
99 74 25 26 87
100 71 86 55 56
101 72 27 28 85
102 70 84 57 58
103 68 82 59 60
104 69 29 30 83
105 66 81 61 62
106 65 31 32 80
107 90 19 20 79
108 83 65 37 38
109 91 76 45 46
110 81 7 8 68
111 85 69 39 40
112 89 15 16 75
113 82 9 10 70
114 88 74 43 44
115 87 72 41 42
116 84 11 12 71
117 86 13 14 73
118 79 78 47 48
119 92 17 18 77
120 34 93 33 2
121 64 94 63 4
122 36 67 96 35
123 6 66 95 5
124 64 1 5 95
125 34 3 35 96
126 67 93 34 96
127 66 94 64 95
$EndElements
# --------------------------------------------------------------------------------------------------------------------