Dear all,
i tried to implement a method to solve Boussinesq equation for square cavity, but i I can not understand the type of error.
import os
import shutil
import numpy as np
from dolfin import *
from RBniCS import *
from sys import argv
from distutils.version import LooseVersion as version
import scipy.io as sio
Mesh
mesh = UnitSquareMesh(20,20)
Function space
if version(dolfin_version()) < version("1.7.0"):
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
E = FunctionSpace(mesh, "Lagrange", 1)
W = MixedFunctionSpace([V, Q, E])
else:
V_element = VectorElement('Lagrange', mesh.ufl_cell(), 2)
Q_element = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
E_element = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
W_element = MixedElement(V_element, Q_element, E_element)
W = FunctionSpace(mesh, W_element)
Boundary conditions
noslip = DirichletBC(W.sub(0), (0, 0), "x[0] < DOLFIN_EPS || x[0] > 1.0 - DOLFIN_EPS || x[1] < DOLFIN_EPS || x[1] > 1 - DOLFIN_EPS")
pref = DirichletBC(W.sub(1), 0, "x[0] < DOLFIN_EPS && x[1] < DOLFIN_EPS", "pointwise")
t_hot = Constant(1.0)
t_cold = Constant(0.0)
Th = DirichletBC(W.sub(2), t_hot, "x[0] < DOLFIN_EPS")
Tc = DirichletBC(W.sub(2), t_cold, "x[0] > 1.0 - DOLFIN_EPS")
Collect BCs
bc_list = [noslip, pref, Th, Tc]
Define unknown and test function(s)
(v, q, z) = TestFunctions(W) # v: velocity, q: pressure, z: temperature
Define trial functions
w = TrialFunction(W)
(u, p, e) = split(w)
current unknown time step
w1 = Function(W)
(u1, p1, e1) = split(w1)
previous known time step
w0 = Function(W)
(u0, p0, e0) = split(w0)
Fixed point iterations
w_prev = Function(W)
(u_prev, p_prev, e_prev) = split(w_prev)
Dimentionless numbers
Re = Constant(30)
Ra = Constant(1000)
Pr = Constant(0.71)
dt = Constant(0.1)
T = Constant(0.6)
uy = Expression(('0.0','1'))
eps = 1.0
tol = 1.0e-4
niter = 0
maxit = 1000
counter = 0
cost1 = Ra/(PrReRe)
cost2 = sqrt(Ra*Pr)
f = cost1uye0
Variational forms
Fa = ( Constant(1./dt)inner(u, v)dx
+ Constant(1./dt)ezdx
+ Constant(1./Re)inner(grad(u), grad(v))dx
+ Constant(cost2)inner(grad(e), grad(z))dx
+ inner(grad(u)u_prev, v)dx
- div(v)pdx
+ div(u)qdx
+ inner(u, grad(e))z*dx
)
FL = ( inner(f, v)dx
+ Constant(1./dt)inner(e0, z)*dx
)
F = Fa + FL
a = lhs(F)
L = rhs(F)
print type(lhs)
print type(rhs)
Define inner product matrices for stopping criterion
x_u = inner(u, v)*dx
X_u = assemble(x_u)
XXX Fixed-point loop
norm_increment_u = 2tol
norm_increment_p = 2tol
w = Function(W)
delta_upt = Function(W) # Function in the mixed space to store the increment
it = 0
for it in range(maxit):
# Solve
solve(a == L, w, bc_list)
# Extract solutions:
(u, p, e) = w.split()
# Compute the (relative) norm of the increment
delta_upt.vector().zero()
delta_upt.vector().add_local(w.vector().array())
delta_upt.vector().add_local(- w_prev.vector().array())
delta_upt.vector().apply("")
print type(delta_upt)
#delta_upt.vector()[:] = upt.vector()[:] - upt_prev.vector()[:]
norm_increment_u = delta_upt.vector().inner(X_u*delta_upt.vector())
print type(norm_increment_u)
print norm_increment_u
# Update the previous value
assign(w_prev, w)
# Print to screen
print "Iteration", it, \
"|| delta_u||_V =", norm_increment_u
# Check stopping criterion
if (norm_increment_u < tol and norm_increment_u < tol):
print "Fixed point loop has converged after", it, "iterations."
break
Best reguards.
Francesco
