Dear Community,
i am trying to implement the following coupled, scalar, “directed” (by b\in\mathbb{R}^d) advection/transport system with Dirichlet BC on both unknowns:
\begin{aligned}
\partial_t V +b\cdot \nabla \phi &= 0\\
\partial_t \phi + b\cdot \nabla V &= 0
\end{aligned}
The weak form for the problem on H^1\times H^1 is
\begin{aligned}
\left \langle \partial_t V ,W \right \rangle- \left \langle \phi b, \nabla W \right \rangle
&= - \left ( \phi_\Gamma b\cdot e_n , W \right )_\Gamma \quad \forall W\in H^1\\
\left \langle \partial_t \phi, \psi \right \rangle - \left \langle Vb, \nabla \psi \right \rangle
&= - \left ( V_\Gamma b\cdot e_n, \psi \right )_\Gamma \quad \forall \psi\in H^1
\end{aligned}
Or fully discrete, with implicit Euler:
\begin{aligned}
\begin{pmatrix}
M &
E \\
F &
N
\end{pmatrix}
\begin{pmatrix}
\vec V ^{n+1}\\
\vec\phi^{n+1}
\end{pmatrix} =
\begin{pmatrix}
M \vec V^n - \vec r^{n+1} \\
F\vec V^n -\vec s^{n+1}
\end{pmatrix}
\end{aligned}
with:
\begin{aligned}
M_{ij} &= \left \langle W_j, W_i \right \rangle , \qquad &&N_{ij} = \left \langle \psi_j,\psi_i \right \rangle \\
E_{ij} &= \left \langle \psi_j b, \nabla W_i \right \rangle , \qquad &&F_{ij} = \left \langle W_j b,\nabla \psi_i\right \rangle \\
r_i^{n+1} &= \left ( \phi_{\Gamma,j}^{n+1} b \cdot e_n, W_i \right )_\Gamma , \qquad &&
s_i^{n+1} = \left ( V_{\Gamma,j}^{n+1} b\cdot e_n, \psi_i \right )_\Gamma
\end{aligned}
I am following the Navier-Stokes example, BUT: I use MixedElement() in the following!
import numpy as np
import ufl
from basix.ufl import element
from dolfinx import plot
from dolfinx.fem import (Constant, Function, dirichletbc, form,
functionspace, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block,assemble_vector, create_vector_block,assemble_vector_block, create_vector
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.mesh import *
from ufl import div, dx, ds, grad, inner, FacetNormal, dot,
import pyvista
from mpi4py import MPI
from petsc4py import PETSc
# parameters
dt = 0.1
t = 0
T_end = 1*dt
num_steps = int(T_end/dt)
# mesh, spaces, functions
msh = create_unit_square(MPI.COMM_WORLD, 20,20)
x = ufl.SpatialCoordinate(msh)
n = FacetNormal(msh)
P1_elem = element("Lagrange", "triangle", 1)
P1_mixed = ufl.MixedElement(P1_elem,P1_elem)
X = functionspace(msh, P1_mixed) # mixed space P1 x P1
XV = X.sub(0).collapse()[0]
Xp = X.sub(1).collapse()[0]
W,q = ufl.TestFunctions(X)
V,p = ufl.TrialFunctions(X)
V_old = Function(XV)
V_new = Function(XV)
p_old = Function(Xp)
p_new = Function(Xp)
bfield = Constant(msh, (+1.0,0.0))
# init cond
def initcond_V(x):
return (x[0]-0.6)**2 + (x[1]-0.3)**2 < 0.1**2
def initcond_p(x):
return (x[0]-0.2)**2 + (x[1]-0.7)**2 < 0.1**2
V_old.interpolate(initcond_V)
p_old.interpolate(initcond_p)
# BC
facets = locate_entities_boundary(msh, dim=1,marker=lambda x: np.logical_or.reduce((
np.isclose(x[0], 0.0),
np.isclose(x[0], 1.0),
np.isclose(x[1], 0.0),
np.isclose(x[1], 1.0))))
V_dofs = locate_dofs_topological((X.sub(0), XV), 1, facets)
p_dofs = locate_dofs_topological((X.sub(1), Xp), 1, facets)
def inflow_BC_expression(x):
return np.zeros(X.dofmap.index_map.size_global)
BC_V = Function(XV)
BC_V.interpolate(inflow_BC_expression)
BC_p = Function(Xp)
BC_p.interpolate(inflow_BC_expression)
BCs = [dirichletbc(BC_V, V_dofs), dirichletbc(BC_p, p_dofs)] #?
# weak form implicit euler
M = inner(V,W)*dx
E = inner(bfield*p,grad(W))*dx
F = inner(V*bfield,grad(q))*dx
N = inner(p,q)*dx
a = form([[M, -dt/2*E],
[-dt/2*F, N]])
# assembly
A = assemble_matrix_block(a, bcs=BCs)
A.assemble()
Vp_vec = A.createVecRight()
# solver
solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
for i in range(num_steps):
t += dt
# Update the right hand side
l = form([inner(grad(V_old), grad(W)) * dx,
inner(grad(p_old), grad(q)) * dx])
L = assemble_vector_block(l,a,bcs=BCs)
# Solve linear problem
solver.solve(L, Vp_vec)
# extract solution
offset = XV.dofmap.index_map.size_global
V_new.x.array[:] = Vp_vec.array_r[:offset]
V_new.x.scatter_forward()
p_new.x.array[:(len(Vp_vec.array_r) - offset)] = Vp_vec.array_r[offset:]
p_new.x.scatter_forward()
# Update solution at previous time step
V_old.x.array = V_new.x.array
p_old.x.array = p_new.x.array
However, it outputs the following error for the “apply lifting”-command, (after tracing back to dolfinx/fem/assemble.py):
Traceback (most recent call last):
File ".../transport-system-mixed.py", line 157, in <module>
p_new.x.array[:] = Vp_vec.array_r[offset:]
~~~~~~~~~~~~~^^^
ValueError: could not broadcast input array from shape (1323,) into shape (441,)
I have questions concerning the following:
- Why does
createVecRight()create an array that is much too large? (also, shouldnt it be “createVecLeft()”?) - How do I assemble the matrix correctly by taking into account Dirichlet BC on both unknowns? (I think I cant just assemble the submatrices individually. E.g. assembling submatrix M with BC would not enforce the essential BC on the essential DOFs correctly, as there is submatrix E messing up the respective equations in the resulting system.) I assume the command
BCs = [dirichletbc(BC_V, V_dofs), dirichletbc(BC_p, p_dofs)]is wrong. Can i somehow concatenate the Dirichlet values before? - When do i need
apply_lifting(L, [a], [BCs])orset_bc(L,[BCs])on the right hand side vector?
Thank you very much for your answers!
Best regards