Thanks for the changelog but I couldn’t use it to figure out my issue. It’s getting recreated on all my scripts, so I’ve linked a M(ish)WE, that I hope helps locate the error. It seems to stem from the solver now expecting a sequence of Forms, rather than just a Form. But I can’t find much on creating Form sequences.
from scipy.interpolate import RegularGridInterpolator as rgi
import ufl
import numpy as np
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, io
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
import dolfinx as df
k = 100
box_dim = [0.1, 0.1, 0.01]
box_res = 20
cube = mesh.create_box(MPI.COMM_WORLD, [[0.,0.,0.], box_dim], [box_res]*3, mesh.CellType.hexahedron)
tdim = cube.topology.dim
fdim = tdim - 1
cube.topology.create_connectivity(fdim, tdim)
# PARAMETERS
t_0 = fem.Constant(cube, 10.0)
V = fem.functionspace(cube, ("Lagrange", 1))
def left(x):
return np.isclose(x[2], 0)
def right(x):
return np.isclose(x[2], box_dim[2])
boundary_facets = mesh.exterior_facet_indices(cube.topology)
left_facets = mesh.locate_entities_boundary(cube, 2, left)
right_facets = mesh.locate_entities_boundary(cube, 2, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tags = mesh.meshtags(cube, 2, marked_facets[sorted_facets], marked_values[sorted_facets])
left_dbc_value = 0
left_dbc = np.array(left_dbc_value, dtype=PETSc.ScalarType)
left_dofs = fem.locate_dofs_topological(V, 2, facet_tags.find(1))
right_dbc_value = 0
right_dbc = np.array(right_dbc_value, dtype=PETSc.ScalarType)
right_dofs = fem.locate_dofs_topological(V, 2, facet_tags.find(2))
bc = []
#bc.append(fem.dirichletbc(left_dbc, left_dofs, V))
bc.append(fem.dirichletbc(right_dbc, right_dofs, V))
# Create initial condition
def initial_condition(x):
return np.zeros_like(x).flatten()
def heat_load(coord):
rangex = np.linspace(0, box_dim[0], 200)
rangey = np.linspace(0, box_dim[1], 200)
rangez = np.linspace(0, box_dim[2], 2)
heatload = np.ones(shape=(2, 200, 200)) * 20
heat_interpolator = rgi([rangez, rangey, rangex], heatload, "linear", False, 99)
out_array = []
for x, y, z in coord.T:
a = heat_interpolator([x, y, z])[0]
out_array.append(a)
return np.array(out_array)
heat_function = fem.Function(V, dtype=df.default_scalar_type)
heat_function.interpolate(heat_load)
u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)
uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], box_dim[0])),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], box_dim[1])),
(5, lambda x: np.isclose(x[2], 0)),
(6, lambda x: np.isclose(x[2], box_dim[2])),]
facet_indices, facet_markers = [], []
fdim = cube.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(cube, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(cube, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure('ds', domain=cube, subdomain_data=facet_tags, metadata={"quadrature_degree":4})
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
x = ufl.SpatialCoordinate(cube)
g = -4 * x[0]
f = fem.Constant(cube, PETSc.ScalarType(-6))
neuman_bc = g * v * ufl.ds
t_surr = fem.Constant(cube, PETSc.ScalarType(10))
r = 100 * x[0]
tdiff = u - t_surr
robin_coeff = r*tdiff
robin_bc = robin_coeff * v * ds(5)
heat_input = heat_function * v * ufl.dx
F = (
u * v * ufl.dx + k * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
- heat_input
)
a, L = ufl.lhs(F), ufl.rhs(F)
bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
b = create_vector(linear_form)
solver = PETSc.KSP().create(cube.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, linear_form)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [bc])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bc)
# Solve linear problem
solver.solve(b, uh.x.petsc_vec)
uh.x.scatter_forward()
Thanks in advanced!