Dear all,
As a test, I would like to solve a simple Stokes flow problem on a rectangular domain, where the aspect ratio is small. Particularly, the length is 1/asp times the width, where asp is the aspect ratio. The velocity field is similarly scaled, so the transversal velocity u[0] is asp smaller than the longitudinal one u[1].
In order to make sure that I am doing the scaling correctly, I solved the simplest problem with aspect ratio scaling and again without any scaling, and used the Transform filter in Paraview to see if I can go from one to the other. I let my frustration get the best of me and ended up doing lots of trial and error until I got it right. Now I have a solution that I can rescale using Transform in Paraview and it looks just like the rectangular version. However, looking at the code, I don’t understand why… As far as I understand it, the velocity scaling cancels out with the divergence operator scaling, and those terms should just be div(v). However, unless I rescale the operator, it doesn’t work. This should be wrong, I should not have any asp terms in the divergence of the velocity. I don’t understand it at all. I would really appreciate if someone could take a look at the code below:
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
# Define mesh and geometry - We solve for half of the domain we need, and impose symmetry
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)
n = FacetNormal(mesh)
# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))
# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))
# Define the viscosity and bcs
mu = Constant(1.0)
u_in = Constant(-2.0)
u_c = Constant(-1.0)
L = 5.0
R = 1
asp = R/L
# Note, x[0] is r and x[1] is x, and x[1] == 0 is the bottom.
inflow = 'near(x[1], 1.0) && x[0]<=0.5'
wall = 'near(x[0], 1.0)'
centre = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(0), (0.0, u_in), inflow)
bcu_wall = DirichletBC(W.sub(0), (0.0, u_c), wall)
bcu_outflow = DirichletBC(W.sub(0), (0.0, u_c), outflow)
bcu_symmetry = DirichletBC(W.sub(0).sub(0), Constant(0.0), centre)
bcs = [bcu_inflow, bcu_wall, bcu_symmetry, bcu_outflow]
# Define stress tensor
x = SpatialCoordinate(mesh)
# Here, I am rescaling asp*x[0], so dx(0) -> 1/asp * dx(0).
def epsilon(v):
return (as_tensor([[(1/asp)*v[0].dx(0), v[0].dx(1) + (1/asp)*v[1].dx(0), 0],
[(1/asp)*v[1].dx(0) + v[0].dx(1), 2*v[1].dx(1), 0],
[0, 0, 0]])
# Define the variational form
f = Constant((0, -1)) # forcing term
# The vectors defined in Fenics are automatically dimensional. We introduce the aspect ratio here, turning the dimensional
# We defined rescaled u and v that are actually rescalings for the divergence operator. So v[0]*(1/asp).dx(0) + v[1].dx(1) ... is the same as div(vsc)
vsc = as_vector([(1/asp)*v[0], v[1]])
usc = as_vector([(1/asp)*u[0], u[1]])
colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0) # default to zero
# We match the colours to the defined sketch in the Fenics chapter
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 5)
CompiledSubDomain("near(x[1], 1.0) && x[0]<=0.5").mark(colors, 1)
CompiledSubDomain("near(x[1], 1.0) && x[0]>=0.5").mark(colors, 2)
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 3) # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 4) # outflow
# Create the measure
ds = Measure("ds", subdomain_data=colors)
# For the governing equations, we have to also multiply by the determinant of the Jacobian, which is asp. All of the
# terms have to be multiplied by asp, which, as they are equal to zero, means we can just remove asp. Note as well that
# the Jacobian for this rescaling is symmetric, so inv(J)*Transp(J) = I, so everything is also multiplied by I.
a1 = -(inner(mu*epsilon(u), epsilon(v))) * dx - p*div(vsc) * dx
a2 = (- div(usc) * q - dot(f, vsc)) * dx # why should I rescale velocity for the divergence here and above? This operation should be asp*u[0]*(1/asp).dx(0) + u[1].dx(1) so the aspect ratio cancels out! but it doesn't work without this ...
F = a1 + a2
# Solve problem
solve(F == 0, w, bcs)
# Solutions
(u, p) = w.split()
What this code produces is the right answer:
If I do what I think should be right thing, which is NOT rescaling the divergence operator (as it cancels out with the velocity rescaling), I get something clearly distorted:
If anyone could tell me why I have to apply the rescaling to div, I would really appreciate it.