Hi @dokken ! thanks for your answer ! The NaN was caused by || used instead of and in boundary conditions definition. Then I didn’t get the fact that ds element uses a marker meshfunction of size n-1 instead of n, thus the code wasn’t finding the facets needed to apply boundary conditions. Took me a while to figure that out ! All sorted out now.
For information, this code works:
'''
this code simulates 2D heat transfer between flowing water and solid stainless steel (ss).
'''
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
T = 10.0 # final time
num_steps = 10 # number of time steps
dt = T / num_steps # time step size
#tensors needed for flow equation:
def epsilon(u):
return sym(nabla_grad(u))
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))
#domain and subdomains definition as in demo_block-assembly-2D2D-nonlinear.py : https://bitbucket.org/fenics-project/dolfin/src/39253449ee544fe04eb8bf387eab115e094d325e/python/demo/undocumented/block-assembly-2D2D/demo_block-assembly-2D2D-nonlinear.py?at=cecile%2Fmixed-dimensional
domain = Rectangle(Point(0, 0), Point(16,8))
domain_water = Rectangle(Point(0, 0), Point(16,4))
domain_ss = domain-domain_water
domain.set_subdomain(1, domain_water)
domain.set_subdomain(2, domain_ss)
#mesh definition
mesh = generate_mesh(domain,16)
# subdomain markers, facet markers and integration measures
markers = MeshFunction('size_t', mesh, mesh.topology().dim(), mesh.domains())
plot(markers)
facetmarkers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
dx = Measure("dx", domain=mesh, subdomain_data=markers)
ds = Measure("ds", domain=mesh, subdomain_data=facetmarkers)
print(markers.array(), min(markers.array()), max(markers.array()))
print(facetmarkers.array(), min(facetmarkers.array()), max(facetmarkers.array()))
#physical properties defined as in fenics tutorial ft11 updated : https://github.com/hplgit/fenics-tutorial/issues/65 :
#water viscosity
mu=1 #Pa.s
#Density :
class Density(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = 1. #water, no units
else:
values[0] = 7.9 #stainless steel, no units
def value_shape(self):
return ()
rho = Density(markers, degree=1)
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
#DEFINITION OF BOUNDARY CONDITIONS USING SECOND MARKER FUNCTION---------------------------------
#source:https://fenicsproject.discourse.group/t/no-slip-boundary-condition-not-being-applied/3924/2
#inflow = 'near(x[0], 0, 1E-10) and x[1]<=4+1E-10'
#outflow = 'near(x[0], 0, 1E-10) and x[1]<=4+1E-10'
#walls = 'near(x[1], 0, 1E-10) || near(x[1], 4, 1E-10)'
# Sub domain for walls (mark whole boundary, inflow and outflow will overwrite)
class Walls(SubDomain):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def inside(self, x, on_boundary):
return on_boundary or near(x[1], 4, 1E-10)
# Sub domain for inflow (right)
class Inflow(SubDomain):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def inside(self, x, on_boundary):
return near(x[0], 0, 1E-10) and x[1]<=4+1E-14 and on_boundary
# Sub domain for outflow (left)
class Outflow(SubDomain):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def inside(self, x, on_boundary):
return x[0] >16-1E-14 and x[1]<=4+1E-14 and on_boundary
# Mark all facets as facetmarkers sub domain 3
facetmarkers.set_all(3)
# Mark no-slip facets as facetmarkers sub domain 0,
walls = Walls(markers)
walls.mark(facetmarkers, 0)
# Mark inflow as facetmarkers sub domain 1,
inflow = Inflow(markers)
inflow.mark(facetmarkers, 1)
# Mark outflow as facetmarkers sub domain 2
outflow = Outflow(markers)
outflow.mark(facetmarkers, 2)
print(facetmarkers.array(), min(facetmarkers.array()), max(facetmarkers.array()))
#inflow = 'near(x[0], 0, 1E-10) and x[1]<=4+1E-10'
#outflow = 'near(x[0], 0, 1E-10) and x[1]<=4+1E-10'
#walls = 'near(x[1], 0, 1E-10) and near(x[1], 4, 1E-10)'
bcu_noslip = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow = DirichletBC(Q, Constant(8), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip]
bcp = [bcp_inflow, bcp_outflow]
#FUNCTIONS DEFINITION-----------------------------------------------------------------------
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
f = Constant((0, 0))
k = Constant(dt)
mu = Constant(mu)
#IPCS IMPLEMENTATION---------------------------------------------------------------------------
#Source :fenics tutorial demo ft07: https://github.com/gdmcbain/fenics-tuto-in-skfem/commit/3d4803e2fd61a7672cd7d1a4fb4f4ccfedde6084
# Define variational problem for step 1
F1 =rho*dot((u - u_n) / k, v)*dx(1) + \
rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx(1) \
+ inner(sigma(U, p_n), epsilon(v))*dx(1) \
+ dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
- dot(f, v)*dx(1)
a1 = lhs(F1)
L1 = rhs(F1)
a2 = dot(nabla_grad(p), nabla_grad(q))*dx(1)
L2 =dot(nabla_grad(p_n), nabla_grad(q))*dx(1) - (1/k)*div(u_)*q*dx(1)
a3 = dot(u, v)*dx(1)
L3 = dot(u_, v)*dx (1)- k*dot(nabla_grad(p_ - p_n), v)*dx(1)
A1 = assemble(a1, keep_diagonal=True)
A1.ident_zeros()
A2 = assemble(a2, keep_diagonal=True)
A2.ident_zeros()
A3 = assemble(a3, keep_diagonal=True)
A3.ident_zeros()
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
xdmffile_u = XDMFFile('example/navstok.xdmf')
t = 0
for n in range(num_steps):
t += dt
b1 = assemble(L1) # Step 1: Tentative velocity step
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1)
print(min(u_.vector().get_local()), max(u_.vector().get_local()))
b2 = assemble(L2) # Step 2: Pressure correction step
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2)
print(min(u_.vector().get_local()), max(u_.vector().get_local()))
b3 = assemble(L3) # Step 3: Velocity correction step
solve(A3, u_.vector(), b3)
print(min(u_.vector().get_local()), max(u_.vector().get_local()))
u_n.assign(u_)
p_n.assign(p_)
xdmffile_u.write(u_, t)
plot(u_)
plt.show()
link to the Github
