I am writting the poisson code with variable co-efficient in fenicsx
a\Delta u=f(x)
where a=a(x) variable coefficient. Could you please suggest something
Can there be any simple case where a(x)= x?
In Expression this was working
a = Expression(x[0], degree=1)
It is resolved, right ?
No, this I was using in fenics. Now I am writting in fenicsx. So, this Expression is not working there.
Could you please help?
This is covered many times in the dolfinx demos and the dolfinx tutorial The FEniCSx tutorial — FEniCSx tutorial
A simple edit of https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_poisson.py adding the lines
A = x[0] + 1
a = inner(A * grad(u), grad(v)) * dx
fulfils your example.
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import ds, dx, grad, inner
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
cell_type=mesh.CellType.triangle)
V = fem.FunctionSpace(msh, ("Lagrange", 1))
facets = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
np.isclose(x[0], 2.0)))
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
g = ufl.sin(5 * x[0])
A = x[0] + 1
a = inner(A * grad(u), grad(v)) * dx
L = inner(f, v) * dx + inner(g, v) * ds
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
with io.XDMFFile(msh.comm, "out_poisson/poisson.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(uh)
try:
import pyvista
cells, types, x = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = uh.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
if pyvista.OFF_SCREEN:
pyvista.start_xvfb(wait=0.1)
plotter.screenshot("uh_poisson.png")
else:
plotter.show()
except ModuleNotFoundError:
print("'pyvista' is required to visualise the solution")
print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")
If you want something not so easily expressible in symbolic algebra you can interpolate that function in an appropriate finite element space, e.g.,
def my_interesting_function(x):
value = ... # Do interesting things
return value
V_coefficient = dolfinx.fem.FunctionSpace(mesh, ("DG", 1)) # Choose appropriate space
A = dolfinx.fem.Function(V_coefficient)
A.interpolate(my_interesting_function)
1 Like
I have tried both the ways that you suggest. This is the minimal 1D code I have written to compute L2 norm. For constant coefficient I am getting results. But not for variable. Could you please check?
from mpi4py import MPI
from dolfinx import mesh
import numpy as np
from math import *
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import (CellType, compute_midpoints, create_unit_cube,
create_unit_square, meshtags)
from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction, TestFunction,
div, dx, ds, grad, inner)
from dolfinx.fem import Function, FunctionSpace, assemble_scalar, form
def poisson(N):
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(xL, xR),nx = N )
hm = 1/N
from dolfinx.fem import FunctionSpace
V = FunctionSpace(msh, ("CG", 1))
from dolfinx import fem
f = fem.Function(V)
# f.interpolate(lambda x: (pi**2)*np.sin(pi*x[0]) ) #### with constant coefficient \nabla u
f.interpolate(lambda x: (x[0] + 1)* ((pi**2)*np.sin(pi*x[0]) ) ) #### with variable coefficient A*\nabla u
import numpy
# Create facet to cell connectivity required to determine boundary facets
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(msh.topology)
dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
import ufl
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
# A = 1 #### with constant coefficient
A = x[0] + 1 #### with variable coefficient
##### or
# def my_interesting_function(x):
# return x[0] + 1
# V_coefficient = dolfinx.fem.FunctionSpace(msh, ("DG", 1)) # Choose appropriate space
# A = dolfinx.fem.Function(V_coefficient)
# A.interpolate(my_interesting_function)
a = inner(A * grad(u), grad(v)) * dx #### with variable coefficient
L = f * v * ufl.dx
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
V2 = fem.FunctionSpace(msh, ("CG", 2))
u_exact = fem.Function(V2)
u_exact.interpolate(lambda x: np.sin(pi*x[0]) )
L2_error = fem.form(ufl.inner(uh - u_exact, uh - u_exact) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = numpy.sqrt(msh.comm.allreduce(error_local, op=MPI.SUM))
return error_L2, hm
if __name__ == '__main__':
from math import log as ln
# Run the convergence test
for N in [4, 8, 16, 24, 48]:
e, h = poisson(N)
if N != 4:
rate = ln(e/e_)/ln(h/h_)
print(rate)
e_, h_ = e, h
Thank you. It is working now with variable coefficient, but my doubt is in 2D coefficient A has be to be in matrix form. How to define that A as matrix here?
import ufl
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem import (Expression, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities_boundary
from ufl import SpatialCoordinate, TestFunction, TrialFunction, div, dot, dx, grad, inner
def u_ex(mod):
return lambda x: mod.cos(2*mod.pi*x[0])*mod.cos(2*mod.pi*x[1])
u_numpy = u_ex(np)
u_ufl = u_ex(ufl)
def solve_poisson(N=10, degree=1):
mesh = create_unit_square(MPI.COMM_WORLD, N, N)
x = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, ("CG", degree))
u = TrialFunction(V)
v = TestFunction(V)
A = x[0]
f = -div(A*grad(u_ufl(x)))
a = inner(A*grad(u), grad(v)) * dx
L = f * v * dx
u_bc = Function(V)
u_bc.interpolate(u_numpy)
facets = locate_entities_boundary(mesh, mesh.topology.dim -1, lambda x: np.full(x.shape[1], True))
dofs = locate_dofs_topological(V, mesh.topology.dim-1, facets)
bcs = [dirichletbc(u_bc, dofs)]
default_problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
return default_problem.solve(), u_ufl(x)
uh, u_ex = solve_poisson(10)
def error_L2(uh, u_ex, degree_raise=3):
# Create higher order function space
degree = uh.function_space.ufl_element().degree()
family = uh.function_space.ufl_element().family()
mesh = uh.function_space.mesh
W = FunctionSpace(mesh, (family, degree+ degree_raise))
# Interpolate approximate solution
u_W = Function(W)
u_W.interpolate(uh)
# Interpolate exact solution, special handling if exact solution
# is a ufl expression or a python lambda function
u_ex_W = Function(W)
if isinstance(u_ex, ufl.core.expr.Expr):
u_expr = Expression(u_ex, W.element.interpolation_points)
u_ex_W.interpolate(u_expr)
else:
u_ex_W.interpolate(u_ex)
# Compute the error in the higher order function space
e_W = Function(W)
e_W.x.array[:] = u_W.x.array - u_ex_W.x.array
# Integrate the error
error = form(ufl.inner(e_W, e_W) * ufl.dx)
error_local = assemble_scalar(error)
error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
return np.sqrt(error_global)
Ns = [4, 8, 16, 32, 64]
Es = np.zeros(len(Ns), dtype=PETSc.ScalarType)
hs = np.zeros(len(Ns), dtype=np.float64)
for i, N in enumerate(Ns):
uh, u_ex = solve_poisson(N, degree=1)
comm = uh.function_space.mesh.comm
# One can send in either u_numpy or u_ex
# For L2 error estimations it is reccommended to send in u_numpy
# as no JIT compilation is required
Es[i] = error_L2(uh, u_numpy)
hs[i] = 1./Ns[i]
if comm.rank == 0:
print(f"h: {hs[i]:.2e} Error: {Es[i]:.2e}")
rates = np.log(Es[1:]/Es[:-1])/np.log(hs[1:]/hs[:-1])
if comm.rank == 0:
print(f"Rates: {rates}")
As you are using dolfinx, consider for instance Speeding up time-dependent diffusion equation - #4 by dokken