### Applying Boundary Conditions to 1D Problem with Nonlinear Solver

139

views

2

I previously switched to a different solver because I was having trouble getting with my solutions not converging,

https://www.allanswered.com/post/lkbwa/newton-solver-not-converging-with-nonlinear-system/

I am now trying to implement nonzero Neumann boundary conditions on my system, but they appear to be working better as zero Neumann boundary conditions. What am I doing wrong in implementing these boundary conditions?

In the future I may also want to apply Dirichlet boundary conditions, how would I go about doing this?

A minimal working example is below.

```
from fenics import *
import matplotlib.pyplot as plt
T = 0.1
num_steps = 50
L = 15
n = 50
dt = T / num_steps
mesh = IntervalMesh(n, -L/2, L/2)
V = VectorFunctionSpace(mesh, 'CG', 1, dim = 3)
v1, v2, v3 = TestFunctions(V)
rho = Function(V)
rho_n = Function(V)
rho1, rho2, rho3 = split(rho)
rho1_n, rho2_n, rho3_n = split(rho_n)
rho_ic = Expression(['0', '0.1', '0.2'], degree = 1)
rho20 = rho_ic[1]
rho30 = rho_ic[2]
rho.assign(rho_ic)
rho_n.assign(rho_ic)
_dt = Constant(dt)
dr = dx
F = ((rho1 - rho1_n) * v1 / _dt
+ (rho2 - rho2_n) * v2 / _dt
+ (rho3 - rho3_n) * v3 / _dt) * dr
class Edge(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[0], -L/2) or near(x[0], L/2))
Edge().mark(FacetFunction("size_t", mesh, 0), 1)
G1 = Constant(1)
G2 = Constant(-1)
G3 = Constant(2)
F += (v1 * G1 + v2 * G2 + v3 * G3) * ds(1)
J = derivative(F, rho)
class NLP(NonlinearProblem):
def F(self, b, x):
assemble(F, tensor=b)
def J(self, A, x):
assemble(J, tensor=A)
solver = PETScSNESSolver()
PETScOptions.set("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
PETScOptions.set("pc_factor_mat_solver_type", "mumps")
PETScOptions.set("snes_atol", 1e-8)
solver.set_from_options()
for i in range(1, num_steps + 1):
solver.solve(NLP(), rho.vector())
rho_n.assign(rho)
plt.figure()
plt.subplot(3,1,1)
plot(rho1)
plt.subplot(3,1,2)
plot(rho2)
plt.subplot(3,1,3)
plot(rho3)
plt.show()
```

Community: FEniCS Project

written
4 months ago by
ricopicone

### 1 Answer

5

I never used the PETSc solver but apparently it also uses the NonlinearProblem interface. So you can apply the Dirichlet boundary conditions in the definition of your NonlinearProblem:

```
class NLP(NonlinearProblem):
def F(self, b, x):
assemble(F, tensor=b)
[bc.apply(b,x) for bc in dbcList]
def J(self, A, x):
assemble(J, tensor=A)
[bc.apply(A) for bc in dbcList]
```

where `dbcList`

is a python list of DirichletBC objects.

Thank you. This looks like exactly what I need to solve half of my problem. However, sometimes I would like to apply Neumann boundary conditions. How would I go about doing this?

written
4 months ago by
Cameron Devine

Neumann BCs are part of the weak formulation. These will be included in your residual

`F = ...`

written
4 months ago by
Nate

This was what I was attempting to do at the line starting with F += ...

written
4 months ago by
Cameron Devine

Please login to add an answer/comment or follow this question.