### Applying source on surface does not work in case of 3D geometry, 2 physics; but works in case of 1 physic

182

views

0

Dear all

I am solving a 3D static problem with 2 physics, namely "phi" and "n", in a cube. On the top surface, a source "S" related to "n" is applied and on the bottom surface, n=0. When I run my code, the solver can not converge. However, when I test the problem with only one physic "n", it works correctly.

If you have any idea, please let me know.

Best regards

My code with "n" and "phi"

I am solving a 3D static problem with 2 physics, namely "phi" and "n", in a cube. On the top surface, a source "S" related to "n" is applied and on the bottom surface, n=0. When I run my code, the solver can not converge. However, when I test the problem with only one physic "n", it works correctly.

If you have any idea, please let me know.

Best regards

My code with "n" and "phi"

```
from __future__ import print_function
from fenics import *
import numpy as np
from mshr import *
import time
import math
# Space parameters
xLength = 100.0
yLength = 100.0
zLength = 100.0
tol = 1E-14
# Define 3D geometry
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(xLength, yLength, zLength), 30, 30, 10)
# Define function space
Q = FiniteElement('Lagrange', tetrahedron, 1) #phi space
V = FiniteElement('Lagrange', tetrahedron, 1) #n space
element = Q*V
W = FunctionSpace(mesh, element)
# Variational problem
unks = Function(W)
tests = TestFunction(W)
phi, n = split(unks)
v, w = split(tests)
# Define boundary conditions
class GammaTop(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], zLength, tol)
class GammaBottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0.0, tol)
class GammaLateral(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[0], 0.0, tol) or near(x[0], xLength, tol) or near(x[1], 0.0, tol) or near(x[1], yLength, tol))
gamma_top = GammaTop()
gamma_bottom = GammaBottom()
gamma_lateral = GammaLateral()
# Define initial value
sigma_x = 5.0
sigma_y = 5.0
X_spot = xLength/2.0
Y_spot = yLength/2.0
S_str = 'exp(-(pow((x[0] - X_spot)/sigma_x, 2.0) + pow((x[1] - Y_spot)/sigma_y, 2.0)))'
S = Expression( S_str, degree = 2,
X_spot = X_spot, Y_spot = Y_spot,
sigma_x = sigma_x, sigma_y = sigma_y, xLength = xLength,
yLength = yLength)
# Define boundary condition
phi_bcBottom = DirichletBC(W.sub(0), Constant(0.0), gamma_bottom)
n_bcBottom = DirichletBC(W.sub(1), Constant(0.0), gamma_bottom)
n_bcTop = DirichletBC(W.sub(1), S, gamma_top)
bcs = [ phi_bcBottom, n_bcBottom, n_bcTop, ]
F1 = dot(grad(phi), grad(w))*dx - n*w*dx
F2 = dot(grad(n), grad(v))*dx
F = F1 + F2
solve(F == 0, unks, bcs)
plot(n)
# Hold plot
interactive()
```

My code with "n" only

```
from __future__ import print_function
from fenics import *
import numpy as np
from mshr import *
import math
# Space parameters
xLength = 100.0
yLength = 100.0
zLength = 100.0
tol = 1E-14
# Define 3D geometry
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(xLength, yLength, zLength), 30, 30, 10)
# Define function space
V = FunctionSpace (mesh, 'Lagrange', 1) #n space
# Variational problem
n = Function(V)
v = TestFunction(V)
# Define boundary conditions
class GammaTop(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], zLength, tol)
class GammaBottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0.0, tol)
class GammaLateral(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[0], 0.0, tol) or near(x[0], xLength, tol) or near(x[1], 0.0, tol) or near(x[1], yLength, tol))
gamma_top = GammaTop()
gamma_bottom = GammaBottom()
gamma_lateral = GammaLateral()
# Define initial value
sigma_x = 5.0
sigma_y = 5.0
X_spot = xLength/2.0
Y_spot = yLength/2.0
S_str = 'exp(-(pow((x[0] - X_spot)/sigma_x, 2.0) + pow((x[1] - Y_spot)/sigma_y, 2.0)))'
S = Expression( S_str, degree = 2,
X_spot = X_spot, Y_spot = Y_spot,
sigma_x = sigma_x, sigma_y = sigma_y, xLength = xLength,
yLength = yLength)
# Define boundary condition
n_bcBottom = DirichletBC(V, Constant(0), gamma_bottom)
n_bcLateral = DirichletBC(V, Constant(1.0), gamma_lateral)
n_bcTop = DirichletBC(V, S, gamma_top) #Constant(1.0)
bcs = [ n_bcBottom, n_bcTop, ]
F1 = dot(grad(n), grad(v))*dx
F = F1
solve(F == 0, n, bcs)
plot(n)
# Hold plot
interactive()
```

result

Community: FEniCS Project

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