Hello,
I am solving an unsteady, coupled multiphysics phase-transition problem with remeshing in every single time step. I have one problem and several questions, and if desired I can also split them into different posts. The MWE is faily exhaustive though because of the complicated setting. I tried minimizing it as far as possible.
I am working in the official docker container using dolfinx0.8.0.
The reduced setting of the MWE is the following:
My domain is the 2D unit square, vertically split (at x = 0.5). The left domain part is solid, the right one is gas. I am solving an unsteady mixed Poisson (heat transport) in both domains, weakly coupled at the vertical split (interface). I am setting essential BCs for the heat flux on the top walls (zero flux), and natural temperature conditions on the left and right side. I have two temperatures and two heat fluxes, both defined on the whole domain, but only relevant in the corresponding subdomains. I simulate the remeshing by changing the grid resolution in every time step. In my actual application the mesh/the vertical split is moving arbitrarily.
The equations I am solving in detail:
All constants from the heat equation are set to one and I am using a big time step so stability will not be an issue, also not for my choice of elements (first order for T and q).
The MWE:
# Physical Markers (Facets and Cells)
# (5) (4)
# -------------------------
# | | |
# | solid | gas |
# | | |
# (3)| |(1) |(6)
# | | |
# | | |
# | 2 | 1 |
# | | |
# -------------------------
# (9) (8)
import numpy as np
import scipy as sp
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, la, io, geometry
from dolfinx.fem import petsc
from dolfinx.io import gmshio
from basix.ufl import element
import ufl
import gmsh
import os
import shutil
import random
random.seed(123)
def create_mesh():
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # Disables Output
mesh_comm = MPI.COMM_WORLD
model_rank = 0 # Create mesh on rank 0
if mesh_comm.rank == model_rank:
gmsh.model.geo.addPoint(0.0, 0.0, 0.1, 1, tag=1)
gmsh.model.geo.addPoint(0.0, 1.0, 0.1, 1, tag=2)
gmsh.model.geo.addPoint(0.5, 1.0, 0.1, 1, tag=3)
gmsh.model.geo.addPoint(1.0, 1.0, 0.1, 1, tag=4)
gmsh.model.geo.addPoint(1.0, 0.0, 0.1, 1, tag=5)
gmsh.model.geo.addPoint(0.5, 0.0, 0.1, 1, tag=6)
gmsh.model.geo.addLine(1, 2, tag=1)
gmsh.model.geo.addLine(2, 3, tag=2)
gmsh.model.geo.addLine(3, 4, tag=3)
gmsh.model.geo.addLine(4, 5, tag=4)
gmsh.model.geo.addLine(5, 6, tag=5)
gmsh.model.geo.addLine(6, 1, tag=6)
gmsh.model.geo.addLine(3, 6, tag=7) # Seperation line
gmsh.model.geo.addCurveLoop([1, 2, 7, 6], tag=1)
gmsh.model.geo.addCurveLoop([3, 4, 5, -7], tag=2)
gmsh.model.geo.addPlaneSurface([1], tag=1)
gmsh.model.geo.addPlaneSurface([2], tag=2)
gmsh.model.geo.synchronize()
gmsh.option.setNumber("Mesh.Algorithm", 5)
gmsh.option.setNumber("Mesh.MinimumCurveNodes", 4)
mesh_size = random.uniform(0.04, 0.06)
gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_size)
gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_size)
gmsh.model.mesh.generate(2)
gmsh.model.addPhysicalGroup(dim=1, tags=[7], tag=1, name="interface") # the separation line marker
gmsh.model.addPhysicalGroup(dim=1, tags=[1], tag=3, name="left_solid")
gmsh.model.addPhysicalGroup(dim=1, tags=[3], tag=4, name="top_gas")
gmsh.model.addPhysicalGroup(dim=1, tags=[2], tag=5, name="top_solid")
gmsh.model.addPhysicalGroup(dim=1, tags=[4], tag=6, name="right_gas")
gmsh.model.addPhysicalGroup(dim=1, tags=[5], tag=8, name="bot_gas")
gmsh.model.addPhysicalGroup(dim=1, tags=[6], tag=9, name="bot_solid")
gmsh.model.addPhysicalGroup(dim=2, tags=[2], tag=1, name="gas")
gmsh.model.addPhysicalGroup(dim=2, tags=[1], tag=2, name="solid")
gmsh.model.mesh.setOrder(1)
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()
return domain, cell_markers, facet_markers
n_steps = 10
dt = 50
t = 0
# Initial condition
domain, cell_markers, facet_markers = create_mesh()
P1 = element("Lagrange", domain.topology.cell_name(), 1)
P1_vec = element("Lagrange", domain.topology.cell_name(), 1, shape=(domain.geometry.dim, ))
V_T1 = fem.functionspace(domain, P1)
V_q1 = fem.functionspace(domain, P1_vec)
V_T2 = fem.functionspace(domain, P1)
V_q2 = fem.functionspace(domain, P1_vec)
# Set temperature to 3.0 everywhere (initial condition)
T1_prev = fem.Function(V_T1)
T1_prev.x.array[:] = 3.0
T2_prev = fem.Function(V_T2)
T2_prev.x.array[:] = 3.0
# For plotting only
q1_prev = fem.Function(V_q1)
q1_prev.x.array[:] = 0.0
q2_prev = fem.Function(V_q2)
q2_prev.x.array[:] = 0.0
T1_prev.name = "T_right"
T2_prev.name = "T_left"
q1_prev.name = "q_right"
q2_prev.name = "q_left"
# Output
folder_path = 'out_MWE'
if os.path.exists(folder_path):
try:
shutil.rmtree(folder_path)
print(f"{folder_path} and its contents have been deleted.")
except OSError as e:
print(f"Error: {e}")
else:
print(f"{folder_path} does not exist.")
vtk_outfile = io.VTKFile(domain.comm, folder_path + '/mixed_poisson.pvd', "w")
vtk_outfile.write_mesh(domain, t)
vtk_outfile.write_function(T1_prev, t)
vtk_outfile.write_function(T2_prev, t)
vtk_outfile.write_function(q1_prev, t)
vtk_outfile.write_function(q2_prev, t)
for n in range(n_steps):
t += dt
T1 = ufl.TrialFunction(V_T1)
q1 = ufl.TrialFunction(V_q1)
T2 = ufl.TrialFunction(V_T2)
q2 = ufl.TrialFunction(V_q2)
Theta1 = ufl.TestFunction(V_T1)
p1 = ufl.TestFunction(V_q1)
Theta2 = ufl.TestFunction(V_T2)
p2 = ufl.TestFunction(V_q2)
dx = ufl.Measure("cell", domain, subdomain_data=cell_markers)
ds = ufl.Measure("exterior_facet", domain, subdomain_data=facet_markers)
dS = ufl.Measure("interior_facet", domain, subdomain_data=facet_markers)
n = ufl.FacetNormal(domain)
# Bilinear form
a_00 = 1.0 / dt * T1 * Theta1 * dx(1)
a_01 = ufl.div(q1) * Theta1 * dx(1)
a_10 = -T1 * ufl.div(p1) * dx(1)
a_11 = ufl.inner(q1, p1) * dx(1)
a_22 = 1.0 / dt * T2 * Theta2 * dx(2)
a_23 = ufl.div(q2) * Theta2 * dx(2)
a_34 = -T2 * ufl.div(p2) * dx(2)
a_33 = ufl.inner(q2, p2) * dx(2)
# Linear form (incl. natural temperature BCs for walls to left and right)
L_0 = 1.0 / dt * T1_prev * Theta1 * ds
L_1 = -4.0 * ufl.inner(p1, n) * ds(6)
L_2 = 1.0 / dt * T2_prev * Theta2 * ds
L_3 = -2.0 * ufl.inner(p2, n) * ds(3)
# Essential BCs for the heat flux, set y-component of top and bot walls to zero
essential_bcs = []
# Gas
Q, _ = V_q1.sub(1).collapse()
bc_fun = fem.Function(Q)
bc_fun.interpolate(lambda x: x[0] * 0.0)
dofs_top = fem.locate_dofs_topological((V_q1.sub(1), Q), domain.topology.dim - 1, facet_markers.find(4))
dofs_bot = fem.locate_dofs_topological((V_q1.sub(1), Q), domain.topology.dim - 1, facet_markers.find(8))
dofs = [np.hstack((dofs_top[0], dofs_bot[0])), np.hstack((dofs_top[1], dofs_bot[1]))]
essential_bcs.append(fem.dirichletbc(bc_fun, dofs, V_q1.sub(1)))
# Solid
Q, _ = V_q2.sub(1).collapse()
bc_fun = fem.Function(Q)
bc_fun.interpolate(lambda x: x[0] * 0.0)
dofs_top = fem.locate_dofs_topological((V_q2.sub(1), Q), domain.topology.dim - 1, facet_markers.find(5))
dofs_bot = fem.locate_dofs_topological((V_q2.sub(1), Q), domain.topology.dim - 1, facet_markers.find(9))
dofs = [np.hstack((dofs_top[0], dofs_bot[0])), np.hstack((dofs_top[1], dofs_bot[1]))]
essential_bcs.append(fem.dirichletbc(bc_fun, dofs, V_q2.sub(1)))
# 'Inner' conditions with the facet normal from ufl
# (gives strange artifacts, probably because it is not correct direction)
a_12 = T2('+') * ufl.inner(p1("-"), n("-")) * dS(1) # weak coupling term
i_3 = -3.0 * ufl.inner(p2("+"), n("-")) * dS(1)
# 'Known' outer normal for first domain (solid): (1, 0), for gas: (-1, 0) (gives correct results)
custom_normal = fem.Function(V_q1)
custom_normal.x.array[:] = 0.0
custom_normal.x.array[0::2] = 1.0
#a_12 = T2('+') * ufl.inner(p1("-"), custom_normal('-')*(-1.0)) * dS(1) # weak coupling term
#i_3 = -3.0 * ufl.inner(p2("+"), custom_normal('-')) * dS(1)
L = [L_0, L_1, L_2, L_3 + i_3]
a = [[a_00, a_01, None, None],
[a_10, a_11, a_12, None],
[None, None, a_22, a_23],
[None, None, a_34, a_33]]
# This matrix has zero rows: in the upper left block for subdomain 2, in the lower left one for subdomain 1. This is fixed with dirichlet-bcs:
a[0][0] += fem.Constant(domain, PETSc.ScalarType(0)) * T1 * Theta1 * dx(2)
a[1][1] += fem.Constant(domain, PETSc.ScalarType(0)) * ufl.inner(q1, p1) * dx(2)
a[2][2] += fem.Constant(domain, PETSc.ScalarType(0)) * T2 * Theta2 * dx(1)
a[3][3] += fem.Constant(domain, PETSc.ScalarType(0)) * ufl.inner(q2, p2) * dx(1)
# Get all DOFs that we want to remove and make them a Dirichlet value.
bcs_fleshout = []
interface_dofs = fem.locate_dofs_topological(V_T1, domain.topology.dim - 1, facet_markers.find(1))
domain.topology.create_connectivity(2, 2) # This is required since 0.8.0?
dofs_to_remove_from_1 = fem.locate_dofs_topological(V_T1, domain.topology.dim, cell_markers.find(2))
dofs_to_remove_from_1 = np.setdiff1d(dofs_to_remove_from_1, interface_dofs, assume_unique=True) # Don't remove interface
dofs_to_remove_from_2 = fem.locate_dofs_topological(V_T2, domain.topology.dim, cell_markers.find(1))
dofs_to_remove_from_2 = np.setdiff1d(dofs_to_remove_from_2, interface_dofs, assume_unique=True)
bcs_fleshout.append(fem.dirichletbc(fem.Constant(domain, PETSc.ScalarType(1)), dofs_to_remove_from_1, V_T1))
bcs_fleshout.append(fem.dirichletbc(np.array([fem.Constant(domain, PETSc.ScalarType(1)), fem.Constant(domain, PETSc.ScalarType(1))], dtype=PETSc.ScalarType), dofs_to_remove_from_1, V_q1))
bcs_fleshout.append(fem.dirichletbc(fem.Constant(domain, PETSc.ScalarType(1)), dofs_to_remove_from_2, V_T2))
bcs_fleshout.append(fem.dirichletbc(np.array([fem.Constant(domain, PETSc.ScalarType(1)), fem.Constant(domain, PETSc.ScalarType(1))], dtype=PETSc.ScalarType), dofs_to_remove_from_2, V_q2))
# Assemble matrix and vector
bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = petsc.assemble_matrix_nest(bilinear_form, bcs=bcs_fleshout + essential_bcs)
A.assemble()
b = petsc.assemble_vector_nest(linear_form)
petsc.apply_lifting_nest(b, bilinear_form, bcs=bcs_fleshout + essential_bcs)
for b_sub in b.getNestSubVecs():
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
# I do not care about the value of fleshed-out DOFs, so I only apply the essential_bcs to the RHS
bcs0 = fem.bcs_by_block(fem.extract_function_spaces(linear_form), essential_bcs)
petsc.set_bc_nest(b, bcs0)
ksp = PETSc.KSP()
ksp.create(domain.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.getPC().setFactorSetUpSolverType()
ksp.setFromOptions()
Th1 = fem.Function(V_T1)
qh1 = fem.Function(V_q1)
Th2 = fem.Function(V_T2)
qh2 = fem.Function(V_q2)
x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(Th1.x),
la.create_petsc_vector_wrap(qh1.x),
la.create_petsc_vector_wrap(Th2.x),
la.create_petsc_vector_wrap(qh2.x)])
ksp.solve(b, x)
Th1.name = "T_right"
Th2.name = "T_left"
qh1.name = "q_right"
qh2.name = "q_left"
vtk_outfile.write_mesh(domain, t)
vtk_outfile.write_function(Th1, t)
vtk_outfile.write_function(Th2, t)
vtk_outfile.write_function(qh1, t)
vtk_outfile.write_function(qh2, t)
# New domain for next time step
new_domain, new_cell_markers, new_facet_markers = create_mesh()
P1 = element("Lagrange", new_domain.topology.cell_name(), 1)
P1_vec = element("Lagrange", new_domain.topology.cell_name(), 1, shape=(new_domain.geometry.dim, ))
V_T1 = fem.functionspace(new_domain, P1)
V_q1 = fem.functionspace(new_domain, P1_vec)
V_T2 = fem.functionspace(new_domain, P1)
V_q2 = fem.functionspace(new_domain, P1_vec)
# Interpolate solution to new domain for next time step (domain changes every timestep!)
T1_prev = fem.Function(V_T1)
interpolation_data = fem.create_nonmatching_meshes_interpolation_data(T1_prev.function_space.mesh,
T1_prev.function_space.element,
Th1.function_space.mesh,
padding=1.0e-10)
T1_prev.interpolate(Th1, nmm_interpolation_data=interpolation_data)
T2_prev = fem.Function(V_T2)
interpolation_data = fem.create_nonmatching_meshes_interpolation_data(T2_prev.function_space.mesh,
T2_prev.function_space.element,
Th2.function_space.mesh,
padding=1.0e-10)
T2_prev.interpolate(Th2, nmm_interpolation_data=interpolation_data)
domain = new_domain
cell_markers = new_cell_markers
facet_markers = new_facet_markers
vtk_outfile.close()
Questions:
(I) Interior Facet Orientation
The interior coupling conditions
a_12 = T2('+') * ufl.inner(p1("-"), n("-")) * dS(1) # weak coupling term
i_3 = -3.0 * ufl.inner(p2("+"), n("-")) * dS(1)
produce strange results for some time steps. I guess the orientation is wrong/inconsistent (it changes because the mesh changes). I also experimented with the ‘+’ and ‘-’ operators and tried every single combination, but nothing was successful. When I am using the correct inner facet normal (only known because this example is easy), I get the expected results. You can try this out, the relevant code is commented in the MWE. I read about this problem already here and also had a look at the pull request that made into fenicsx, but could not find a way to finally use it for my problem for general connected subsets of interior facets. Any idea? Or is the only chance to calculate the normal myself for each relevant interior facet?
(II) Essential Flux Conditions
I am setting the zero-flux condition by forcing the y-direction of the heat flux to zero at the boundary on top and bot. This is possible here, but when my boundary is not longer aligned with the coordinate system, what is the way to go? MPC or is there something in plain fenicsx?
(III) Remove Zero-Rows from Nested System Matrix
I have zero-rows in my system (because my equations are only defined on a subdomain, but the unknowns are defined everywhere). Currently, I am adding a zero-form to the bilinear form on the diagonal. This way I can set dirichletbc and make the diagonal zero where I want. In another project of mine I directly accessed the sparsity pattern and allocated what I needed, see here. This way I did not need to assemble the matrix where the values are irrelevant. I cannot access the sparsity pattern of the nested matrix directly. Is there an equivalently efficient way using the nested structure? I am aware of multiphenicsx but would like to stay with plain fenicsx.
Thank you for your time. If anything is unclear, please let me know.