Adding unknown to the nonlinear system

variational formulation
1

I have a problem which involves nonlinear Neumann boundary condition q(u,k) controlled by a parameter k. I want to solve the system and find k for which the integral at the boundary reach a specific value A. I want to add the equation q(u,k)*ds - A = 0 to the system and k as an additional dof. Is it possible to do so?

Following is an example of the problem for a simple Poisson equation on a square. Changing k will result in monotonically varying values of the flux between 0 and 1. How can I adjust the code so it solves for k to get specified flux in addition to the dofs of uh?

import ufl
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver 

def q(u,k):
    return k*u**2

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

V = fem.functionspace(domain, ("Lagrange", 1))

u_D = fem.Function(V)
u_D.x.array[:] = 1

def left_bnd(x):
    return np.isclose(x[0], 0)
def right_bnd(x):
    return np.isclose(x[0], 1)

fdim = domain.topology.dim - 1
left_boundary_facets = mesh.locate_entities_boundary(domain, fdim, left_bnd)
left_dofs = fem.locate_dofs_topological(V, fdim, left_boundary_facets)
bcs = [fem.dirichletbc(default_scalar_type(1), left_dofs, V)]

facets = mesh.locate_entities(domain, fdim, right_bnd)
facet_markers = np.full_like(facets, 0)
facet_tag = mesh.meshtags(domain, fdim, facets, facet_markers)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

k = 10
uh = fem.Function(V)
v = ufl.TestFunction(V)
F = ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx + v*q(uh,k)*ds(0)
problem = NonlinearProblem(F, uh, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
idx, converged = solver.solve(uh)

flux = fem.assemble_scalar(fem.form(q(uh,k) * ds(0)))
print(flux)
2

I think using real spaces would be an option, see

3

Dear Dokken, it seems like this should be the way. I was able to prepare the code using nightly build, but I have a problem. If I use ufl.MixedFunctionSpace, how can I create a Function object? I tried:

V = fem.functionspace(domain, ("Lagrange", 1))
R = create_real_functionspace(domain)
ME = basix.ufl.mixed_element([V.ufl_element(), R.ufl_element()]) 
W = fem.functionspace(domain, ME)
U = fem.Function(W)

which works, but probably wrong (instead of 122 dofs I get 321 and the solver gives inf in first iteration) and

W = ufl.MixedFunctionSpace(V, R)
U = fem.Function(W)

which does not work.

The code that runs but diverges:

import ufl
from basix.ufl import mixed_element
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, default_scalar_type, log
from dolfinx.fem.petsc import NonlinearProblem,NewtonSolverNonlinearProblem
from dolfinx.nls.petsc import NewtonSolver 
from scifem import create_real_functionspace, assemble_scalar
from petsc4py import PETSc

def q(u,k):
    return k*u**2

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = fem.functionspace(domain, ("Lagrange", 1))
R = create_real_functionspace(domain)
ME = mixed_element([V.ufl_element(), R.ufl_element()]) 

# Probably wrong ...
W = fem.functionspace(domain, ME)
# ... but following does not allow to create Function object
# W = ufl.MixedFunctionSpace(V, R)
U = fem.Function(W)
c, r = ufl.split(U)
c_test, r_test = ufl.TestFunctions(W)
c_tag = 0
r_tag = 1
c_sub = W.sub(c_tag)
r_sub = W.sub(r_tag)

def left_bnd(x):
    return np.isclose(x[0], 0)
def right_bnd(x):
    return np.isclose(x[0], 1)

fdim = domain.topology.dim - 1
left_boundary_facets = mesh.locate_entities_boundary(domain, fdim, left_bnd)
left_dofs = fem.locate_dofs_topological(c_sub, fdim, left_boundary_facets)
bcs = [fem.dirichletbc(default_scalar_type(1), left_dofs, c_sub)]

facets = mesh.locate_entities(domain, fdim, right_bnd)
facet_markers = np.full_like(facets, 0)
facet_tag = mesh.meshtags(domain, fdim, facets, facet_markers)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

zero = fem.Constant(domain, default_scalar_type(0.0))

A00 = ufl.dot(ufl.grad(c), ufl.grad(c_test)) * ufl.dx + c_test*q(c,r)*ds(0)
A01 = ufl.inner(zero, c_test)* r * ufl.dx
A10 = ufl.inner(zero, r_test)* c * ufl.dx
A11 = r_test*q(c,r)*ds(0) 
L0 = ufl.inner(zero, c_test) * ufl.dx
L1 = r_test * 0.5 * ufl.dx # expected integral: mesh_area*0.5, with mesh_area=1 

A = A00+A01+A10+A11
L = L0+L1
F = A - L

problem = NewtonSolverNonlinearProblem(F, U, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
U.x.array[:]=1
solver.solve(U)

flux_n = fem.assemble_scalar(fem.form(q(c,r) * ds(0)))
print(flux_n)