Poor precision on tutorial with different function

mathematics
2

Your code appears to compute a discrete approximation to the manufactured solution

u_{exact}(x) = x^3 y + x^2 y + 2y^2 + 10y

However, this solution is not the solution to the problem stated in your question:

I have provided two codes below: one to solve the problem stated in your question, and one to solve the manufactured solution present in your code.

1. Solving the problem stated in the question

If you wish to solve this problem:

you will need to modify your boundary conditions to have homogeneous Dirichlet, Neumann, and Robin bcs:

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)]

(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:

  1. 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)

  2. The definition of rob, which should be:

    rob = 2*A*D*ufl.inner(ufl.grad(u_ex(x)),n) + u_ex(x)

  3. The definition of the Robin boundary condition, which should be:

    boundary_conditions = [BoundaryCondition("Dirichlet", 1, u_ex),
                       BoundaryCondition("Dirichlet", 2, u_ex),
                       BoundaryCondition("Robin", 3, (1/(2*A), rob)),
                       BoundaryCondition("Neumann", 4, g)]

The complete code to replicate the manufactured solution is below:

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]

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, u_ex),
                       BoundaryCondition("Dirichlet", 2, u_ex),
                       BoundaryCondition("Robin", 3, (1/(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()
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}")

Producing:

Error_L2 : 8.65e-03
Error_max : 5.59e-03
4 Likes