A problem of numerical anomaly during calculate

Ian · February 20, 2025, 3:08pm

Hello! I am doing a simulation of a physic model through fenicsx. The model is about the electromagnetic field. I haven’t finished PML outside the domain I set, and I tried to run my semi-finished project. Some of PDEs are time dependent and it’s a mixed finite element problems. I calculate a integration of a variable I need to get on the domain. At first, the integration didn’t change for a few steps, after few steps of calculation, the value of a variable seem to tend to infinity(nan). I guess this is due to the insufficient boundary conditions I set and the I don’t finish PML yet. I have some problems from this:
1.Does the insufficient boundary conditions could cause a numerical anomaly when solve a matrix? How to judge the number of boundary conditions in a system?
2.Implicit iteration is better than explicit iteration, right?
So I want set

F_continuity = (
    ((n - n_n)/dt * v_n
    + ufl.div(n * velocity) * v_n) * ufl.dx
)

rather than

F_continuity = (
    ((n - n_n)/dt * v_n
    + ufl.div(n_n * velocity_n) * v_n) * ufl.dx
)

But

a_compiled = dolfinx.fem.form(a)

doesn’t accept this. Are there some solutions?

I tried to simplify my code, it doesn’t have a nan problem ,but still has the integration didn’t change.

import numpy as np
from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import basix.ufl
import ufl
from dolfinx import mesh, fem
import matplotlib.pyplot as plt
numberofmesh = 2
domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, numberofmesh, numberofmesh, numberofmesh)
T_element = basix.ufl.element("Lagrange", "tetrahedron", 1)
E_element = basix.ufl.element("N1curl", "tetrahedron", 3)
Mix = basix.ufl.mixed_element([T_element, E_element])
W = dolfinx.fem.functionspace(domain, Mix)
(T, E) = ufl.TrialFunctions(W)
(v_T,  v_E) = ufl.TestFunctions(W)
t = dolfinx.fem.Constant(domain, 0.0)
dt = dolfinx.fem.Constant(domain, 1e-15)
w = dolfinx.fem.Function(W)
w_n = dolfinx.fem.Function(W)
T_n, E_n = w_n.split()
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
boundary_dofs_T = dolfinx.fem.locate_dofs_topological(W.sub(0), domain.topology.dim - 1, boundary_facets)
bc_T = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, 300.0), boundary_dofs_T, W.sub(0))
bcs = [bc_T]
w_n.sub(0).interpolate(lambda x: np.full(x.shape[1], 300.0))  # 温度场
w_n.sub(1).interpolate(lambda x: np.zeros((3, x.shape[1])))  # 电场
def top_surface(x):
    return np.isclose(x[2], 1, atol=1)
marker_value=1
top_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, top_surface)
facet_tags = mesh.meshtags(domain, 2, top_facets, np.full(len(top_facets), marker_value, dtype=np.int32))
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
F_heat =  (
        ((T - T_n) / dt * v_T
         + ufl.div(T_n * E_n) * v_T) * ufl.dx
)
F_ampere =  (
        ufl.inner((E - E_n) / dt, v_E) * ufl.dx
        +  ufl.inner(ufl.grad(2 / 3 * T_n), v_E) * ufl.dx
        +  ufl.inner(E_n, v_E) * ufl.dx
        +  ufl.inner(E_n, v_E) * ufl.dx
)
F = F_heat +  F_ampere
a, L = ufl.system(F)
a_compiled = dolfinx.fem.form(a)
L_compiled = dolfinx.fem.form(L)
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps"
}
problem = dolfinx.fem.petsc.LinearProblem(a_compiled, L_compiled, u=w, bcs=bcs, petsc_options=petsc_options)
T_final = 1e-12
num = 0
while t.value < T_final:
    t.value += dt.value
    problem.solve()
    integrate = fem.assemble_scalar(fem.form(T * ufl.dx))
    print(integrate)
    w_n.x.array[:] = w.x.array
    num += 1

The whole code is this, it has a nan problem:

import numpy as np
from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import basix.ufl
import ufl
from dolfinx import mesh, fem
import matplotlib.pyplot as plt
import scipy.sparse
import scipy.sparse.linalg
from dolfinx.io import VTXWriter

numberofmesh = 2
domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, numberofmesh, numberofmesh, numberofmesh)

T_element = basix.ufl.element("Lagrange", "tetrahedron", 1)
n_element = basix.ufl.element("Lagrange", "tetrahedron", 1)
v_element = basix.ufl.element("Lagrange", "tetrahedron", 3, shape=(3,))
E_element = basix.ufl.element("N1curl", "tetrahedron", 3)
H_element = basix.ufl.element("N1curl", "tetrahedron", 3)
Mix = basix.ufl.mixed_element([T_element, n_element, v_element, E_element, H_element])
W = dolfinx.fem.functionspace(domain, Mix)

(T, n, velocity, E, H) = ufl.TrialFunctions(W)
(v_T, v_n, v_v, v_E, v_H) = ufl.TestFunctions(W)

t = dolfinx.fem.Constant(domain, 0.0)

w = dolfinx.fem.Function(W)
w_n = dolfinx.fem.Function(W)
wh = dolfinx.fem.Function(W)
T_n, n_n, velocity_n, E_n, H_n = w_n.split()

D = 1.1e-4
nu = 1e15
epsilon_F = 5.5
m_e = 9.1e-31
e = 1.6e-19
omega_p = 1.4e16
epsilon_0 =  8.854e12
mu_0 = 4 * np.pi * 1e-7
c0 = 299792458.0

Q = dolfinx.fem.functionspace(domain, ("Lagrange", 1))

E_energy = dolfinx.fem.Function(Q, name="Electric_Energy")
H_energy = dolfinx.fem.Function(Q, name="Magnetic_Energy")

total_energy = dolfinx.fem.Function(Q, name="Total_Energy")

bc_value = np.array([0.0, 0.0, 0.0],dtype=float)

domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
boundary_dofs_T = dolfinx.fem.locate_dofs_topological(W.sub(0), domain.topology.dim - 1, boundary_facets)
boundary_dofs_v = dolfinx.fem.locate_dofs_topological(W.sub(2).sub(0), domain.topology.dim - 1, boundary_facets)

bc_T = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, 300.0), boundary_dofs_T, W.sub(0))
bc_v = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, 0.0), boundary_dofs_v, W.sub(2).sub(0))

bcs = [bc_T, bc_v]

w_n.sub(0).interpolate(lambda x: np.full(x.shape[1], 300.0))
w_n.sub(1).interpolate(lambda x: np.full(x.shape[1], 8.8e22))
w_n.sub(2).interpolate(lambda x: np.zeros((3, x.shape[1])))
w_n.sub(3).interpolate(lambda x: np.zeros((3, x.shape[1])))
w_n.sub(4).interpolate(lambda x: np.zeros((3, x.shape[1])))

def laser_function(t):
    tau_p = 100e-15
    return np.exp(-(t - 2*tau_p)**2 / tau_p**2)
laser_source = dolfinx.fem.Constant(domain, 0.0)
dt = dolfinx.fem.Constant(domain, 1e-15)
t.value = 0

def top_surface(x):
    return np.isclose(x[2], 1, atol=1e-8)
marker_value=1
top_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, top_surface)
facet_tags = mesh.meshtags(domain, 2, top_facets, np.full(len(top_facets), marker_value, dtype=np.int32))
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)

F_heat = (
    (T - T_n)/dt * v_T * ufl.dx
    + D * ufl.inner(ufl.grad(T), ufl.grad(v_T)) * ufl.dx
    - (nu * e**2 * 100 / (2 * m_e * omega_p**2)) * laser_source * v_T * ds(1)
)

F_momentum = (
    ufl.inner((velocity - velocity_n)/dt, v_v) * ufl.dx
    + (1/(m_e * n_n)) * ufl.inner(ufl.grad(2/3 * n_n * T_n), v_v) * ufl.dx
    + (e/m_e) * ufl.inner(E_n, v_v) * ufl.dx
    + nu * ufl.inner(velocity_n, v_v) * ufl.dx
)

F_continuity = (
    ((n - n_n)/dt * v_n
    + ufl.div(n_n * velocity_n) * v_n) * ufl.dx
)

F_ampere = (
    ufl.inner(ufl.curl(H_n), v_E) * ufl.dx
    - epsilon_0 * ufl.inner((E - E_n)/dt, v_E) * ufl.dx
    + (4 * np.pi * e/c0) * n_n * ufl.inner(velocity_n, v_E) * ufl.dx
)

F_faraday = (
    ufl.inner(ufl.curl(E_n), v_H) * ufl.dx
    + mu_0 * ufl.inner((H - H_n)/dt, v_H) * ufl.dx
)

F = F_heat + F_momentum + F_continuity + F_ampere + F_faraday

T_final = 1e-12  # 1 ps
num = 0
while t.value < T_final:
    t.value += dt.value
    laser_source.value = laser_function(t.value)
    a, L = ufl.system(F)
    a_compiled = dolfinx.fem.form(a)
    L_compiled = dolfinx.fem.form(L)
    A = dolfinx.fem.create_matrix(a_compiled)
    b = dolfinx.fem.create_vector(L_compiled)
    A_scipy = A.to_scipy()
    dolfinx.fem.assemble_matrix(A, a_compiled, bcs=bcs)
    dolfinx.fem.assemble_vector(b.array, L_compiled)
    dolfinx.fem.apply_lifting(b.array, [a_compiled], [bcs])
    b.scatter_reverse(dolfinx.la.InsertMode.add)
    [bc.set(b.array) for bc in bcs]
    A_inv = scipy.sparse.linalg.splu(A_scipy)
    wh.x.array[:] = A_inv.solve(b.array)

    total_em_energy = fem.assemble_scalar(fem.form(T_n * ufl.dx))
    w_n.x.array[:] = wh.x.array
    num += 1
    print(f"Step {num} | Time: {t.value:.2e} s | Energy: {total_em_energy:.4e} J")

dokken · February 20, 2025, 6:37pm

First of all, this should be replaced with:

w = dolfinx.fem.Function(W)
w_n = dolfinx.fem.Function(W)
T_n, E_n = ufl.split(w_n)

and

Ian:

while t.value < T_final:
    t.value += dt.value
    problem.solve()
    integrate = fem.assemble_scalar(fem.form(T * ufl.dx))
    print(integrate)
    w_n.x.array[:] = w.x.array
    num += 1

with

while t.value < T_final:
    t.value += dt.value
    problem.solve()
    integrate = fem.assemble_scalar(fem.form(w.sub(0) * ufl.dx))
    print(integrate)
    w_n.x.array[:] = w.x.array
    num += 1

For your second issue, please simplify your example and adjust it as done with the first one.

Ian · February 21, 2025, 5:16am

Thank you Prof. dokken.
I simplify the code behind.
It’s strange that integration doesn’t change.

import numpy as np
from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import basix.ufl
import ufl
from dolfinx import mesh, fem
import matplotlib.pyplot as plt
import scipy.sparse
import scipy.sparse.linalg
from dolfinx.io import VTXWriter
numberofmesh = 2
domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, numberofmesh, numberofmesh, numberofmesh)

T_element = basix.ufl.element("Lagrange", "tetrahedron", 1)
E_element = basix.ufl.element("N1curl", "tetrahedron", 3)
Mix = basix.ufl.mixed_element([T_element, E_element])
W = dolfinx.fem.functionspace(domain, Mix)

(T, E) = ufl.TrialFunctions(W)
(v_T,v_E) = ufl.TestFunctions(W)

t = dolfinx.fem.Constant(domain, 0.0)
dt = dolfinx.fem.Constant(domain, 1e-15)
w = dolfinx.fem.Function(W)
w_n = dolfinx.fem.Function(W)
wh = dolfinx.fem.Function(W)
T_n , E_n = ufl.split(w_n)

domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
boundary_dofs_T = dolfinx.fem.locate_dofs_topological(W.sub(0), domain.topology.dim - 1, boundary_facets)
bc_T = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, 300.0), boundary_dofs_T, W.sub(0))
bcs = [bc_T]
w_n.sub(0).interpolate(lambda x: np.full(x.shape[1], 300.0))
w_n.sub(1).interpolate(lambda x: np.zeros((3, x.shape[1])))
def top_surface(x):
    return np.isclose(x[2], 1, atol=1e-8)
marker_value=1
top_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, top_surface)
facet_tags = mesh.meshtags(domain, 2, top_facets, np.full(len(top_facets), marker_value, dtype=np.int32))
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
F_heat = (
    (T - T_n)/dt * v_T * ufl.dx + 10*T*v_T*ufl.dx
)
F_ampere = (
    ufl.inner(ufl.grad(T), v_E) * ufl.dx - ufl.inner((E - E_n)/dt, v_E) * ufl.dx
)
F = F_heat + F_ampere
T_final = 1e-12  # 1 ps
num = 0
a, L = ufl.system(F)
a_compiled = dolfinx.fem.form(a)
L_compiled = dolfinx.fem.form(L)
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps"
}
problem = dolfinx.fem.petsc.LinearProblem(a_compiled, L_compiled, u=wh, bcs=bcs, petsc_options=petsc_options)
while t.value < T_final:
    t.value += dt.value
    problem.solve()
    total_em_energy = fem.assemble_scalar(fem.form(w_n.sub(0) * ufl.dx))
    w_n.x.array[:] = wh.x.array
    num += 1
    print(f"Step {num} | Time: {t.value:.2e} s | integration of T in domain: {total_em_energy:.4e} ")

dokken · February 21, 2025, 5:52am

I guess this should be wh.sub(0).

You are solving for an incredibly small time interval. Imagine multiplying your PDE by st. Then you have
T=T_n + dt*f(T,T_n).

If f(T,T_n) is massive you will still not get a large average chance when multiplied by 1e-15, which is close to machine precision.

Ian · February 21, 2025, 12:49pm

Absolutely a precise insight. I forget I remove the small parameters to simply the code. Then I will check my code again to see what causes integration resulting “nan”.

Topic		Replies	Views
A problem of numerical anomaly（infinity） General	0	6	February 25, 2025
Integral Constraints on PDEs with FEniCSx - Real Finite Elements General	0	296	May 29, 2023
Solution does not meet boundary conditions General	5	380	October 14, 2023
Imposition boundary condition problem General	11	633	December 19, 2022
Error occuring in fem.form dolfinx	9	993	October 2, 2023

A problem of numerical anomaly during calculate

Related topics