Neumann BC Problem


114
views
0
3 months ago by

I am applying Neumann Boundary Conditions to a 1d problem. However the result does not appear to be abiding by these boundary conditions. What needs to be changed for these boundary conditions to work?

A minimal working example is below.

from fenics import *
import math
import numpy as np

# Defining an invese hyperbolic tangent because fenics doesn't have one
def atanh(x): # Only good for abs(x) < 1!
	return ln((1 + x) / (1 - x)) / 2

T = 0.1
num_steps = 50
L = 15
n = 50

temp = 10
kB = 1.3806488e-23
Grad = 44e3
B0 = 0.0893

γe = 1.760859708e11
γp = 2.675222005e8
γ2 = γp
γ3 = γe
γb = γ3/γ2

NA = 6.02214129e23
MwPS = 12*0 + 1*8
MwDPPH = 394.32
dDPPH = 1.4
nAMDPPH = 1
concDPPH = .04
den3 = 1/MwDPPH*dDPPH*1e6*NA*nAMDPPH
δ3 = concDPPH*den3
dPS = 1.047
nAMPS = 2
concPS = 1 - concDPPH
den2 = 1/MwPS*dPS*1e6*NA*nAMPS
δ2 = concPS*den2
δb = δ3/δ2

hb = 1.054571726e-34
μe = -hb*γe/2 
μp = hb*γp/2

π = math.pi
μ = 4*π*1e-7
Bd2 = μ/(4 * π)*hb*γ2*δ2
Bd3 = μ/(4 * π)*hb*γ3*δ3
Bd = Bd2 + Bd3

Bb = 1
cb = Bb*(1+δb)/(1+γb*δb)

Γp = μ/(4 * π)*hb*γ2**2*δ2**(1/3)
Γe = μ/(4 * π)*hb*γ3**2*δ3**(1/3)
Γ2 = Γp
Γ3 = Γe
Γb = Γ3/Γ2

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)

Bt = Expression('b0 + bd * x[0]', b0 = B0, bd = Bd, degree = 1)

ic = (
		'tanh(((1 + gamma * Delta) / (gamma * (1 + Delta))) * (mu3 * Bd / (kB * Temp)))',
		'tanh(mu2 * Bt / (kB * Temp))',
		'tanh(mu3 * Bt / (kB * Temp))')

rho_ic = Expression(ic,
	gamma = γb,
	Delta = δb,
	mu2 = μp,
	mu3 = μe,
	Bd = Bd,
	Temp = temp,
	kB = kB,
	grad = Grad,
	Bt = Bt,
	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

F += -((cb**2 / (1 + δb)) * ((1 - rho2**2) + Γb * δb * γb**2 * (1 - rho3**2)) * v1 * atanh(rho1)
  - (cb / (1 + δb) * (grad(rho2)[0] - Γb * δb * γb * grad(rho3)[0]) * v1)) * dr \
  + ((-cb * ((1 - rho2**2) / (1 - rho1**2)) * grad(rho1)[0]
  + 2 * cb * rho2 * atanh(rho1) * grad(rho2)[0]) * v2) * dr \
  + ((-cb * ((1 - rho3**2) / (1 - rho1**2)) * grad(rho1)[0]
  + 2 * cb * rho3 * atanh(rho1) * grad(rho3)[0]) * v3) * dr

F += (1 + Γb) * grad(rho1)[0] * grad(v1)[0] * dr \
  + grad(rho2)[0] * grad(v2)[0] * dr \
  + Γb * grad(rho3)[0] * grad(v3)[0] * dr

tol = 1e-5

class LeftEdge(SubDomain):
	def inside(self, x, on_boundary):
		return on_boundary and near(x[0], -L/2, tol)

class RightEdge(SubDomain):
	def inside(self, x, on_boundary):
		return on_boundary and near(x[0], L/2, tol)

boundaryLeft = LeftEdge()
boundaryRight = RightEdge()

boundaries = FacetFunction("size_t", mesh)
boundaryLeft.mark(boundaries,1)
boundaryRight.mark(boundaries,2)
ds = Measure('ds')[boundaries]

G1 = Constant(0)
G2 = Constant(5.98e-8)
G3 = Constant(-3.93e-3)

F += (v1 * G1 + v2 * G2 + v3 * G3) * ds(1)
F -= (v1 * G1 + v2 * G2 + v3 * G3) * ds(2)

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)

x = np.linspace(-L/2, L/2, n + 1)
dx = np.diff(x[:2])[0]
rho = rho.vector().get_local().reshape([x.shape[0], len(rho.split())])
rho1 = rho[::-1,0]
rho2 = rho[::-1,1]
rho3 = rho[::-1,2]

print('rho1 applied Neumann BC', np.diff(rho1, dx)[[0,-1]])
print('rho2 applied Neumann BC', np.diff(rho2, dx)[[0,-1]])
print('rho3 applied Neumann BC', np.diff(rho3, dx)[[0,-1]])

This code outputs,

rho1 applied Neumann BC [ -4.16301723e-15  -3.05960152e-15]
rho2 applied Neumann BC [  1.79462500e-08   1.79462498e-08]
rho3 applied Neumann BC [ -3.29147258e-08  -3.29147258e-08]
Community: FEniCS Project
2
The code is a little too large for someone to just dive in. You could increase your chances of getting a reply by creating a _minimal_ example that highlights the behavior. (In my experience, the process of creating a minimal example often leads to enlightenment already.)
written 3 months ago by Nico Schlömer  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »