Hello,
I am trying to use use fenics for multiphysics problem but i do not understand how to properly use submeshes. How do I define variational form that solves two problem simultaniously?
And perhaps how do I create NonlinearVariationalProblem and Jacobian having mixed function spaces (potentialy with different number of variables) in each submesh? I always get errors and I cannot find a consistent way of doing such stuff.
Are there any good tutorials that use only fenics? Should i use multiphemics instead?
Consider this MWE:
from mshr import *
from fenics import *
import matplotlib.pyplot as plt
tol = 1E-14
# define two parts of the complete-domain
omega1 = Rectangle(Point(0, 0), Point(0.5, 1))
omega2 = Rectangle(Point(0.5, 0), Point(1, 1))
domain = omega1 + omega2
# set left and right domains
domain.set_subdomain(1, omega1)
domain.set_subdomain(2, omega2)
# generate mesh fuction that marks domains
mesh = generate_mesh(domain, 32)
marker = MeshFunction("size_t", mesh, 2, mesh.domains())
# create submeshes
submesh1 = MeshView.create(marker, 1)
submesh2 = MeshView.create(marker, 2)
CG = FiniteElement("CG", mesh.ufl_cell(), 1) # element type
element_1 = MixedElement([CG, CG]) # mixed element for left domain
element_2 = MixedElement([CG, CG]) # mixed element for right
W1 = FunctionSpace(submesh1, element_1) # function space for left domain
W2 = FunctionSpace(submesh2, element_2) # function space for right domain
W = MixedFunctionSpace(W1, W2)
dx_1 = Measure("dx", domain=submesh1) # so we can integrate only over left domain
dx_2 = Measure("dx", domain=submesh2) # so we can integrate only over right domain
u = Function(W)
u1, u2 = u.split() # mixed solutions of each domain
T1, C1 = u1.split() # individual fields of the let domain
T2, C2 = u2.split() # individual fields of the roght
v1, v2 = TestFunctions(W)
T1_v, C1_v = v1
T2_v, C2_v = v2
lmbd_1 = 1
lmbd_2 = 2
# poisson for left side
PDE = inner(lmbd_1*nabla_grad(T1), nabla_grad(T1_v))*dx_1
# poisson for right side
PDE += inner(lmbd_2*nabla_grad(T2), nabla_grad(T2_v))*dx_2
# hello C1 be constant on the left
PDE += C1*C1_v*dx_1 - Constant(1.0)*C1_v*dx_1
# hello C2 be constant on the right but higher
PDE += C2*C2_v*dx_2 - Constant(2.0)*C2_v*dx_2
# poisson needs a boundary conditions
def boundary_left(x, on_boundary):
return on_boundary and near(x[0], 0., tol)
def boundary_mid(x, on_boundary):
# How does on_boundary behave with meshview?
return on_boundary and near(x[0], 0.5, tol)
def boundary_right(x, on_boundary):
return on_boundary and near(x[0], 1., tol)
# not sure about this
bc1 = DirichletBC(W.sub_space(0).sub(0), Constant(1.0), boundary_left)
bc2 = DirichletBC(W.sub_space(0).sub(0), Constant(2.0), boundary_mid)
bc3 = DirichletBC(W.sub_space(1).sub(0), Constant(1.0), boundary_mid)
bc4 = DirichletBC(W.sub_space(1).sub(0), Constant(2.0), boundary_right)
bcs = [bc1, bc2, bc3, bc4]
solve(PDE == 0, u, bcs) #KSP error internal 76
sol_u1, sol_u2 = u.split() # domain solutions
sol_T1, sol_C1 = sol_u1.split() # left side solution
c = plot(sol_T1)
plt.colorbar(c)
plt.show()
Which gives this error:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'KSPSolve'.
*** Reason: PETSc error code is: 76 (Error in external library).
*** Where: This error was encountered inside ./dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------
How do I make the MWE above work?
Thank you for any reply.