This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Free Convection

–1 vote

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 = 2
tol

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

asked Oct 24, 2016 by Francesco_nuke_91 FEniCS Novice (280 points)

Please use proper code formatting and describe the kind of error you are experiencing.

...