Hi, here is the minimal code to reproduce it
file mesh_test.py
import dolfinx
from dolfinx import mesh, cpp
from mpi4py import MPI
import pyvista as pv
import numpy as np
def create_nested_rectangles():
# Create outer rectangle domain
domain = mesh.create_rectangle(
comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (0.4, 0.4)),
n=(32, 32)
)
def inner_domain(x):
return np.logical_and(
np.logical_and(x[0] >= 0.1, x[0] <= 0.3),
np.logical_and(x[1] >= 0.15, x[1] <= 0.25)
)
dim = domain.topology.dim
domain.topology.create_entities(dim)
inner_cells = mesh.locate_entities(domain, dim, inner_domain)
# Create default values for all cells
num_cells = domain.topology.index_map(dim).size_local
markers_array = np.zeros(num_cells, dtype=np.int32)
# Set values for inner domain cells
markers_array[inner_cells] = 1
# Create MeshTags object
markers = mesh.meshtags(domain, dim, np.arange(num_cells), markers_array)
return domain, markers
def plot_mesh_with_subdomains(domain, markers):
# Convert to VTK format for plotting
topology, cells, geometry = dolfinx.plot.vtk_mesh(domain)
grid = pv.UnstructuredGrid(topology, cells, geometry)
# Add markers as cell data
grid.cell_data["markers"] = markers.values
# Create plotter
plotter = pv.Plotter(shape=(1, 2))
# Plot full mesh
plotter.subplot(0, 0)
plotter.add_mesh(grid, show_edges=True)
plotter.add_title('Complete Mesh')
# Plot mesh with subdomains
plotter.subplot(0, 1)
plotter.add_mesh(grid, show_edges=True, scalars="markers")
plotter.add_title('Mesh with Subdomains')
plotter.show()
# Create and plot the mesh
domain, markers = create_nested_rectangles()
plot_mesh_with_subdomains(domain, markers)
file A0_f_test.py
from dolfinx import fem
import ufl
import numpy as np
from ufl import as_ufl
def create_coefficient_matrices_test(p_f, u_f, T_f, domain):
# Define constants
R = fem.Constant(domain, np.float64(287.05)) # Specific gas constant for air
c_v = fem.Constant(domain, np.float64(718.0)) # Specific heat at constant volume
# Define basic thermodynamic expressions directly
rho_f_expr = p_f / (as_ufl(R) * T_f)
v_expr = as_ufl(R) * T_f / p_f
e_expr = as_ufl(c_v) * T_f
h_expr = e_expr + p_f * v_expr
u_squared = 0.5 * (u_f[0] * u_f[0] + u_f[1] * u_f[1])
e_tot_expr = e_expr + u_squared
# Define thermodynamic coefficients directly
beta_T_expr = 1.0 / p_f
alpha_p_expr = 1.0 / T_f
# Create coefficient matrices using ufl.as_matrix
get_A0 = ufl.as_matrix([
[rho_f_expr * beta_T_expr, fem.Constant(domain, np.float64(0)), fem.Constant(domain, np.float64(0)),
-rho_f_expr * alpha_p_expr],
[rho_f_expr * beta_T_expr * u_f[0], rho_f_expr, fem.Constant(domain, np.float64(0)),
-rho_f_expr * alpha_p_expr * u_f[0]],
[rho_f_expr * beta_T_expr * u_f[1], fem.Constant(domain, np.float64(0)), rho_f_expr,
-rho_f_expr * alpha_p_expr * u_f[1]],
[(e_expr - fem.Constant(domain, np.float64(1))), rho_f_expr * u_f[0], rho_f_expr * u_f[1],
fem.Constant(domain, np.float64(0))]
])
return get_A0
file assemble_test.py
from dolfinx import fem
from dolfinx.fem import *
from ufl import *
from dolfinx.fem.petsc import assemble_matrix
def assemble_weak_form_test(V_fluid, dt, A0_f,markers):
"""
Assemble weak form for FSI problem in DOLFINx
"""
# Test and trial functions
W = TestFunction(V_fluid)
Y=TrialFunction(V_fluid)
# Mesh information
mesh = V_fluid.mesh
# Define measures with markers
dx = Measure("dx", domain=mesh)
dx_s = Measure("dx", domain=mesh, subdomain_data=markers, subdomain_id=1)
# Initialize forms
a = Form([]) # Initialize empty bilinear form
L = Form([]) # Initialize empty linear form
## 1. Mass Terms
a += (1 / as_ufl(dt)) * inner(W, dot(A0_f, Y)) * dx
a -= (1 / as_ufl(dt)) * inner(W, dot(A0_f, Y)) * dx_s
form_a = fem.form(a)
form_L = fem.form(L)
# Assemble matrix and vector
A = fem.petsc.assemble_matrix(form_a)
A.assemble()
b = fem.petsc.assemble_vector(form_L)
file main_test.py
import dolfinx
from dolfinx.fem import functionspace, dirichletbc, locate_dofs_topological
from dolfinx.mesh import locate_entities_boundary
from petsc4py import PETSc
import numpy as np
from basix.ufl import element, mixed_element
from dolfinx import fem
from A0_f_test import create_coefficient_matrices_test
from assemble_test import assemble_weak_form_test
from mesh_test import create_nested_rectangles
# Time stepping parameters
num_steps = 100
dt = 1 / num_steps
max_iter = 10
tol = 1e-6
domain, markers = create_nested_rectangles()
# Physical constants
R = fem.Constant(domain, np.float64(287.05))
c_v = fem.Constant(domain, np.float64(718.0))
c_p = fem.Constant(domain, np.float64(1005.0))
mu = fem.Constant(domain, np.float64(1.81e-5))
kappa = fem.Constant(domain, np.float64(0.025))
# Create mesh
length = 0.4
height = 0.4
nx = 40
ny = 40
# Create mesh and markers using the previously defined function
# Define function spaces
# Create elements
P1_pressure = element("Lagrange", domain.basix_cell(), degree=1, dtype=dolfinx.default_real_type)
P1_velocity = element("Lagrange", domain.basix_cell(), degree=1, shape=(domain.geometry.dim,), dtype=dolfinx.default_real_type)
P1_temperature = element("Lagrange", domain.basix_cell(), degree=1, dtype=dolfinx.default_real_type)
#V_u, V_p, V_T = functionspace(domain, P1_velocity), functionspace(domain, P1_pressure),functionspace(domain, P1_temperature)
V_fluid = mixed_element([P1_pressure, P1_velocity, P1_temperature])
V_fluid_fs=functionspace(domain,V_fluid)
p_init = fem.Function(V_fluid_fs.sub(0).collapse()[0])
u_init = fem.Function(V_fluid_fs.sub(1).collapse()[0])
T_init = fem.Function(V_fluid_fs.sub(2).collapse()[0])
# Define expressions for each component
# Define initial condition parameters
center_x = 0.0
center_y = 0.2
radius = 0.06
high_pressure = 6746268.65
low_pressure = 100000.0
high_temp = 1465.0
low_temp = 290.0
def initial_condition_expr(x, high_val, low_val, center_x, center_y, radius):
values = np.zeros(x.shape[1])
r_squared = (x[0] - center_x)**2 + (x[1] - center_y)**2
values = np.where(r_squared < radius**2, high_val, low_val)
return values
# Create functions for individual fields
# Apply initial conditions to individual fields
p_init.interpolate(lambda x: initial_condition_expr(x, high_pressure, low_pressure,
center_x, center_y, radius))
T_init.interpolate(lambda x: initial_condition_expr(x, high_temp, low_temp,
center_x, center_y, radius))
u_init.interpolate(lambda x: np.zeros((2, x.shape[1])))
# Create the mixed function for the complete solution
solution_n = fem.Function(V_fluid_fs)
# Assign initial values to each component of the mixed function
solution_n.sub(0).interpolate(p_init) # Pressure component
solution_n.sub(1).interpolate(u_init) # Velocity component
solution_n.sub(2).interpolate(T_init)
def noslip_boundary(x):
return np.isclose(x[0], 0.0) | np.isclose(x[0], 0.4) | np.isclose(x[1], 0.0)| np.isclose(x[1], 0.4)
noslip = np.zeros(domain.geometry.dim, dtype=PETSc.ScalarType) # type: ignore
facets = locate_entities_boundary(domain, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V_fluid_fs.sub(1).collapse()[0], 1, facets), V_fluid_fs.sub(1).collapse()[0])
for n in range(num_steps):
t = n * dt
# Store current solution as previous iteration
solution_k = fem.Function(V_fluid_fs)
solution_k=solution_n
# Fixed-point iteration loop
for k in range(max_iter):
# Split current iteration solutions
p_k, u_k, T_k = solution_k.split()
# Update coefficient matrices for current iteration
A0_f= create_coefficient_matrices_test(solution_n.sub(0),solution_n.sub(1),solution_n.sub(2), domain)
a, L = assemble_weak_form_test(V_fluid_fs, dt, A0_f,markers)
run main_test.py then you can reproduce the error in the title