Applying Boundary Conditions to 1D Problem with Nonlinear Solver


139
views
2
4 months ago by

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
I'm having a hard time finding documentation on this solver's way of handling boundary conditions ... perhaps someone can suggest where to find that?
written 4 months ago by ricopicone  

1 Answer


5
4 months ago by
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.

Similar posts:
Search »