Hello everyone,
I am trying to solve a simple problem
\begin{aligned}
-\nabla \cdot \sigma \nabla V &= f && \text{in } \Omega, \\
V &= 1 && \text{in } \partial \Omega_{\text{in}}, \\
V &= -1 && \text{in } \partial \Omega_{\text{out}}, \\
\sigma \nabla V \cdot \vec{n} &= 0 && \text{in } \partial \Omega / (\partial \Omega_{\text{in}} \cup \partial \Omega_{\text{out}}), \\
\{V\}_{\Gamma} &= 0 && \text{on } \Gamma(t), \\
\{\sigma \nabla V \cdot \vec{n}\}_{\Gamma} &= \vec{g} && \text{on } \Gamma(t).
\end{aligned}
where f and g are derived from the fact that I choose an analytical solution V = -cos(pi*y). Also I have different material property in the two domain (the conductivity changes) and for the analytical solution I choose the jump on the interface is not zero.
I tried to look at different posts, also tried different solution. But I did not managed to let the code work. For sure I am messing with this jump condition, since if I change the analytical solution to V = cos(2piy) and accordingly the boundary condition, in this case the solution on the interface is zero (sin(pi/2) = 0) and I get the right convergence.
So, can you help me in writing correctly this jump condition on the interface due to the material property? Thanks in advance. I attach also the non-working code:
import matplotlib.pyplot as plt
import pyvista
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, div, FacetNormal,VectorElement, inner
import numpy as np
from dolfinx import default_scalar_type
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, io, plot
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
def V_exact(mode):
#return lambda x: mode.atan(mode.pi * x[1])
# return lambda x: (10000 -8000 * x[1]) * x[1]
return lambda x: -mode.cos(mode.pi*x[1])
V_numpy = V_exact(np) # which will be used for interpolation
V_ufl = V_exact(ufl) # which will be used for defining the source term
def solve_poisson(N, degree=1):
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
[N, N], mesh.CellType.triangle)
Q = fem.FunctionSpace(domain, ("DG", 0))
def Omega_0(x):
return x[1] <= 0.50
def Omega_1(x):
return x[1] >= 0.50
def inerface(x):
return x[1] == 0.50
x = SpatialCoordinate(domain)
sigma = fem.Function(Q)
cells_0 = mesh.locate_entities(domain, domain.topology.dim, Omega_0)
cells_1 = mesh.locate_entities(domain, domain.topology.dim, Omega_1)
def anode_conductivity(T):
return 1. / (5.929e-5 - T * 1.235e-8)
sigma.x.array[cells_0] = np.full_like(cells_0, 210, dtype=default_scalar_type)
sigma.x.array[cells_1] = np.full_like(cells_1, anode_conductivity(800), dtype=default_scalar_type)
# sigma = fem.Constant(domain, PETSc.ScalarType(500))
f = div(-sigma * grad(V_ufl(x)))
W = fem.FunctionSpace(domain, ("Lagrange", 1))
V = TrialFunction(W)
csi = TestFunction(W)
a = dot(sigma * grad(V), grad(csi)) * dx
# Boundary conditions
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], 1)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], 1)),
(5, lambda x: np.isclose(x[1], 0.5))]
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, 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(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
L = f * csi * dx + ufl.pi*(sigma('+')*csi('+') - sigma('-')*csi('-'))*ds(5)
# Dirichlet condition
facets = facet_tag.find(3)
dofs = fem.locate_dofs_topological(W, fdim, facets)
facets2 = facet_tag.find(4)
dofs2 = fem.locate_dofs_topological(W, fdim, facets2)
BCs = [fem.dirichletbc(PETSc.ScalarType(-1), dofs, W), fem.dirichletbc(PETSc.ScalarType(1), dofs2, W)]
default_problem = fem.petsc.LinearProblem(a, L, bcs=BCs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
return default_problem.solve(), V_ufl(x)
def error_L2_func(Vh, V_ex, degree_raise=3):
# Create higher order function space
degree = 1 #Vh.function_space.ufl_element().degree
family = Vh.function_space.ufl_element().family_name
mesh = Vh.function_space.mesh
Q = fem.functionspace(mesh, (family, degree + degree_raise))
# Interpolate approximate solution
V_W = fem.Function(Q)
V_W.interpolate(Vh)
# Interpolate exact solution, special handling if exact solution
# is a ufl expression or a python lambda function
V_ex_W = fem.Function(Q)
if isinstance(V_ex, ufl.core.expr.Expr):
u_expr = fem.Expression(V_ex, Q.element.interpolation_points)
V_ex_W.interpolate(u_expr)
else:
V_ex_W.interpolate(V_ex)
# Compute the error in the higher order function space
e_W = fem.Function(Q)
e_W.x.array[:] = V_W.x.array - V_ex_W.x.array
# Integrate the error
error = fem.form(ufl.inner(e_W, e_W) * ufl.dx)
error_local = fem.assemble_scalar(error)
error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
return np.sqrt(error_global)
N = [40, 80, 160, 320, 640]
error_L2 = []
error_H1 = []
h = []
mpi_rank = 5
for i in range(len(N)):
Vh, Vex = solve_poisson(N[i])
comm = Vh.function_space.mesh.comm
error_L2 += [error_L2_func(Vh, V_numpy)]
eh = Vh - Vex
error_H10 = fem.form(dot(grad(eh), grad(eh)) * dx)
E_H10 = np.sqrt(comm.allreduce(fem.assemble_scalar(error_H10), op=MPI.SUM))
error_H1 += [E_H10]
h += [1. / N[i]]
if comm.rank == 0:
mpi_rank = comm.rank
print(f"h: {h[i]:.2e} Error L2: {error_L2[i]:.2e}")
print(f"h: {h[i]:.2e} Error H1: {error_H1[i]:.2e}")
if mpi_rank == 0:
plt.figure(figsize=(10, 6))
plt.loglog(N, error_L2, label='$L^2$ error')
plt.loglog(N, error_H1, label='$H^1$ error')
plt.loglog(N, h)
h_square = [x**2 for x in h]
plt.loglog(N, h_square)
plt.xlabel('N')
plt.ylabel('Error')
plt.legend()
plt.grid(True)
plt.show()