Hello everyone. I’ve been trying to simulate the flow of two fluids with different dynamic viscosity. I’ve found how to create subdomains, so that each domain has a different viscosity, however, when I attempt to assemble the matrix for the navier stokes problem, it fails to do so, giving this error:
216 # Call C++ assemble
217 assembler.assemble(tensor, dolfin_form)
218
219 # Convert to float for scalars
RuntimeError:
It appears to fail when I try to assemble matrix A
The code I have looks like this:
from dolfin import *
#Creación de un mallado. Se va a refinar en los bordes inferiores y superiores
n = 100
mesh = UnitSquareMesh(n,n)
#La malla puede refinarse si es necesario
#Define Taylor-Hood elements (For now, they work best for Navier-Stokes problems)
V = VectorFunctionSpace(mesh, 'DG', 2) #Function space that solves for velocity
Q = FunctionSpace(mesh, 'DG', 1) #Function space that solves for pressure
#Define Two sets of Test/Trial functions, as we are working with two unknowns.
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
class Omega0(SubDomain):
def inside(self, x, on_boundary):
return True if x[0] <= 0.5 + tol else False
class Omega1(SubDomain):
def inside(self, x, on_boundary):
return True if x[0] >= 0.5 - tol else False
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomain0 = Omega0().mark(subdomains, 0)
subdomain1 = Omega1().mark(subdomains, 1)
class K(UserExpression):
def __init__(self, subdomains, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.subdomains = subdomains
self.k_0 = k_0
self.k_1 = k_1
mu1 = Constant(0.1)
mu2 = Constant(1)
mu = K(subdomains=subdomains, k_0 = mu1, k_1 = mu2, degree=0)
def up_boundary(x, on_boundary):
tol = 1E-8
return(on_boundary and abs(x[1]-1)< tol)
def down_boundary(x, on_boundary):
tol = 1E-8
return(on_boundary and abs(x[1])< tol)
def left_boundary(x, on_boundary):
tol = 1E-8
return(on_boundary and abs(x[0])< tol)
def right_boundary(x, on_boundary):
tol = 1E-8
return(on_boundary and abs(x[0]-1)< tol)
def pressureCond(x, on_boundary):
tol = 1E-8
return(on_boundary)
DBup = DirichletBC(V, Constant((-1.0,0)) , up_boundary) #Note that constants must be two dimentional
DBdown = DirichletBC(V, Constant((0,0)) , down_boundary)
DBleft = DirichletBC(V, Constant((0,0)) , left_boundary)
DBright = DirichletBC(V, Constant((0,0)) , right_boundary)
DBcond = DirichletBC(Q, Constant(0) , pressureCond)
bcs = [DBup, DBdown, DBleft, DBright] #Add all Dirichlet BC into a list
bcs2 = [DBcond]
# Define functions for solutions at previous and current time steps
u_n = Function(V) #Es importante definirlas puesto que voy a iterar en el tiempo
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
#Formulation of the variational form of NavierStokes
#Parameter Creation
T = 1 # final time
num_steps = 500 # number of time steps
dt = T / num_steps # time step size
#mu = 0.01 # dynamic viscosity
rho = 1 # density
ucond = 1
reyndods = (rho * ucond)/mu
# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh) #Esto me crea un vector normal a la superficie planteada
f = Constant((0,0)) #Fuerza constante, la idea es ir cambiando este parametro a ver que ocurre
k = Constant(dt)
mu = UserExpression(mu)
rho = Constant(rho) #Integers must be changed into constant type
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))
#Variational form due to Chorin's Projection.
#Velocity without pressure calculation
F1 = rho*dot((u - u_n) / k, v)*dx + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx + inner(sigma(U, p_n), epsilon(v))*dx + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
#Pressure Correction
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx
#Velocity correction due to pressure
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
[bc.apply(A1) for bc in bcs] #Apply boundary conditions to A
# Create XDMF files for visualization output
xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')
# Save mesh to file (for use in reaction_system.py)
File('navier_stokes_cylinder/cylinder.xml.gz') << mesh
# Time-stepping
import matplotlib.pyplot as plt
t = 0
i = 0
for n in range(num_steps):
# Update current time
t += dt
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcs]
solve(A1, u_.vector(), b1 ,"gmres", "ilu",) #Preconditioned settings for linear solving methods that accelerate process
# Step 2: Pressure correction step
b2 = assemble(L2)
#[bc.apply(b2) for bc in bcs2]
solve(A2, p_.vector(), b2, "gmres", "amg",)
# Step 3: Velocity correction step
b3 = assemble(L3)
#[bc.apply(b3) for bc in bcs2]
solve(A3, u_.vector(), b3, "gmres", "ilu",)
#If we wish to plot magnitude of velocity, we must plot norm (Only works on 2D)
x1,y1=u_.split(deepcopy=True)
norm = (x1*x1 + y1*y1)**(0.5)
name = "output" + "a"*i + ".jpg"
name = "check" + str(i) + ".jpg"
a = plot(u_ , title = "Trial simulation")
plt.colorbar(a)
#plt.draw()
plt.savefig(name) #save as jpg
#plt.pause(0.02)
plt.clf()
# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
i = i+1