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 »