Hi all,
The following problem is a 2d simplification of a field theory calculation.
Consider the following problem: calculate a field u(x,y) in \Omega=[0,1]\times[0,1] that has obeys Dirichelt BCs u_{bc}=1+2x^2+3y^2 while also making sure it obeys the integral constraint \int\limits_0^{1/2} {dy} u\left( {\frac{1}{2},y} \right) = 1.
The constraint can be written as (\delta is the delta function and \theta is the Heaviside step function):
\int\limits_0^{1/2} {dy} u\left( {\frac{1}{2},y} \right) = \int {d\Omega } \delta \left( {x - \frac{1}{2}} \right)\theta \left( {\frac{1}{2} - y} \right)u(x,y) \equiv \int {d\Omega } g(x,y)u(x,y) = 1
The action is (\lambda is a Lagrange Multiplier):
S = \int {d\Omega } {\left| {\nabla u(x,y)} \right|^2} - \int {d\Omega } f(x,y)u(x,y) - \lambda \left[ {\int\limits_0^{1/2} {dy} u\left( {\frac{1}{2},y} \right) - 1} \right] = \int {d\Omega } \left[ {{{\left| {\nabla u} \right|}^2} - fu - \lambda gu} \right]
Minimization towards the LM gives back the integral constraint while minimizing towards the field gives the PDE:
- {\nabla ^2}u(x,y) = f(x,y) + \lambda g(x,y).
To solve the problem, I suggest the following procedure:
- Find the variational form of the PDE. Say that the trial function is \phi_i while the test function is \sum\limits_j {{u_j}{\varphi _j}} where \phi are proper elements for this problem. The variational form, thus, is:
\sum\limits_j {{A_{ij}}{u_j}} = {F_i} + \lambda {G_i}
where {A_{ij}} = \int {d\Omega } \nabla {\varphi _i} \cdot \nabla {\varphi _j},~~~{F_i}=\int {d\Omega } f{\varphi _i},~~~{G_i}=\int {d\Omega } g{\varphi _i} - Do the same in the integral constraint: \sum\limits_j {{G_j}{u_j}} \equiv G \cdot u = 1
- Write a matrix system as follows:\left[ {\begin{array}{*{20}{c}} A&G\\ {{G^ \bot }}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} u\\ { - \lambda } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} F\\ { - 1} \end{array}} \right]
- Make the usual assumptions: \delta (x - a) \to \frac{1}{{\sigma \sqrt {2\pi } }}\exp \left[ { - \frac{{{{\left( {x - a} \right)}^2}}}{{2{\sigma ^2}}}} \right],\theta (x - a) \to \frac{1}{{1 + {e^{ - \frac{{x - a}}{\varepsilon }}}}}
- Solve the problem in FEniCSx.
However, I cannot implement the problem correctly in FEniCSx. A code that presents the above idea but doesn’t work and the error will be attached in the end.
Therefore I have the following questions for discussion. Firstly, are the above correct? How can I do a proper implementation in FEniCSx? Secondly, can you suggest another viable option? I see the discussions in other topics about null spaces but I not familiar with them and my problem is not singular.
Finally, can you think of a “better” way to simulate the delta and step functions? Numerical tests in simple problems show high dependance of their parameters in the computational grid.
Thanks!!
import dolfinx
import ufl
import numpy as np
from mpi4py import MPI
import petsc4py.PETSc
import dolfinx.fem.petsc
# Create a mesh for the unit square domain
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 50, 50, dolfinx.mesh.CellType.quadrilateral)
V = dolfinx.fem.FunctionSpace(domain, ("CG", 3))
uD = dolfinx.fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
spatial_dimension = domain.topology.dim
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, domain.topology.dim-1, boundary_facets)
dof_coordinates = V.tabulate_dof_coordinates()[boundary_dofs]
bc = dolfinx.fem.dirichletbc(uD, boundary_dofs)
# Define test and trial functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Define the source term f(x, y)
f = dolfinx.fem.Constant(domain, petsc4py.PETSc.ScalarType(-6)) # Define your source term here, e.g., ufl.Constant(0)
x = ufl.SpatialCoordinate(domain)
ic_term = dolfinx.fem.Constant(domain, petsc4py.PETSc.ScalarType(0))
sigma_d=10**(-3)
sigma_s=10**(-5)
mean_d=0.5
mean_s=0.5
spatial_dimension = 2
dx_temp = ufl.exp(-((x[0] - mean_d)**2) / (2 * sigma_d**2)) / (sigma_d * ufl.sqrt(2 * ufl.pi)) * ((1 / (1 + ufl.exp(-(mean_s - x[1]) / sigma_s)))-(1 / (1 + ufl.exp(-(mean_s - x[1]) / sigma_s))))
# Define and solve the variational problem without constraint
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem_unconstrained = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc])
u_unconstrained = problem_unconstrained.solve()
u_unconstrained_array=np.array(u_unconstrained.x.array)
print(u_unconstrained_array)
# Define the integral constraint as a linear form
g = dx_temp
constraint_form = ufl.inner(g, v) * ufl.dx # Define the domain of the line integral
constraint_form_transpose = ufl.inner(v, g) * ufl.dx
# Stiffness Matrix
a_ic = dolfinx.fem.form([[a, constraint_form], [constraint_form_transpose, None]])
# Load Vector
L_ic = dolfinx.fem.form([[L], [1]])
problem_constrained = dolfinx.fem.petsc.LinearProblem(a_ic, L_ic, bcs=[bc])
u_constrained = problem_constrained.solve()
u_constrained_array=np.array(u_constrained.x.array)
Lagrange_multiplier = u_constrained_array[-1]
print(Lagrange_multiplier)
Error:
File "/Tests.py", line 54, in <module>
problem_constrained = dolfinx.fem.petsc.LinearProblem(a_ic, L_ic, bcs=[bc])
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 589, in __init__
self._A = create_matrix(self._a)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 139, in create_matrix
return _cpp.fem.petsc.create_matrix(a._cpp_object)
AttributeError: 'list' object has no attribute '_cpp_object'
Exception ignored in: <function LinearProblem.__del__ at 0x7f4d7e899900>
Traceback (most recent call last):
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 627, in __del__
self._solver.destroy()
AttributeError: 'LinearProblem' object has no attribute '_solver'