Multidimensional coupled equations

General
1

I’m playing with the mixed-dimensional branch by @cdaversin, specifically exploring the example script in /demo. After fixing some of the outdated syntax, I attempted to adapt the problem so that instead of u = c along Gamma, I tried to set u = x[1] along Gamma, so that u would linearly increase towards the top.
This gave the following error

*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 10000 iterations (PETSc reason DIVERGED_ITS, residual norm ||r|| = 7.869592e-05).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.

I can’t understand why this problem would be ill-posed or be difficult to solve using the linear solver. Is something poorly defined? Am I giving UFL conflicting constraints? Is the problem with the initial guess for the nonlinear solver?

Any assistance would be greatly appreciated. :blush:

Below is my code:

from dolfin import *

### Step 1 : Build the meshes ###

# Parent mesh (Omega) : very basic unit square
n = 30
mesh = UnitSquareMesh(n, n)

# Submesh (Gamma, for the Lagrange multiplier) is the line x=0.5
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for f in facets(mesh):
    marker[f] = 0.5 - DOLFIN_EPS < f.midpoint().x() < 0.5 + DOLFIN_EPS
gamma_mesh = MeshView.create(marker, 1)

### Step 2 : Define the function spaces ###

V = FunctionSpace(mesh, "CG", 1)
LM = FunctionSpace(gamma_mesh, "DG", 0)
W = MixedFunctionSpace(V,LM)
u,l = TrialFunctions(W)
v,e = TestFunctions(W)

### Step 3 : The variational formulation ###

# Homogeneous Dirichlet on left and right (x=0 and x=1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
bc = DirichletBC(V, Constant(0.0), boundary)

# Define explicitely the integration domains we will use in the bilinear and linear forms
dV = Measure("dx", domain=W.sub_space(0).mesh()) # To integrate on the 2D domain
dL = Measure("dx", domain=W.sub_space(1).mesh()) # To integrate on the 1D submesh (gamma)
# Function used in the RHS

fu = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
##fu = Expression("10*exp(-50*(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)))", degree=2)
fl = Expression('x[1]', degree = 1) # u = c on Gamma
# Variational formulation
a = inner(grad(u),grad(v))*dV + v*l*dL + u*e*dL
L = fu*v*dV + fl*e*dL

### Step 4 : Solve the problem ###
sol = Function(W)
solve(a == L, sol, bc)
# Note : The solver can be configured through the solver_parameters={}, as usual....
# solve(a == L, sol, bc, solver_parameters={"krylov_solver":{"absolute_tolerance":1e-6, "monitor_convergence":True}})
2

Hi,

You are using the default iterative solver (gmres), you could try to use another iterative solver (for ex : minres), or use a direct solver since your problem is small :

# Direct solver
solve(a == L, sol, bc, solver_parameters={"linear_solver":"direct"})
# Minres solver
solve(a == L, sol, bc, solver_parameters={"linear_solver":"minres"})

I have tested those two in the latest Docker container of FEniCS, and it seems to work. You could try other linear solvers as well (https://fenicsproject.org/pub/tutorial/html/._ftut1017.html).

2 Likes
3

Great! Thanks this was quite helpful, especially in reminding me about the tutorial page on the linear solvers!