jw3126
July 31, 2020, 10:00am
1
I would like to solve the following:
∇(u) = 0
∫u = Q*vol
Solution:
u = Q
Here is what I tried:
from fenics import *
mesh = UnitIntervalMesh(10)
vol = assemble(Constant(1)*dx(domain=mesh))
Q = 1
V = FunctionSpace(mesh, "P", 1)
v = TestFunction(V)
u = Function(V)
F = inner(grad(u)[0],v) * dx + (u - Q/vol) * dx
solve(F == 0, u)
However I get an error:
ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_0',).
How to formulate this problem correctly?
The volume of a mesh can be computed by
vol = assemble(Constant(1) * dx)
1 Like
jw3126
July 31, 2020, 11:13am
3
Thanks! I edited this into my original post.
dokken
July 31, 2020, 12:14pm
4
I would suggest having a look at this Poisson with only Neumann conditions demo
Which shows you how to use lagrange multipliers in the case Q=0. It is very easy to add a non-zero Q to this constraint
1 Like
jw3126
July 31, 2020, 2:19pm
6
from fenics import *
mesh = UnitIntervalMesh(10)
vol = assemble(Constant(1)*dx(domain=mesh))
cell = mesh.ufl_cell()
fe_rho = FiniteElement("P", cell, 1)
fe_lagrange = FiniteElement("R", cell, 0)
fe = MixedElement([fe_rho, fe_lagrange])
V = FunctionSpace(mesh, fe)
Qbar = project(Constant(1/vol), V.sub(0).collapse())
v_rho, v_l = TestFunctions(V)
u = Function(V)
rho, _ = split(u)
F = inner(grad(rho)[0], v_rho) * dx + inner(rho - Qbar, v_l) * dx
solve(F == 0, u)
Seems to work. What looks a bit strange to me is that along with the “test scalar” l
comes a new unknown rho, _ = split(u)
that I just ignore. Does this confuse solvers? Seems bad for condition numbers etc?
dokken
July 31, 2020, 2:23pm
7
Add it to the system as done in the Poisson demo. Otherwise you end up with a Zero entry.
As you have a linear system, i suggest trying to solve it with an LU solver first, as it should then become clear that you Need an additional entry.
jw3126
August 2, 2020, 9:45pm
9
Thanks! The following works for me:
from fenics import *
import numpy as np
mesh = UnitIntervalMesh(4)
vol = assemble(Constant(1)*dx(domain=mesh))
cell = mesh.ufl_cell()
fe_rho = FiniteElement("P", cell, 1)
fe_lagrange = FiniteElement("R", cell, 0)
fe = MixedElement([fe_rho, fe_lagrange])
V = FunctionSpace(mesh, fe)
Qbar = project(Constant(1/vol), V.sub(0).collapse())
v_rho, v_l = TestFunctions(V)
u = TrialFunction(V)
rho, l = split(u)
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0)
facedim = mesh.topology().dim() - 1
boundary_markers = MeshFunction("size_t", mesh, facedim, 0)
LeftBoundary().mark(boundary_markers, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
F = inner(grad(rho)[0], v_rho) * dx + inner(rho - Qbar, v_l) * dx + inner(rho -l, v_rho) * ds(1)
sol = Function(V)
A = assemble(lhs(F)).array()
print(f"A = {A}")
print(f"cond(A) = {np.linalg.cond(A)}") # better way to compute cond without going to dense matrix?
solve(lhs(F) == rhs(F), sol)
rho, l = split(sol)
plot(rho)
It feels slightly weird to make up the boundary condition rho = l
though.