Pde results explode

[eugene13240]

Hi, I’m new to Fenics. I have the following equations
\Delta u(x) - \nabla U(x) \cdot \nabla u(x) = 0, x\in \Omega_{ \setminus A\cup B}
u(x) = 0, x\in \partial A, u(x) = 1, x\in \partial B
\nabla u(x) \vec n = 0, x \in \partial \Omega
\Omega = \{x,y:x^2+y^2\leq4\}, A = \{x,y:(x+1)^2+y^2\leq0.5^2\}
B = \{x,y:(x-1)^2+y^2\leq0.5^2\}

And I obtain the variation formulation like this:
0=\int_{\Omega}(\Delta u - \nabla U \cdot \nabla u)e^{-U}v dx = -\int_{\Omega} \nabla u \nabla v e^{-U} dx
After I setting up the problem, the result of u I obtained is quite exploded from what I expected(say u \leq 1 for all x,y)

Here is my code. Could anyone help me figure it out. Many thanks.

import gmsh
from mpi4py import MPI
from dolfinx.io import gmshio
from gmsh import *

import numpy as np
from dolfinx.io import XDMFFile, gmshio
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import ds, dx, grad, inner
import pyvista

from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, locate_dofs_geometrical,
                         locate_dofs_topological)

from dolfinx import default_scalar_type
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner

gmsh.initialize()
gmsh.clear()
markerId = 1
disk = gmsh.model.occ.addDisk(0, 0, 0, 2, 2)
hole1 = gmsh.model.occ.addDisk(-1, 0, 0, 0.5, 0.5)
hole2 = gmsh.model.occ.addDisk(1, 0, 0, 0.5, 0.5)
membrane1 = gmsh.model.occ.cut([(2, disk)], [(2, hole1), (2,hole2)])

gmsh.model.occ.synchronize()

volumes = gmsh.model.getEntities(dim=2)
gdim = 2
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], 11)

boundary = gmsh.model.getBoundary(membrane1[0], oriented=False)
boundary_ids = [b[1] for b in boundary]
gmsh.model.addPhysicalGroup(2, boundary_ids, tag=1)

meshSize = 0.05
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",meshSize)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",meshSize)
gmsh.model.mesh.generate(gdim)

gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)

def outer(x):
    return np.isclose(x[0]**2+x[1]**2, 4)
def inner1(x):
    return np.isclose((x[0]+1)**2+x[1]**2, 0.25)
def inner2(x):
    return np.isclose((x[0]-1)**2+(x[1])**2, 0.25)

boundary_facets_outer = locate_entities_boundary(domain, domain.topology.dim-1, outer)
boundary_facets_inner1 = locate_entities_boundary(domain, domain.topology.dim-1, inner1)
boundary_facets_inner2 = locate_entities_boundary(domain, domain.topology.dim-1, inner2)

boundary_dofs_outer = locate_dofs_topological(V, domain.topology.dim - 1,boundary_facets_outer)

boundary_dofs_inner1 = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_inner1)

boundary_dofs_inner2 = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_inner2)

bcs_inner1 = dirichletbc(default_scalar_type(0), boundary_dofs_inner1, V)

bcs_inner2 = dirichletbc(default_scalar_type(1), boundary_dofs_inner2, V)

bcs = [bcs_inner1, bcs_inner2]

V = FunctionSpace(domain, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(domain)
uD = fem.Function(V)
class MyExpression:
    def __init__(self):
        self.t = 0.0
    def eval(self, x):
        return np.full(x.shape[1], np.exp((-x[0]**4)/4 + (x[0]**2)/2 - (x[1]**2)/2))

g = MyExpression()
uD.interpolate(g.eval)

a = dot(grad(u), grad(v))*uD * dx

f = Constant(domain, default_scalar_type(0))
L = f * v * dx

problem = LinearProblem(a,  L,bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

pyvista_cells, cell_types, geometry = vtk_mesh(V)

grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry)
grid.point_data["u"] = uh.x.array 
grid.set_active_scalars("u")

plotter = pyvista.Plotter()
plotter.add_text("uh", position="upper_edge", font_size=14, color="black")
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
plotter.show()

[dokken]

Dear @eugene13240,
Please include all imports to ensure that the code is reproducible.
Secondly, I would check that the boundary conditions are applied correctly.

Print

@eugene13240

boundary_dofs_inner1

@eugene13240

boundary_dofs_inner2

to see if there are any dofs the condition is applied to.

It is also unclear to me why there is an e^{-U} in your weak form.

[eugene13240]

Dear @dokken,
I have included all the imports.
And the boundary conditions is seem to be applied.

e^{-U} applies because \nabla \cdot (e^{-U} \nabla u) = -\nabla U e^{-U}\cdot \nabla u + e^{-U} \Delta u = 0
I also tried variation \int_{\Omega}(v\nabla U + \nabla v)\cdot \nabla u dx with

v_x = (x[0]**4)/4 - (x[0]**2)/2 + (x[1]**2)/2
uD.interpolate(fem.Expression(v_x, V.element.interpolation_points()))
f = grad(v_x)
a = dot(grad(u), grad(v)) * dx + dot(v*f,grad(u)) * dx
f = Constant(domain, default_scalar_type(0))
L = f * v * dx

which gives inf for each entry of the result.

[eugene13240]

I didn’t quite understand the petsc_options, but it seems that by changing from

problem = LinearProblem(a,  L, bcs=bcs, petsc_options={"ksp_type": "preonly","pc_type": "lu"})
uh = problem.solve()

towards

problem = LinearProblem(a,  L, bcs=bcs)
uh = problem.solve()

The result is much more tolerant like this.

But it’s still not much as I expected. Is it possible to constraint the value of u \in [0,1]?

[eugene13240]

I found that the potential constraint problems/explode problems comes from the fact that the boundary conditions is not successfully enforced.

essentially when using the default solver:

problem = LinearProblem(a,  L, bcs=bcs)
uh = problem.solve()
result = uh.x.array
print(result[boundary_dofs_inner1])

print(result[boundary_dofs_inner2])

one gives

with petsc_options one has

where neither satisfies the boundary conditions(and explosion occurs for the whole results with petsc_options)

Is there a solid way to enforce the boundary conditions to be satisfied? Many thanks for anyone can help me out

Here attached the code for reproducibility

import gmsh
from mpi4py import MPI
from dolfinx.io import gmshio
from gmsh import *

import numpy as np
from dolfinx.io import XDMFFile, gmshio
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import ds, dx, grad, inner
import pyvista
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, locate_dofs_geometrical,
                         locate_dofs_topological)
from dolfinx.plot import vtk_mesh
from dolfinx import default_scalar_type
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner

gmsh.initialize()
gmsh.clear()
markerId = 1
disk = gmsh.model.occ.addDisk(0, 0, 0, 2, 2)
hole1 = gmsh.model.occ.addDisk(-1, 0, 0, 0.5, 0.5)

hole2 = gmsh.model.occ.addDisk(1, 0, 0, 0.5, 0.5)
membrane1 = gmsh.model.occ.cut([(2, disk)], [(2, hole1), (2,hole2)])

gmsh.model.occ.synchronize()

volumes = gmsh.model.getEntities(dim=2)
gdim = 2
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], 11)

boundary = gmsh.model.getBoundary(membrane1[0], oriented=False)
boundary_ids = [b[1] for b in boundary]
gmsh.model.addPhysicalGroup(2, boundary_ids, tag=1)

meshSize = 0.1
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",meshSize)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",meshSize)
gmsh.model.mesh.generate(gdim)

gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)

def outer(x):
    return np.isclose(x[0]**2+x[1]**2, 4)
def inner1(x):
    return np.isclose((x[0]+1)**2+x[1]**2, 0.25)
def inner2(x):
    return np.isclose((x[0]-1)**2+(x[1])**2, 0.25)

boundary_facets_outer = locate_entities_boundary(domain, domain.topology.dim-1, outer)
boundary_facets_inner1 = locate_entities_boundary(domain, domain.topology.dim-1, inner1)
boundary_facets_inner2 = locate_entities_boundary(domain, domain.topology.dim-1, inner2)

boundary_dofs_outer = locate_dofs_topological(V, domain.topology.dim - 1,boundary_facets_outer)
boundary_dofs_inner1 = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_inner1)
boundary_dofs_inner2 = locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets_inner2)

bcs_inner1 = dirichletbc(default_scalar_type(0), boundary_dofs_inner1, V)

bcs_inner2 = dirichletbc(default_scalar_type(1), boundary_dofs_inner2, V)

bcs = [bcs_inner1, bcs_inner2]

V = FunctionSpace(domain, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(domain)
uD = fem.Function(V)

#v_x = 0
#v_x = (x[0]**4)/4 - (x[0]**2)/2 + (x[1]**2)/2

#uD.interpolate(fem.Expression(v_x, V.element.interpolation_points()))
#g = grad(v_x)
#a = -0.5*dot(grad(u), grad(v)) * dx - dot(g, grad(u)) * v * dx
#a = -0.5*dot(grad(u), grad(v)) * dx 

class MyExpression:
    def __init__(self):
        self.t = 0.0
    def eval(self, x):
        return np.full(x.shape[1], np.exp(-(x[0]**4)/4 + (x[0]**2)/2 - (x[1]**2)/2))
g = MyExpression()
uD.interpolate(g.eval)
a = -0.5*dot(grad(u), grad(v))*uD * dx

f = Constant(domain, default_scalar_type(0))
L = f * v * dx


problem = LinearProblem(a,  L, bcs=bcs)
uh = problem.solve()
result = uh.x.array
print(result[boundary_dofs_inner1])
print(result[boundary_dofs_inner2])

"""
problem = LinearProblem(a,  L, bcs=bcs, petsc_options={"ksp_type": "preonly","pc_type": "lu"})
uh = problem.solve()
result = uh.x.array
print(result[boundary_dofs_inner1])
print(result[boundary_dofs_inner2])
"""