Dear FEniCSx Community,
I’m solving the heat equation using FEniCSx. While it works fine on small meshes, I encounter memory issues with larger meshes (e.g., 12+ million elements). Interestingly, solving the potential equation on the same large mesh works without issues.
Is it feasible to run time-dependent simulations on such large meshes in FEniCSx? Are there any recommended strategies to optimize memory usage?
Thank you for your help!
Best regards,
Min
It is feasible to run time dependent simulations.
The memory issue might be down to how you are using FEniCS, and thus to give you any further guidance you would have to provide a reproducible code.
Thank you for your response. I understand the need for a reproducible example. Unfortunately, I cannot share my mesh due to certain constraints, but I can share my code if that helps. Would it still be possible for you to provide guidance based on that? Let me know if I need to include any specific details to make the issue clearer.
Then I would suggest adapting the code to a built in mesh (which you can adjust to appropriate dimensions, ie a unit cube with NxMxO elements in each direction.
Thank you for the suggestion. The code does work with a simple built-in mesh, such as a unit cube, but the issue arises when I run it on the actual larger mesh I am working with. The problem seems to be specific to the size or complexity of the mesh. I can share the code if you’d like to review it for potential optimizations or issues.
Please share the code (and check if having a large enough number of elements in the cube can reproduce your issue).
Certain problems are harder to solve on complex meshes, where iterative solvers might struggle.
But again, a code is needed for any guidance.
i upload my code and mesh there please find it there
many thanks
What version of DOLFINx are you using?
I do not see any memory creep with the mesh that you have supplied (I do note that it does not have 12 million elements).
I’ve slightly modified the code to run with the latest version of DOLFINx:
import dolfinx
import numpy as np
# from dolfinx.fem import assemble
import ufl
from dolfinx.fem import ( functionspace,
assemble_scalar, dirichletbc, form,Constant, locate_dofs_geometrical)
from petsc4py import PETSc
from mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, dx, ds, grad, Measure
from petsc4py.PETSc import ScalarType
from dolfinx import mesh, fem, io, log
from dolfinx.io import XDMFFile
comm = MPI.COMM_WORLD
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
log.set_log_level(log.LogLevel.INFO)
if MPI.COMM_WORLD.rank == 0:
log.log(log.LogLevel.INFO,'Hello')
# Read mesh from file
infile = io.XDMFFile(MPI.COMM_WORLD, "Domain.xdmf", "r")# Read mesh from file
msh = infile.read_mesh(name="Grid")
# Read mesh tags associated with the mesh
meshTags =infile.read_meshtags(msh,name="Grid")
infile.close()
newValues = [-i for i in meshTags.values]
meshTags = mesh.meshtags(msh, meshTags.dim, meshTags.indices, newValues)
msh.topology.create_connectivity(3, 2)
boundaryInfile = io.XDMFFile(MPI.COMM_WORLD, "Domain_boundary.xdmf", "r")
# Read boundary mesh tags associated with the boundary mesh
boundaryMeshTags = boundaryInfile.read_meshtags(msh, name="Grid")
boundaryInfile.close()
# Invert the values of boundary mesh tags in positive sign
newValues = [-i for i in boundaryMeshTags.values]
# Apply the inverted values to boundary mesh tags
boundaryMeshTags = mesh.meshtags(msh, boundaryMeshTags.dim, boundaryMeshTags.indices, newValues)
# Set polynomial degree for the function space
p = 1
V = functionspace(msh, ("Lagrange", p))
DomainTags = {"blood": [6,7.03E-1,1050,3071,0.52],
"tissue": [7, 2.15E-1,1076,3017,0.56],
"board": [8, 0.1,1076,3017,0.518],
"electrode1": [9, 4.6e6,215000,132,71],
"thermorstate": [10,0.2,32,835,0.038] #now i am puttingbne cortical#pervious8.67E-2
}
k_0 = 0.518
c_0 = 3017
sigma_t0 = 0.54
sigma_b0 = 0.1
sigma_ele_0 = 4.6e6
rho = 1040
def C(T):
return c_0 * (1 - 0.0042 * (T - 37))
def k(T):
return k_0 * (1 - 0.0005 * (T - 37))
def sigma_non(T):
return sigma_t0 * (1 + 0.015 * ( T - 37))
def b(T):
return 10e8 * (T - 37 )
def elec(T):
return 10e8 * ( T - 33 )
# solving poteinal equation
Tissue_tags = {"tissue"}
blood_tag = {"blood"}
electrode = {"electrode1"}
boundaryTags = {"face": 20, "bootom": 14 , "left" : 11, "right" : 12, "mid":15, "3rd":16, "up": 13}
patchFacets = boundaryMeshTags.find(boundaryTags["bootom"])
msh.topology.create_connectivity(2, 3)
patchDofs = fem.locate_dofs_topological(V = V, entity_dim = 2, entities = patchFacets)
bcs_1 = fem.dirichletbc(value=ScalarType(0), dofs=patchDofs, V=V)
faceFacets = boundaryMeshTags.find(boundaryTags["face"])
faceDofs = fem.locate_dofs_topological(V = V, entity_dim = 2, entities = faceFacets)
voltage = Constant(msh, ScalarType(2)) # in this way i dont to define boundary condition agsin and again it will work
bcs_2 = fem.dirichletbc(value = voltage, dofs = faceDofs, V = V)
# Define Dirichlet boundary conditions for the patch with zero voltage
bcs =[bcs_1,bcs_2]
f = Constant(msh, 0.0)
dx = ufl.Measure("dx", msh, subdomain_data=meshTags)
ds = ufl.Measure("ds", msh, subdomain_data=boundaryMeshTags)
# Define the trial and test functions
T_initial = 37
phi_n = fem.Function(V, name = "Phi")
# phi_h = fem.Function(V, name="phi_H")
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
T_n = fem.Function(V, name = "T2" )
T_n.x.array[:] = T_initial
# Initialize lists for linear and nonlinear forms
tissue_linear_forms = []
linear_forms = []
# Loop through DomainTags to process each domain
for domain_name, properties in DomainTags.items():
tag = properties[0] # Mesh tag
sigma = properties[1] # Sigma value
if domain_name in Tissue_tags:
# Nonlinear form for tissue
tissue_form = sigma_non(T_n) * dot(grad(u), grad(v)) * dx(tag)
tissue_linear_forms.append(tissue_form)
else:
# Linear form for non-tissue regions
a_phi = sigma * dot(grad(u), grad(v)) * dx(tag)
linear_forms.append(a_phi)
# Combine all forms into a single bilinear form
# here we solve problem 1
a = sum(tissue_linear_forms) + sum(linear_forms)
L = dot(f, v) * dx
problem = fem.petsc.LinearProblem(a, L, bcs = bcs, petsc_options={"ksp_type": "cg",
"pc_type": "hypre",
"ksp_rtol" : 1e-12,
"ksp_atol" : 1e-10,
"ksp_monitor" : None,
# "ksp_converged_reason" : ""
})
if MPI.COMM_WORLD.rank == 0:
print("Starting to solve linear system of size {}".format(problem.b.size))
uh = problem.solve()
uh.name = "u"
uh.x.scatter_forward()
# start doing problem to find total power then we can get desired voltage
# !to get the desired power
voltage_square = voltage.value * voltage.value
remaing = 0
tissue_power = 0
# Loop through DomainTags to process each domain
for domain_name, properties in DomainTags.items():
tag = properties[0] # Mesh tag
sigma = properties[1] # Sigma value
if domain_name in Tissue_tags:
tissue_power += MPI.COMM_WORLD.allreduce(
dolfinx.fem.assemble_scalar(
fem.form(sigma_non(T_n) * dot(grad(uh), grad(uh)) * dx(tag))
),
MPI.SUM
)
else:
remaing += MPI.COMM_WORLD.allreduce(
dolfinx.fem.assemble_scalar(
fem.form(sigma * dot(grad(uh), grad(uh)) * dx(tag))
),
MPI.SUM
)
# Calculate total power and resistance
total_power = remaing + tissue_power
R_ALL = voltage_square / total_power
if MPI.COMM_WORLD.rank == 0:
log.log(log.LogLevel.INFO, f"power_allsubdomain: {total_power}")
log.log(log.LogLevel.INFO, f"resistances of all domain: {R_ALL}")
# we find the total power to get the desired voltage here we wants to do high power shourt duration so total power should be 90 W
p_known = 5
new_voltage = np.sqrt (p_known/ total_power) * voltage.value
voltage.value = new_voltage
log.log(log.LogLevel.INFO, f"new_voltage1: {new_voltage}")
uh = problem.solve()
uh.name = "u"
uh.x.scatter_forward()
phi_n.x.array[:] = uh.x.array[:]
# this part just for checking is it that we want to get the powe did we achived get or not
remaing = 0
tissue_power = 0
# voltage_square = voltage.value * voltage.value
for domain_name, properties in DomainTags.items():
tag = properties[0] #mesh tag
sigma = properties[1]
if domain_name in Tissue_tags:
tissue_power+= MPI.COMM_WORLD.allreduce(dolfinx.fem.assemble_scalar(fem.form( sigma_non(T_n)* dot(grad(uh),grad(uh)) * dx(tag))),MPI.SUM)
else:
remaing+= MPI.COMM_WORLD.allreduce(dolfinx.fem.assemble_scalar(fem.form(sigma * dot(grad(uh),grad(uh)) * dx(tag))),MPI.SUM)
totalpower1 = remaing + tissue_power
R_ALL1 = voltage_square / totalpower1
if MPI.COMM_WORLD.rank == 0:
log.log(log.LogLevel.INFO, f"power_allsubdomain: {totalpower1}")
# start bioheat equation
dt = 0.1
E = grad(phi_n)
v_T = ufl.TestFunction(V)
T_h = fem.Function(V, name="T")
# boundary condition for intaila boundary condition
# this why apply some many boundary 37 becuasue every boundary should have intial temp 37
bio_1Facets = boundaryMeshTags.find(boundaryTags["bootom"])
bio_1Dofs = fem.locate_dofs_topological(V=V, entity_dim=2, entities=bio_1Facets)
bio_bc1 = fem.dirichletbc(value=ScalarType(T_initial), dofs=bio_1Dofs, V=V)
bio_2Facets = boundaryMeshTags.find(boundaryTags["left"])
bio_2Dofs = fem.locate_dofs_topological(V=V, entity_dim=2, entities=bio_2Facets)
bio_bc2 = fem.dirichletbc(value=ScalarType(T_initial), dofs=bio_2Dofs, V=V)
bio_3Facets = boundaryMeshTags.find(boundaryTags["right"])
bio_3Dofs = fem.locate_dofs_topological(V=V, entity_dim=2, entities=bio_3Facets)
bio_bc3 = fem.dirichletbc(value=ScalarType(T_initial), dofs=bio_3Dofs, V=V)
bio_4Facets = boundaryMeshTags.find(boundaryTags["mid"])
bio_4Dofs = fem.locate_dofs_topological(V=V, entity_dim=2, entities=bio_4Facets)
bio_bc4 = fem.dirichletbc(value=ScalarType(T_initial), dofs=bio_4Dofs, V=V)
bio_5Facets = boundaryMeshTags.find(boundaryTags["3rd"])
bio_5Dofs = fem.locate_dofs_topological(V=V, entity_dim=2, entities=bio_5Facets)
bio_bc5 = fem.dirichletbc(value=ScalarType(T_initial), dofs=bio_5Dofs, V=V)
bio_6Facets = boundaryMeshTags.find(boundaryTags["up"])
bio_6Dofs = fem.locate_dofs_topological(V=V, entity_dim=2, entities=bio_6Facets)
bio_bc6 = fem.dirichletbc(value=ScalarType(T_initial), dofs=bio_6Dofs, V=V)
bcs_b = [bio_bc1,bio_bc2,bio_bc3,bio_bc4,bio_bc5,bio_bc6]
# weak formulation for bioheat equation
linear_F = 0
non_linear_F = 0
blood_F = 0
electrode_F = 0
# Loop through domain tags to process each region
for domain_name, properties in DomainTags.items():
tag = properties[0] # Mesh tag
sigma = properties[1] # Conductivity
rho = properties[2] # Density
c = properties[3] # Specific heat capacity
k_t = properties[4] # Thermal conductivity
if domain_name in Tissue_tags:
# Form for tissue tags (dependent on temperature)
non_linear_F += (
C(T_h) * rho * ((T_h - T_n) / dt) * v_T * dx(tag)
+ k(T_h) * dot(grad(T_h), grad(v_T)) * dx(tag)
- sigma_non(T_n) * dot(E, E) * v_T * dx(tag)
)
elif domain_name in blood_tag:
# Form for blood tag (may include blood-specific terms)
blood_F += (
rho * c * ((T_h - T_n) / dt) * v_T * dx(tag)
+ k_t * dot(grad(T_h), grad(v_T)) * dx(tag) + b(T_h) * v_T * dx(tag)
- sigma * dot(E, E) * v_T * dx(tag)
)
elif domain_name in electrode:
# Form for electrode tag (specific to electrode)
electrode_F += (
rho * c * ((T_h - T_n) / dt) * v_T * dx(tag)
+ k_t * dot(grad(T_h), grad(v_T)) * dx(tag) + elec(T_h) * v_T * dx(tag)
- sigma * dot(E, E) * v_T * dx(tag)
)
else:
# Linear form for non-temperature-dependent regions
linear_F += (
rho * c * ((T_h - T_n) / dt) * v_T * dx(tag)
+ k_t * dot(grad(T_h), grad(v_T)) * dx(tag)
- sigma * dot(E, E) * v_T * dx(tag)
)
# Combine all contributions
F2 = non_linear_F + linear_F + blood_F + electrode_F
problem2 = NonlinearProblem(F2, T_h, bcs = bcs_b)
solver1 = NewtonSolver(MPI.COMM_WORLD, problem2)
solver1.convergence_criterion = "incremental"
solver1.rtol = 1e-6
solver1.report = True
ksp = solver1.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}ksp_rtol"] = 1.0e-8
opts[f"{option_prefix}pc_type"] = "asm"
# opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
opts[f"{option_prefix}ksp_monitor"] = None
t = 0
with io.XDMFFile(msh.comm, "ioheat1.xdmf", "w") as file:
file.write_mesh(msh)
# file.write_function(uh)
T = 10.2 #high power shourt duration
while t < T:
t += dt
# ist bio heat then
n , converged = solver1.solve(T_h)
assert (converged)
T_h.x.scatter_forward()
#bioheat equation # for updating tempe
T_n.x.array[:] = T_h.x.array[:]
log.log(log.LogLevel.INFO,'istproblem finish')
# poteinal equation for p1
uh = problem.solve()
uh.x.scatter_forward()
uh.name = "u"
total_tissue_power = 0
remaing = 0
# for geting the desired power
remaing = 0
tissue_power = 0
# Loop through DomainTags to process each domain
for domain_name, properties in DomainTags.items():
tag = properties[0] # Mesh tag
sigma = properties[1] # Sigma value
if domain_name in Tissue_tags:
# Compute and accumulate power for tissue-related tags
tissue_power += MPI.COMM_WORLD.allreduce(
dolfinx.fem.assemble_scalar(
fem.form(sigma_non(T_n) * dot(grad(uh), grad(uh)) * dx(tag))
),
MPI.SUM
)
else:
# Compute and accumulate power for non-tissue regions
remaing += MPI.COMM_WORLD.allreduce(
dolfinx.fem.assemble_scalar(
fem.form(sigma * dot(grad(uh), grad(uh)) * dx(tag))
),
MPI.SUM
)
totalpower2 = remaing + total_tissue_power
new_voltage1 = np.sqrt ( p_known / totalpower2 ) * voltage.value
voltage.value = new_voltage1
if MPI.COMM_WORLD.rank == 0:
log.log(log.LogLevel.INFO,'2ndproblem finish')
log.log(log.LogLevel.INFO, f"power_allsubdomain2nd: {totalpower2}")
# log.log(log.LogLevel.INFO, f"resistances of all domain2nd: {R_ALL2}")
uh = problem.solve()
uh.x.scatter_forward()
phi_n.x.array[:] = uh.x.array[:]
if MPI.COMM_WORLD.rank == 0:
log.log(log.LogLevel.INFO, f"Total Power: {totalpower2}, New Voltage: {new_voltage1}")
log.log(log.LogLevel.INFO, f"Time step: {t}, Temperature: {T_h.x.array[:]}")
# print(f"Number of interations: {n:d}")
file = XDMFFile(comm,"ioheat1.xdmf" , "a")
file.write_function(T_h,t)
file.write_function(uh,t)
file.close()
However, please do note that your example is rather complicated and takes a lot of effort to parse for anyone. Have you tried running just a temporal evolving heat equation (time dependent diffusion eq) on the same mesh and seen if you experience the same memory creep.
Thank you for your reply. However, I did not share my actual mesh. The mesh I shared with you was working with the code I sent, but my actual mesh is not working when even I try to run a simple heat equation.
It would then be useful to have the actual mesh with a simple heat-equation, illustrating the memory issue. I cannot help any further without a code that reproduces the issue.
I have added my actual mesh with a very simple code. Could you please help me identify the issue? It works for the potential equation but not for the time-dependent problem.
You have supplied a mesh with 100 million cells and 16 millions degrees of freedom.
I do not have a computing resource that I can use (without paying for it) to run that example.
For me, when trying to run the code on the resources I have available, the first solve uses about 69.5 GB RAM.
I get the output:
[2025-01-31 16:52:23.054] [info] power_allsubdomain: 0.022023461823258966
[2025-01-31 16:52:23.054] [info] resistances of all domain: 181.6244890154191
[2025-01-31 16:52:23.054] [info] new_voltage1: 30.135070019449028
The next solve uses the same amount of memory.
You have yourself not provided any number on the memory creep that you are seeing, and your code involves way more than just a single temporal problem (as it involves a linear problem and a nonlinear problem).
The first solve within the temporal loop consumes 86.4 GB RAM on my system.
Within the loop memory rises to my limit (252 GB Ram).
It is unclear where in the process this happens (likely within the PETSc solver).
For any further advice, you would need to do some proper profiling yourself, and work on isolating the issue you are having.
I am fully aware of the details regarding DoFs and element definitions—they are identical to the first problem, which works perfectly. The issue is specifically with the second problem, where memory usage grows uncontrollably.
If you cannot provide any meaningful help in diagnosing this issue, then please let me know where I should ask instead. Would the PETSc community be the right place, or is there another more relevant forum?
I’d appreciate a clear answer.
Are you sure it is not the solver setting for the second problem that is the issue?
I haven’t used asm, but it seems like it should be configured PCASM — PETSc 3.22.3 documentation
See for instance https://petsc.org/main/src/ksp/ksp/tutorials/ex8.c.html
Frequently Asked Questions (FAQ) — PETSc 3.15.0 documentation
Thanks for your replay.
I don’t think this is a PETSc issue because I changed the preconditioner, but I am still getting the same error. I suspect that FEniCSx is struggling with handling a large mesh or there might be another issue. I’m not sure what the exact problem is. How can I resolve this?
It is quite clear to me that this is a PETSc/solver issue, as the two linear solves works quite well (only uses 80 GB ram).
The huge RAM usage happens when one hits the non-linear solver (as far as I can tell).
The only thing that happens within the non-linear solver, compared to a linear problem, is that one computes the symbolic derivative of the residual to form the jacobian. This computation is independent of the mesh size.
The next thing that happens is the assembly into the Jacobian matrix and residual vector, which uses the exact same routines as a linear-problem (i.e. dolfinx.fem.petsc.assemble_*.
The only place memory usage post this can happen, is the solve at a given Newton step. This is then dependent on the solver, and the sub-preconditioner used for each partition of the multigrid.
In previous posts, you have blamed the time-dependent nature of the problem, but as far as I can tell from trying to run this with the resources I have you don’t get past the first Newton iteration. Have you inspected the output one can get from the DOLFINx Newton solver about the iterations and residuals?