Hi, I would like advice for the community one more time regarding integral constraints.
This is a follow up to this link while using info for the integrations needed from this topic
The problem is the following: consider the 2D problem of the known FEniCSx tutorial:
- {\nabla ^2}u(x,y) = f(x,y) where x \in [0.1] and y \in [0.1] with \Omega = [0,1] \times [0,1], f(x,y)=-6 and Dirichlet BCs: u_{D}(x,y)=1+x^2+2y^2. The analytical solution is u(x,y)=1+x^2+2y^2.
Say that this problem is accopanied with the following integral condition: C=\int\limits_0^{0.5} {dy} u(0.5,y) = \frac{17}{24} = 0.708333 which is the actual result if we integrate the analytical solution. Due to the 2D problem being more complex I could not form a less trivial toy example. If one solves the problem, therefore, the aforementioned result should be recovered.
So, we introduce a Lagrange multiplier \lambda. The action is J[u,\lambda ] = \int\limits_0^1 {dx} \int\limits_0^1 {dy} {\left( {\vec \nabla u(x,y)} \right)^2} - \lambda \left( {\int\limits_0^{0.5} {dy} u(0.5,y) - \frac{17}{24} } \right) or J[u,\lambda] = \int\limits_0^1 {dx} \int\limits_0^1 {dy} \left[ {{{\left( {\vec \nabla u(x,y)} \right)}^2} - \lambda \left( {\delta (x)\theta (0.5 - y)u(x) - \frac{17}{24} } \right)} \right] and the new PDE is - {\nabla ^2}u(x,y) = f(x,y) + \lambda \delta (x)\theta (0.5 - y).
So, in the same sense as the first topic I tried to use dolfinx_mpc to concetrate the dofs in order for \lambda to be constant in the subdomain {x=0.5, y \in [0,0.5]}. To define the integral measurement I used the info from the second link above.
Consider v and u the test and trial functions for u(x,y) and r and m the test and trial functions for \lambda. If I am not mistaken the variational form is \int\limits_0^1 {dx} \int\limits_0^1 {dy} \vec \nabla v \cdot \vec \nabla u - \int\limits_0^{0.5} {dy} mv - \int\limits_0^{0.5} {dy} ur = \int\limits_0^1 {dx} \int\limits_0^1 {dy} fv + \int\limits_0^{0.5} {dy} Cr
Unfortunately, the code does not work. So, I would like to discuss two topics. Firstly, are the above sound from the mathematical point of view and secondly what is the wrong in my code. I think it is either the implementation of dolfinx_mpc or the definition of the integral measures dx. Of course, I am open to any suggestion on how to tackle this problem.
Any input is highly appreciated.
The minimal code and the error follow.
Code:
from mpi4py import MPI
import dolfinx
import dolfinx_mpc
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
comm = MPI.COMM_WORLD
# Analytical Solution
def uu(x): return (1 + x[0]**2 + 2*x[1]**2)
# Analytical Solution Vector
def u_v(x):
xT = x.T
n=x.shape[1]
values = np.zeros(n)
for i in range(n):
xx=xT[i]
values[i] = uu(xx)
return values
### Create mesh and function spaces
domain = dolfinx.mesh.create_unit_square(comm, 10 ,10)
element1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 2)
V = dolfinx.fem.FunctionSpace(domain, ufl.MixedElement([element1, element1]))
V0c, submap = V.sub(0).collapse()
V1c, submap = V.sub(1).collapse()
V0 = V.sub(0)
V1 = V.sub(1)
### Create test and test functions
v, r = ufl.TestFunctions(V)
uh = dolfinx.fem.Function(V)
u, m = ufl.TrialFunctions(V)#
### Mark boundaries
def boundary(x):
return np.logical_or(np.logical_or(np.logical_or(
np.isclose(x[0],0), np.isclose(x[0],1)),
np.isclose(x[1],0)),np.isclose(x[1],1))
# Set BC
u_bc = dolfinx.fem.Function(V0c)
u_bc.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim-1, boundary)
bc = dolfinx.fem.dirichletbc(
u_bc, dolfinx.fem.locate_dofs_topological((V.sub(0), V0c), domain.topology.dim-1, facets), V.sub(0))
### MultiPointConstraint
mpc = dolfinx_mpc.MultiPointConstraint(V)
# Find the master dof ((x,y) = (0,0)) and owning process
nodes = dolfinx.mesh.locate_entities_boundary(domain, 0, lambda x: np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0)))
master_dof = dolfinx.fem.locate_dofs_topological(V1, 0, nodes, remote=False)
master_owner = np.full_like(master_dof, comm.rank)
master_dof = np.concatenate(comm.allgather(master_dof), dtype=np.int64)
master_owner = np.concatenate(comm.allgather(master_owner), dtype=np.int32)
# Check for correct number of master dofs (should be 1)
if master_dof.size != 1:
raise RuntimeError(f"Expected to find 1 master dof; found {master_dof.size}")
else:
print(f"[{comm.rank}]: master_dof = {master_dof}")
print(f"[{comm.rank}]: master_owner = {master_owner}")
# Find the slave nodes
nodes = dolfinx.mesh.locate_entities(domain, 0, lambda x: ~np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0)))
slave_dof = dolfinx.fem.locate_dofs_topological(V1, 0, nodes, remote=False)
# Create arrays for use by mpc.add_constraint
masters = np.broadcast_to(master_dof, slave_dof.shape)
owners = np.broadcast_to(master_owner, slave_dof.shape)
coeffs = np.ones_like(masters, dtype=ScalarType)
offsets = np.arange(masters.size + 1, dtype=np.int32)
# Create constraint and finalize mpc
mpc.add_constraint(V, slave_dof, masters, coeffs, owners, offsets)
mpc.finalize()
# Create custom integral measure for Lagrange multiplier
# 1D Integral i.e. Integrate[1+x**2+2y**2,{y,0,0.5} while x=0.5]
def facet_eq_1D(x):
tol=1e-6
close_to_half = np.isclose(x[0], 0.5)
less_than_half = x[1] < 0.5 + tol
return np.logical_and(close_to_half, less_than_half)
facets = dolfinx.mesh.locate_entities(domain, domain.topology.dim-1, facet_eq_1D)
facet_tags = dolfinx.mesh.meshtags(domain, domain.topology.dim-1, facets, np.full_like(facets, 1))
# Define the measure
dx = ufl.Measure("dS", subdomain_data=facet_tags, subdomain_id=1 )
# Check if the integral is evaluated correctly
u_test = dolfinx.fem.Function(V0c)
u_test.interpolate(lambda x: 1 + x[0]**2 + 2*x[1]**2)
intergal = dolfinx.fem.form(u_test * dx)
print(dolfinx.fem.assemble_scalar(intergal)) # Should print Integrate[1+x**2+2y**3,{y,0,0.5} while x=0.5]=0.708333
### Create variational form and solve
C0 = dolfinx.fem.Constant(domain, ScalarType(0.708333))
# Source
f = dolfinx.fem.Constant(domain, ScalarType(-6))
a = (
ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx
- ufl.inner(m, v)*dx(1)
+ ufl.inner(u, r)*dx(1)
)
L = ufl.inner(f, v)*ufl.dx + ufl.inner(C0, r)*dx(1)
problem = dolfinx_mpc.LinearProblem(
a, L, mpc, bcs=[bc],#l, bc_r],
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
}
)
uh = problem.solve()
uh1=uh.x.array
new_array1 = [uh1[i] for i in range(len(uh1)) if i not in slave_dof]
new_array = [new_array1[i] for i in range(len(new_array1)) if i not in master_dof]
points = domain.geometry.x
nn=len(points)
for i in range (nn):
xx=points[i]
print(points[i],new_array[i],uu(points[i]))
Error:
[0]: master_dof = [437]
[0]: master_owner = [0]
0.7083333333333333
Discontinuous type Argument must be restricted.
ERROR:UFL:Discontinuous type Argument must be restricted.
Traceback (most recent call last):
File "/home/george/codes/Integral_Constraint_Local_2D_V2_Copy_Copy_Copy.py", line 111, in <module>
problem = dolfinx_mpc.LinearProblem(
File "/usr/lib/python3/dist-packages/dolfinx_mpc/problem.py", line 75, in __init__
self._a = _fem.form(a, jit_options=jit_options, form_compiler_options=form_compiler_options)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 176, in form
return _create_form(form)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 171, in _create_form
return _form(form)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 145, in _form
ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py", line 204, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 186, in compile_forms
impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 250, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
File "/usr/lib/python3/dist-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, options)
File "/usr/lib/python3/dist-packages/ffcx/analysis.py", line 86, in analyze_ufl_objects
form_data = tuple(_analyze_form(form, options) for form in forms)
File "/usr/lib/python3/dist-packages/ffcx/analysis.py", line 86, in <genexpr>
form_data = tuple(_analyze_form(form, options) for form in forms)
File "/usr/lib/python3/dist-packages/ffcx/analysis.py", line 164, in _analyze_form
form_data = ufl.algorithms.compute_form_data(
File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 321, in compute_form_data
form = apply_restrictions(form)
File "/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py", line 166, in apply_restrictions
return map_integrand_dags(rules, expression,
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 66, in map_integrand_dags
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 29, in map_integrands
mapped_integrals = [map_integrands(function, itg, only_integral_type)
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 29, in <listcomp>
mapped_integrals = [map_integrands(function, itg, only_integral_type)
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 37, in map_integrands
return itg.reconstruct(function(itg.integrand()))
File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 66, in <lambda>
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 36, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress,
File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 97, in map_expr_dags
r = handlers[v._ufl_typecode_](v)
File "/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py", line 94, in reference_value
g = self(f)
File "/usr/lib/python3/dist-packages/ufl/corealg/multifunction.py", line 97, in __call__
return self._handlers[o._ufl_typecode_](o, *args)
File "/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py", line 53, in _require_restriction
error("Discontinuous type %s must be restricted." % o._ufl_class_.__name__)
File "/usr/lib/python3/dist-packages/ufl/log.py", line 135, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Discontinuous type Argument must be restricted.```