Hi, I tried to use the tutorial on Dirichlet, Neumann and Robin boundary conditions to solve a different pde but the precision is poor and I can’t understand why.
Thanks in advance for the help, it’s been two weeks since I started trying to solve this…
(If it can help, I think the problem is in the Robin bc because by substituting them with Dirichlet/Neumann it works)
The pde is the following:
For simplicity I used a non-homogeneous Robin bc but with the homogeneous condition the precision is still low.
My code (slight modification of the tutorial) is
import numpy as np
import pyvista
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction,
div, dot, dx, grad, inner, lhs, rhs)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))
import ufl
import dolfinx
from petsc4py.PETSc import ScalarType
####
n_temp = 1.4
reff= -1.44*n_temp**-2+0.71*n_temp**-1+0.668+0.0636*n_temp
temp = (1+reff)/(1-reff)
A = Constant(mesh, ScalarType(temp))
####
x = ufl.SpatialCoordinate(mesh)
u_ex = lambda x: x[0]**3*x[1]+x[0]**2*x[1]+2*x[1]**2+10*x[1]
n = ufl.FacetNormal(mesh)
g = -ufl.dot(ufl.grad(u_ex(x)),n)
mu_s = Constant(mesh, ScalarType(1.22))
D = Constant(mesh, ScalarType(1/(3*1.22)))
mu_a = Constant(mesh, ScalarType(0.011))
S = -ufl.div(D*ufl.grad(u_ex(x)))+mu_a*u_ex(x)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = D * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + mu_a * ufl.inner(u,v) * ufl.dx - ufl.inner(S, v) * ufl.dx r
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))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = np.array(facet_tag.indices[facet_tag.values == marker])
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
elif type == "Robin":
self._bc = values[0] * inner(u-values[1], v)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
rob = D*ufl.inner(ufl.grad(u_ex(x)),n)+u_ex(x)/(2*A)
boundary_conditions = [BoundaryCondition("Dirichlet", 1, u_ex),
BoundaryCondition("Dirichlet", 2, u_ex),
BoundaryCondition("Robin", 3, (1/(2*A*D), 2*A*rob)),
BoundaryCondition("Neumann", 4, g)]
bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# Compute L2 error and error at nodes
V_ex = FunctionSpace(mesh, ("CG", 2))
u_exact = Function(V_ex)
u_exact.interpolate(u_ex)
error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(form((uh - u_exact)**2 * dx)), op=MPI.SUM))
u_vertex_values = uh.x.array
uex_1 = Function(V)
uex_1.interpolate(u_ex)
u_ex_vertex_values = uex_1.x.array
error_max = np.max(np.abs(u_vertex_values - u_ex_vertex_values))
error_max = mesh.comm.allreduce(error_max, op=MPI.MAX)
print(f"Error_L2 : {error_L2:.2e}")
print(f"Error_max : {error_max:.2e}")
(By removing the Neumann bc, a homogeneous Neumann bc is implied.) A complete code to solve this problem follows.
import numpy as np
import pyvista
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction,
div, dot, dx, grad, inner, lhs, rhs)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 3))
import ufl
import dolfinx
from petsc4py.PETSc import ScalarType
####
n_temp = 1.4
reff= -1.44*n_temp**-2+0.71*n_temp**-1+0.668+0.0636*n_temp
temp = (1+reff)/(1-reff)
A = Constant(mesh, ScalarType(temp))
####
x = ufl.SpatialCoordinate(mesh)
u_ex = lambda x: x[0]**3*x[1]+x[0]**2*x[1]+2*x[1]**2+10*x[1]
mu_s = Constant(mesh, ScalarType(1.22))
D = Constant(mesh, ScalarType(1/(3*1.22)))
mu_a = Constant(mesh, ScalarType(0.011))
n = ufl.FacetNormal(mesh)
g = -ufl.dot(D*ufl.grad(u_ex(x)),n)
S = -ufl.div(D*ufl.grad(u_ex(x)))+mu_a*u_ex(x)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = D * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + mu_a * ufl.inner(u,v) * ufl.dx - ufl.inner(S, v) * ufl.dx
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))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = np.array(facet_tag.indices[facet_tag.values == marker])
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
elif type == "Robin":
self._bc = values[0] * inner(u-values[1], v)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
rob = 2*A*D*ufl.inner(ufl.grad(u_ex(x)),n) + u_ex(x)
boundary_conditions = [BoundaryCondition("Dirichlet", 1, lambda x: np.zeros(x[0].size)),
BoundaryCondition("Dirichlet", 2, lambda x: np.zeros(x[0].size)),
BoundaryCondition("Robin", 3, (1/(2*A), Constant(mesh, 0.0))),]
#BoundaryCondition("Neumann", 4, g)]
bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
uh.name = 'uh'
# Compute L2 error and error at nodes
V_ex = FunctionSpace(mesh, ("CG", 3))
u_exact = Function(V_ex)
u_exact.interpolate(u_ex)
u_exact.name = 'u_exact'
error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(form((uh - u_exact)**2 * dx)), op=MPI.SUM))
u_vertex_values = uh.x.array
uex_1 = Function(V)
uex_1.interpolate(u_ex)
u_ex_vertex_values = uex_1.x.array
error_max = np.max(np.abs(u_vertex_values - u_ex_vertex_values))
error_max = mesh.comm.allreduce(error_max, op=MPI.MAX)
print(f"Error_L2 : {error_L2:.2e}")
print(f"Error_max : {error_max:.2e}")
Of course, this means that error calculated at the end of the code is meaningless, since you are comparing the discrete solution of your problem to an unrelated manufactured solution:
Error_L2 : 7.45e+00
Error_max : 1.40e+01
2. Computing the manufactured solution
Your code appears to be structured to compute a discrete approximation to the manufactured solution
u_{exact}(x) = x^3 y + x^2 y + 2y^2 + 10y
similarly to the tutorial but with an altered PDE and a different manufactured solution.
If this is your goal, there are three errors to correct in your code:
The definition of g, which should be (you need to add a factor of D):
g = -ufl.dot(D*ufl.grad(u_ex(x)),n)
The definition of rob, which should be:
rob = 2*A*D*ufl.inner(ufl.grad(u_ex(x)),n) + u_ex(x)
The definition of the Robin boundary condition, which should be:
