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

182
views
0
8 months ago by
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"
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.

Similar posts:
Search »