Error and I can't follow up with new fixes

Errors
1

I’m trying to learn python and get used to FenicsX, But I’m encountering a lot of issues, I’m trying to simulate 3d Wave Equation in a Steel Plate

import numpy as np
import dolfinx
from dolfinx import fem, mesh
import ufl
from mpi4py import MPI
import dolfinx.fem.petsc
from petsc4py import PETSc
import basix
import matplotlib.pyplot as plt

# --- 1. Material Properties ---
density = 7850  # kg/m^3
E = 210e9  # Young's modulus in Pascals

# --- 2. Mesh Generation ---
length = 1.0  # Length of the plate in meters
width = 0.5   # Width of the plate in meters
height = 0.01  # Thickness of the plate in meters

# Creating the mesh
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 20, 10, 2, cell_type=dolfinx.mesh.CellType.hexahedron)
mesh.geometry.x[:, 0] *= length  # Scale in x-direction
mesh.geometry.x[:, 1] *= width   # Scale in y-direction
mesh.geometry.x[:, 2] *= height  # Scale in z-direction

# --- 3. Function Space ---
element = basix.ufl.element("Lagrange", basix.CellType.hexahedron, 1, shape=(3,))
V = dolfinx.fem.functionspace(mesh, element)

# --- 4. Boundary Conditions ---
def boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], length)

u_bc = fem.Function(V)
u_bc.x.array[:] = 0.0
bc = fem.dirichletbc(u_bc, fem.locate_dofs_geometrical(V, boundary))

# --- 5. Time Parameters ---
T = 0.5  # Total time in seconds
num_steps = 100  # Number of time steps
dt = T / num_steps  # Time step size

# --- 6. Wave Speed ---
c = np.sqrt(E / density)
c_squared = fem.Constant(mesh, PETSc.ScalarType(c**2).real)
dt_sq = fem.Constant(mesh, PETSc.ScalarType(dt**2).real)

# --- 7. Trial and Test Functions ---
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# --- 8. Initial Conditions ---
u_n = fem.Function(V)
u_n1 = fem.Function(V)

def initial_condition(x):
    values = np.zeros((x.shape[0], 3), dtype=PETSc.ScalarType)
    values[:, 0] = np.sin(np.pi * x[:, 0]) * np.sin(2 * np.pi * x[:, 1]) 
    return values

u_n1.interpolate(initial_condition)

# --- 9. External Force ---
f = fem.Constant(mesh, np.array([0.0, 0.0, 0.0])) 

# --- 10. Variational Formulation ---
dx = ufl.Measure("dx", domain=mesh)

# Define the variational forms
a = ufl.inner(u, v) * dx + dt_sq * c_squared * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
L = dt_sq * ufl.inner(f, v) * dx + ufl.inner(2 * u_n - u_n1, v) * dx


# --- 11. Time-Stepping Loop ---
u = fem.Function(V)
for n in range(num_steps):
    petsc_options = {"ksp_type": "cg", "pc_type": "ilu"} 
    problem = dolfinx.fem.petsc.LinearProblem(a, L, u=u, bcs=[bc], petsc_options=petsc_options)
    problem.solve()

    if n % 10 == 0:
        with dolfinx.io.XDMFFile(MPI.COMM_WORLD, f"displacement_3d_{n}.xdmf", "w") as file:
            file.write_mesh(mesh)
            file.write_function(u, n*dt)

    # Update solution 
    u_n1.x.array[:] = u_n.x.array
    u_n.x.array[:] = u.x.array


# --- 12. Plotting ---
# Get coordinates of the mesh nodes
x = mesh.geometry.x

# Reshape the displacement array (extract the x-component)
u_array = u.x.array.reshape(-1, 3)[:, 0]

# Get the number of nodes in x and y directions
num_nodes_x = 20  
num_nodes_y = 10  
u_array = u_array.reshape(num_nodes_y, num_nodes_x)

num_nodes_x = mesh.topology.index_map(0).size_local
num_nodes_y = mesh.topology.index_map(1).size_local

# Create the plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(np.linspace(0, length, num_nodes_x), np.linspace(0, width, num_nodes_y))
ax.plot_surface(X, Y, u_array, cmap='viridis')
ax.set_xlabel('X coordinate (m)')
ax.set_ylabel('Y coordinate (m)')
ax.set_zlabel('Displacement (m)')
ax.set_title('Displacement of the Steel Plate at Final Time')
plt.show()

And I get this error
TypeError Traceback (most recent call last)
Cell In[5], line 79
77 for n in range(num_steps):
78 petsc_options = {“ksp_type”: “cg”, “pc_type”: “ilu”}
—> 79 problem = dolfinx.fem.petsc.LinearProblem(a, L, u=u, bcs=[bc], petsc_options=petsc_options)
80 problem.solve()
82 if n % 10 == 0:

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/petsc.py:782, in LinearProblem.init(self, a, L, bcs, u, petsc_options, form_compiler_options, jit_options)
744 def init(
745 self,
746 a: ufl.Form,
(…)
752 jit_options: typing.Optional[dict] = None,
753 ):
754 “”“Initialize solver for a linear variational problem.
755
756 Args:
(…)
780 “mumps”})
781 “””
→ 782 self._a = _create_form(
783 a,
784 dtype=PETSc.ScalarType,
785 form_compiler_options=form_compiler_options,
786 jit_options=jit_options,
787 )
788 self._A = create_matrix(self._a)
789 self._L = _create_form(
790 L,
791 dtype=PETSc.ScalarType,
792 form_compiler_options=form_compiler_options,
793 jit_options=jit_options,
794 )

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/forms.py:337, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
334 else:
335 return form
→ 337 return _create_form(form)

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/forms.py:331, in form.._create_form(form)
329 return None
330 else:
→ 331 return _form(form)
332 elif isinstance(form, collections.abc.Iterable):
333 return list(map(lambda sub_form: _create_form(sub_form), form))

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/forms.py:306, in form.._form(form)
303 else:
304 _entity_maps = {msh._cpp_object: emap for (msh, emap) in entity_maps.items()}
→ 306 f = ftype(
307 module.ffi.cast(“uintptr_t”, module.ffi.addressof(ufcx_form)),
308 V,
309 coeffs,
310 constants,
311 subdomains,
312 _entity_maps,
313 mesh,
314 )
315 return Form(f, ufcx_form, code, module)

TypeError: init(): incompatible function arguments. The following argument types are supported:
1. init(self, spaces: collections.abc.Sequence[dolfinx.cpp.fem.FunctionSpace_float64], integrals: collections.abc.Mapping[dolfinx.cpp.fem.IntegralType, collections.abc.Sequence[tuple[int, int, ndarray[dtype=int32, shape=(), order=‘C’, writable=False], ndarray[dtype=int8, shape=(), order=‘C’, writable=False]]]], coefficients: collections.abc.Sequence[dolfinx.cpp.fem.Function_complex128], constants: collections.abc.Sequence[dolfinx.cpp.fem.Constant_complex128], need_permutation_data: bool, entity_maps: collections.abc.Mapping[dolfinx.cpp.mesh.Mesh_float64, ndarray[dtype=int32, shape=(), order=‘C’, writable=False]], mesh: dolfinx.cpp.mesh.Mesh_float64 | None) → None
2. init(self, form: int, spaces: collections.abc.Sequence[dolfinx.cpp.fem.FunctionSpace_float64], coefficients: collections.abc.Sequence[dolfinx.cpp.fem.Function_complex128], constants: collections.abc.Sequence[dolfinx.cpp.fem.Constant_complex128], subdomains: collections.abc.Mapping[dolfinx.cpp.fem.IntegralType, collections.abc.Sequence[tuple[int, ndarray[dtype=int32, shape=(), order=‘C’, writable=False]]]], entity_maps: collections.abc.Mapping[dolfinx.cpp.mesh.Mesh_float64, ndarray[dtype=int32, shape=(*), order=‘C’, writable=False]], mesh: dolfinx.cpp.mesh.Mesh_float64 | None) → None

Invoked with types: dolfinx.cpp.fem.Form_complex128, _CDataBase, list, list, list, dict, dict, dolfinx.cpp.mesh.Mesh_float64

2

I would change all dtype definitions to dolfinx.default_scalar_type, i.e.

c_squared = fem.Constant(mesh, dolfinx.default_scalar_type(c**2))
dt_sq = fem.Constant(mesh, dolfinx.default_scalar_type(dt**2))

# --- 7. Trial and Test Functions ---
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# --- 8. Initial Conditions ---
u_n = fem.Function(V)
u_n1 = fem.Function(V)

def initial_condition(x):
    values = np.zeros((x.shape[0], 3), dtype=dolfinx.default_scalar_type)
    values[:, 0] = np.sin(np.pi * x[:, 0]) * np.sin(2 * np.pi * x[:, 1]) 
    return values

u_n1.interpolate(initial_condition)

# --- 9. External Force ---
f = fem.Constant(mesh, np.array([0.0, 0.0, 0.0], dtype=dolfinx.default_scalar_type))

that will bring you all the way down to postprocessing
where

this doesn’t make sense to me.

3

For the displacement data (u_array), and reorganizes it into a grid (num_nodes_y rows by num_nodes_x columns) that matches the structure of the mesh and especially for the plotting. Please correct me if I’m wrong

6

I have made a few changes but, There’s no displacement values

import numpy as np
import dolfinx
from dolfinx import fem, mesh
import ufl
from mpi4py import MPI
import dolfinx.fem.petsc
from petsc4py import PETSc
import basix
import matplotlib.pyplot as plt

# --- 1. Material Properties ---
density = 7850  # kg/m^3
E = 210e9  # Young's modulus in Pascals

# --- 2. Mesh Generation ---
length = 1.0  # Length of the plate in meters
width = 0.5   # Width of the plate in meters
height = 0.01  # Thickness of the plate in meters

# Creating the mesh
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 20, 10, 2, cell_type=dolfinx.mesh.CellType.tetrahedron)
mesh.geometry.x[:, 0] *= length  # Scale in x-direction
mesh.geometry.x[:, 1] *= width   # Scale in y-direction
mesh.geometry.x[:, 2] *= height  # Scale in z-direction

# Create edge entities (dimension 1)
mesh.topology.create_entities(1)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "my_mesh.xdmf", "w") as file:
    file.write_mesh(mesh)

# --- 3. Function Space ---
element = basix.ufl.element("Lagrange", basix.CellType.tetrahedron, 1, shape=(3,))
V = dolfinx.fem.functionspace(mesh, element)

# --- 4. Boundary Conditions ---
def boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], length)

u_bc = fem.Function(V)
u_bc.x.array[:] = 0.0
bc = fem.dirichletbc(u_bc, fem.locate_dofs_geometrical(V, boundary))

# --- 5. Time Parameters ---
T = 1.0  # Total time in seconds
num_steps = 100  # Number of time steps
dt = T / num_steps  # Time step size

# --- 6. Wave Speed ---
c = np.sqrt(E / density)
c_squared = fem.Constant(mesh, dolfinx.default_scalar_type(c**2))
dt_sq = fem.Constant(mesh, dolfinx.default_scalar_type(dt**2))

# --- 7. Trial and Test Functions ---
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# --- 8. Initial Conditions ---
u_n = fem.Function(V)
u_n1 = fem.Function(V)

def initial_condition(x):
    values = np.zeros((x.shape[0], 3), dtype=dolfinx.default_scalar_type)
    center_x = length / 2
    center_y = width / 2
    sigma = 0.1  
    amplitude = 0.01  
    values[:, 0] = amplitude * np.exp(-((x[:, 0] - center_x)**2 + (x[:, 1] - center_y)**2) / (2 * sigma**2))
    return values

u_n1.interpolate(initial_condition)

# --- 9. External Force ---
f = fem.Constant(mesh, np.array([0.0, 0.0, -9.81 * density], dtype=dolfinx.default_scalar_type)) 

# --- 10. Variational Formulation ---
dx = ufl.Measure("dx", domain=mesh)

a = ufl.inner(u, v) * dx + dt_sq * c_squared * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
L = dt_sq * ufl.inner(f, v) * dx + ufl.inner(2 * u_n - u_n1, v) * dx

# --- 11. Time-Stepping Loop ---
u = fem.Function(V)
for n in range(num_steps):
    petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}

    # Initialize LinearProblem.
    problem = dolfinx.fem.petsc.LinearProblem(a, L, u=u, bcs=[bc], petsc_options=petsc_options)
    
    problem.solve()

    # Check displacement values every few steps to debug.
    if n % 10 == 0:
        print(f"Step {n}, Displacement max: {np.max(u.x.array)}, min: {np.min(u.x.array)}")

        with dolfinx.io.XDMFFile(MPI.COMM_WORLD, f"displacement_3d_{n}.xdmf", "w") as file:
            file.write_mesh(mesh)
            file.write_function(u, n*dt)

    # Update solution 
    u_n1.x.array[:] = u_n.x.array
    u_n.x.array[:] = u.x.array

# --- 12. Plotting ---
x = mesh.geometry.x
u_array = u.x.array.reshape(-1, 3)[:, 0]

num_nodes_x = mesh.topology.index_map(0).size_local  
total_nodes = u_array.size  
num_nodes_y = total_nodes // num_nodes_x  

if total_nodes % num_nodes_x != 0:
    raise ValueError("Total nodes are not divisible by number of nodes in x direction.")

X, Y = np.meshgrid(np.linspace(0, length, num_nodes_x), np.linspace(0, width, num_nodes_y))
Z_array = u_array.reshape(num_nodes_y, num_nodes_x)  

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z_array, cmap='viridis')
ax.set_xlabel('X coordinate (m)')
ax.set_ylabel('Y coordinate (m)')
ax.set_zlabel('Displacement (m)')
ax.set_title('Displacement of the Steel Plate at Final Time')
plt.show()

Step 0, Displacement max: (nan+nanj), min: (nan+nanj)
Step 10, Displacement max: (nan+nanj), min: (nan+nanj)
Step 20, Displacement max: (nan+nanj), min: (nan+nanj)
Step 30, Displacement max: (nan+nanj), min: (nan+nanj)
Step 40, Displacement max: (nan+nanj), min: (nan+nanj)
Step 50, Displacement max: (nan+nanj), min: (nan+nanj)
Step 60, Displacement max: (nan+nanj), min: (nan+nanj)
Step 70, Displacement max: (nan+nanj), min: (nan+nanj)
Step 80, Displacement max: (nan+nanj), min: (nan+nanj)
Step 90, Displacement max: (nan+nanj), min: (nan+nanj)

Any help?

7

Try adding
"pc_factor_mat_solver_type": "mumps"
and
"ksp_error_if_not_converged": True
to petsc_options get some output from the petsc ksp solver (and use a more stable LU solver

8

Tried that

import numpy as np
import dolfinx
from dolfinx import fem, mesh
import ufl
from mpi4py import MPI
import dolfinx.fem.petsc
from petsc4py import PETSc
import basix
import matplotlib.pyplot as plt

# --- 1. Material Properties ---
density = 7850  # kg/m^3
E = 210e9  # Young's modulus in Pascals

# --- 2. Mesh Generation ---
length = 1.0  # Length of the plate in meters
width = 0.5   # Width of the plate in meters
height = 0.01  # Thickness of the plate in meters

# Creating the mesh
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 20, 10, 2, cell_type=dolfinx.mesh.CellType.tetrahedron)
mesh.geometry.x[:, 0] *= length  # Scale in x-direction
mesh.geometry.x[:, 1] *= width   # Scale in y-direction
mesh.geometry.x[:, 2] *= height  # Scale in z-direction

# Create edge entities (dimension 1)
mesh.topology.create_entities(1)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "my_mesh.xdmf", "w") as file:
    file.write_mesh(mesh)

# --- 3. Function Space ---
element = basix.ufl.element("Lagrange", basix.CellType.tetrahedron, 1, shape=(3,))
V = dolfinx.fem.functionspace(mesh, element)

# --- 4. Boundary Conditions ---
def boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], length)

u_bc = fem.Function(V)
u_bc.x.array[:] = 0.0  # Ensure real values for boundary conditions
bc = fem.dirichletbc(u_bc, fem.locate_dofs_geometrical(V, boundary))

print(f"Boundary condition values: {u_bc.x.array}")

# --- 5. Time Parameters ---
T = 1.0  # Total time
num_steps = 700  # Increase the number of steps for smaller time step size
dt = T / num_steps  # Smaller time step

# --- 6. Wave Speed ---
c = np.sqrt(E / density)
c_squared = fem.Constant(mesh, dolfinx.default_scalar_type(c**2))
dt_sq = fem.Constant(mesh, dolfinx.default_scalar_type(dt**2))
c = np.sqrt(E / density)
print(f"Wave speed: {c}")

# --- 7. Trial and Test Functions ---
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# --- 8. Initial Conditions ---
u_n = fem.Function(V)
u_n1 = fem.Function(V)

def initial_condition(x):
    values = np.zeros((x.shape[0], 3), dtype=PETSc.ScalarType)
    center_x = length / 2
    center_y = width / 2
    sigma = 0.1
    amplitude = 0.01
    values[:, 2] = amplitude * np.exp(-((x[:, 0] - center_x)**2 + (x[:, 1] - center_y)**2) / (2 * sigma**2))  # Z-direction only
    return values

u_n.interpolate(initial_condition)
u_n1.interpolate(initial_condition)

print(f"Initial max displacement (u_n): {np.max(u_n.x.array)}, min displacement (u_n): {np.min(u_n.x.array)}")

# --- 9. External Force ---
f = fem.Constant(mesh, np.array([0.0, 0.0, -9.81 * density], dtype=dolfinx.default_scalar_type)) 

# --- 10. Variational Formulation ---
dx = ufl.Measure("dx", domain=mesh)

a = ufl.inner(u, v) * dx + dt_sq * c_squared * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
L = ufl.inner(2 * u_n - u_n1, v) * dx  # Temporarily remove the external force for debugging

print(f"L form: {L}")


# --- 11. Time-Stepping Loop ---
u = fem.Function(V)
for n in range(num_steps):
    # Define PETSc solver options with MUMPS and error reporting
    petsc_options = {
        "ksp_type": "preonly",  # Direct solver
        "pc_type": "lu",  # LU decomposition
        "pc_factor_mat_solver_type": "mumps",  # Use MUMPS solver
        "ksp_error_if_not_converged": True  # Raise error if the solver doesn't converge
    }

    # Initialize the LinearProblem.
    problem = dolfinx.fem.petsc.LinearProblem(a, L, u=u, bcs=[bc], petsc_options=petsc_options)
    
    # Solve the linear problem
    problem.solve()

    # Check displacement values every few steps to debug.
    if n % 10 == 0:
        print(f"Step {n}, Displacement max: {np.max(u.x.array)}, min: {np.min(u.x.array)}")

        # Write displacement field to file
        with dolfinx.io.XDMFFile(MPI.COMM_WORLD, f"displacement_3d_{n}.xdmf", "w") as file:
            file.write_mesh(mesh)
            file.write_function(u, n * dt)

    # Update solution for the next time step
    u_n.x.array[:] = np.real(u_n.x.array)
    u_n1.x.array[:] = np.real(u_n1.x.array)

print(f"Initial max displacement (u_n): {np.max(u_n.x.array)}, min displacement (u_n): {np.min(u_n.x.array)}")
print(f"Initial max displacement (u_n1): {np.max(u_n1.x.array)}, min displacement (u_n1): {np.min(u_n1.x.array)}")

# --- 12. Plotting ---
x = mesh.geometry.x
u_array = u.x.array.reshape(-1, 3)[:, 0]

num_nodes_x = mesh.topology.index_map(0).size_local  
total_nodes = u_array.size  
num_nodes_y = total_nodes // num_nodes_x  

print(f"Step {n}, Max displacement: {np.max(u.x.array)}, Min displacement: {np.min(u.x.array)}")

if total_nodes % num_nodes_x != 0:
    raise ValueError("Total nodes are not divisible by number of nodes in x direction.")

X, Y = np.meshgrid(np.linspace(0, length, num_nodes_x), np.linspace(0, width, num_nodes_y))
Z_array = u_array.reshape(num_nodes_y, num_nodes_x)  

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z_array, cmap='viridis')
ax.set_xlabel('X coordinate (m)')
ax.set_ylabel('Y coordinate (m)')
ax.set_zlabel('Displacement (m)')
ax.set_title('Displacement of the Steel Plate at Final Time')
plt.show()

Gives this and there’s nothing in the plot

Screenshot 2024-10-23 010052
Screenshot 2024-10-23 010052509×379 27.8 KB

Boundary condition values: [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]
Wave speed: 5172.194153034851
Initial max displacement (u_n): (4.093e-320+nanj), min displacement (u_n): (4.093e-320+nanj)
L form: { conj(((v_0) : (({ A | A_{i_{58}} = 2 * f[i_{58}] }) + ({ A | A_{i_{59}} = -1 * f[i_{59}] })))) } * dx(<Mesh #19>[everywhere], {})
Step 0, Displacement max: (nan+nanj), min: (nan+nanj)
Step 10, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 20, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 30, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 40, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 50, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 60, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 70, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 80, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 90, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 100, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 110, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 120, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 130, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 140, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 150, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 160, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 170, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 180, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 190, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 200, Displacement max: (1.3635697704911794e+290+0j), min: 0j