### How to include plane averages in variational form

378
views
0
9 months ago by
I'm having some trouble implementing a sub grid scale model to the Boussinesq equations.
$\frac{D\overline{u}}{Dt}=-grad\left(P\right)-\left(f\times\overline{u}+b\right)+tau_{ij}$DuDt =grad(P)(ƒ ×u+b)+tauij
where D/Dt is the material derivative, P is the pressure, f is the Coriolis force and b is a buoyancy only in the z-direction.
The tensor includes the average of the strain rate tensor in the homogeneous directions (x,y), which makes it a function of z. The brackets imply a homogeneous average.
$tau_{ij}=-2v_t\gamma S_{ij}-2v_T$tauij=2vtγSij2vT<Sij>
This vT also depends on other homogeneous averages, namely

$v_T=\left(k_1-<\gamma v_t>-k_2\cdot\sqrt{\left(uw\right)^2+\left(vw\right)^2}\right)$vT=(k1<γvt>k2·(uw)2+(vw)2) *sqrt( 2 <Sij> <Sij> )
where k1 and k2 are just constants and I've defined how to calculate gamma.
My problem is I'm not sure how to incorporate this homogeneous average from the bracket operators into the variational form if it will be constantly changing.

from __future__ import print_function
from fenics import *
from ufl import *
#parameters["plotting_backend"]="matplotlib"
#from mshr import *
import numpy as np
# # TIME SPECS
T, dt = 0.03, 0.01
k = dt # For variational form
# domain dimensions
Lx, Ly, Lz = 320, 320, 96
rsl = 20

# For eddy viscosity delta
d_x = 1.5 * (Lx/rsl)
d_y = 1.5 * (Ly/rsl)
d_z = Lz/rsl
delta = (d_x*d_y*d_z) ** (1/3)

# # PARAMETERS
Coriolis = 0.0000729 # Coriolis magnitude
g= 9.81 #for buoyancy implementation
#Coriolis = n*0.0000729 # Coriolis magnitude
rho0, rhoair = 1000.0, 1.0 # reference H2O density [kg/m^3] , air
#g= -9.81 #for buoyancy implementation
B_T = 2.0 * 10.0 ** -4.0 # H2O thermal expansion coefficient - buoyancy
theta0 = 290.16 # reference temperature [K]
U10 = 5.75 # U_10 - 10 meter wind speed for BC [m/s]
Cdw = 8.5 * 10 ** -4 # Cd - coefficient of drag on water
u_z0 = sqrt ( (rhoair/rho0) * Cdw * U10 ** 2 ) # Surface wind stress boundary condition
Ck, Ceps = 0.1, 0.93 # SGS constant (Moeng and Wyngaard, 1988)
Cs = 0.18 # Smagorinsky constant (Sullivan, 1994)
# Monin-Obukhov stabiity function, from Sullivan, 1994
etta = - float(d_z) / float(Lz) # normalize
phim_z1 = ( 1.0 - 15.0*etta)**(-0.25) # MO function at first grid point
vonk = 0.35 # von Karman constant
MOcoeff= float(vonk * d_z)/float(u_z0 * phim_z1) # coefficient for eddy viscosity

# "Forcing" functions for variational forms
f_1 = ( 0 , 0 , 0 )  # Velocity
f_3 = 0 # Temperature forcing function
f_5 = 0 # Energy forcing function

# Domain
mesh = BoxMesh(Point(0,0,0), Point(Lx,Ly,Lz), rsl,rsl,rsl)

# Initial condition for temperature
class InitialTemperature (Expression):
def eval(self, value, x):
if (x[2] >= 66):
value[0] = 298.15
else:
value[0] = 298.15 - 0.01 * x[2]

# Periodic Boundary conditions
#   Subdomain is for horizontal velocities
class PeriodicDomain(SubDomain):
def inside(self, x, on_boundary):
return bool(( near(x[0], 0) or near(x[1], 0) or near(x[2], 0)) and
(not ((near(x[0], Lx) and near(x[2], 0)) or
(near(x[0], 0) and near(x[2], Lz)) or
(near(x[1], Ly) and near(x[2], 0))or
(near(x[1], 0) and near(x[2], Lz)))) and on_boundary)
#   Mapping surfaces and edge for PBCs
def map(self, x , y ):
if near(x[0],Lx) and near(x[1],Ly):
y[0] = x[0]-Lx
y[1] = x[1]-Ly
y[2] = x[2]
elif near(x[0],Lx):
y[0] = x[0]-Lx
y[1] = x[1]
y[2] = x[2]
elif near(x[1],Ly):
y[0] = x[0]
y[1] = x[1]-Ly
y[2] = x[2]
else:
y[0]= -1000
y[1]= -1000
y[2]= -1000

pbc = PeriodicDomain() # to use in function space definition

# Function Space for:
V = VectorFunctionSpace(mesh,'CG', 1, dim=3, constrained_domain=pbc)   # velocity
Q = VectorFunctionSpace(mesh,"CG",1,dim=1, constrained_domain=pbc)         # pressure
R = VectorFunctionSpace(mesh,'CG',1,dim=1, constrained_domain=pbc) # temperature
E = VectorFunctionSpace(mesh, "CG",1,dim=1, constrained_domain=pbc)       # energy

# Velocity
u, v = TrialFunction(V), TestFunction(V)
# Pressure
p, q = TrialFunction(Q), TestFunction(Q)
# Temperature
theta, r = TrialFunction(R), TestFunction(R)
# Energy
e = TrialFunction(E)
nrg = TestFunction(E)

# Define functions for prev and current solutions
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
theta_n = Function(R)
theta_ = Function(R)
e_n = Function(E)
# Surfaces
top = 'near(x[2],96)'
bottom = 'near(x[2],0)'
# Boundary conditions
bc_top = DirichletBC(V.sub(2) , 0 , top) # No vertical velocity @ top
bc_bottom = DirichletBC(V.sub(2) , 0, bottom) # NVV @ bottom
bc_wind = DirichletBC(V.sub(0) , u_z0, top) # Wind sirface forcing
bcu = [bc_top , bc_bottom , bc_wind] # gather up BCs for function space constrained_domain

# Expression Definition
# Define Coriolis and buoyancy term for CL equations
def CL(u,theta):
return as_vector([Coriolis*u[1]] +  [-Coriolis*u[0]] + [-g*(1.0+B_T*theta0-B_T*theta)] )

# - dot(CL(u_n,theta_n),v)*dx
# TODO check variational forms
# # VARIATIONAL FORM DEFINITION
F1 = dot( (u - u_n) / k , v)*dx + dot(dot(u_n, nabla_grad(u_n)),v)*dx  + inner(tij(u_n),sigma(v))*dx
a1 = lhs(F1)
L1 = rhs(F1)

a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

a4 = theta*r*dx

a5 = e*nrg*dx
L5 = e_n*nrg*dx - k*dot(u_n,nabla_grad(e_n))*nrg*dx + P(u_n)*nrg*dx - eps(e_n, theta_n)*nrg*dx + D(e_n,u_n)*nrg*dx

A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
A4 = assemble(a4)
A5 = assemble(a5)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"
# Initialize temperature field
T0= InitialTemperature(degree=1)
theta_n= interpolate(T0,R)
u0 = Constant((0,0,0))
e0 = Constant(0.0)
u_n = interpolate(u0,V)
e_n = interpolate(e0,E)
progress = Progress('Time-stepping')
set_log_level(PROGRESS)
# Time stepping
t=0
# Vector to move thru z values for slices
for i in range(0,rsl):
z[i]=d_z*i

while t < T + DOLFIN_EPS:
## SULLIVAN SGS MODEL ##
# Create function to obtain plane averages
for x in range(0,rsl):
# Generate the strain rate tensor from previous calculations
# Define strain rate tensor
SijTensor = project(Sij, TensorFunctionSpace(mesh,'DG',0))
# Define velocity fluctuation expressions
UW = project(inner(u_n[0],u_n[2]), FunctionSpace(mesh,'CG',0))
VW = project(inner(u_n[1],u_n[2]), FunctionSpace(mesh,'CG',0))
# Use z-component
hSij = project(Expression("x[2]"),SijTensor)
hUW = project(Expression("x[2]"),UW)
hVW = project(Expression("x[2]"),VW)
# unit square mesh of top 2 in 3d
bmesh = BoundaryMesh(mesh, "exterior")
# choose side with normal in z-direction
z00 = AutoSubDomain(lambda x, on_boundary: x[2] < DOLFIN_EPS)
ff = CellFunction("size_t",bmesh,0)
z00.mark(ff,1)
smesh = SubMesh(bmesh,ff,1)
#move slice to z-location
xyz = smesh.coordinates()
xyz[:,2]=z[x]
# interpolate from h0 to the slice
ISij = FunctionSpace(smesh, 'CG', 1)
Iuw = FunctionSpace(smesh,'CG', 1)
Ivw = FunctionSpace(smesh, 'CG',1)
lp = LagrangeInterpolator()
bracSij = Function(ISij)
bracUW = Function(Iuw)
bracVW = Function(Ivw)
lp.interpolate(bracSij,hSij)
lp.interpolate(bracUW,hUW)
lp.interpolate(bracVW,hVW)
Sij_brac[x] = assemble(bracSij*dx)
uw[x] = assemble(bracUW*dx)
vw[x] = assemble(bracVW*dx)
# isotropy factor
Sprime = sqrt(2*inner((Sij - Sij_brac[x]),(Sij - Sij_brac[x])))
Sb = sqrt(2*inner(Sij_brac[x],Sij_brac[x]))
gamma = Sprime / (Sprime+Sb)
# fluctuating eddy viscosity
nut = (Cs*delta)**2)*sqrt(2*inner(Sij,Sij))
ntg = inner(nut,gamma)
NTG = project(ntg, FunctionSpace(mesh,'CG',0))
hNTG = project(Expression("x[2]"),NTG)
Ing = FunctionSpace(smesh,'CG',1)
bracNTG = Function(Ing)
lp.interpolate(bracNTG,hNTG)
nu_tgamma[x]=assemble(bracNTG*dx)
nu_Tstar = (u_z0 * vonk * d_z) / phim_z1 - nu_tgamma[x] - MOcoeff*sqrt(uw[x]**2 + vw[x]**2)
if (x == 0):
nu_T[x] = nu_Tstar # first grid point exception
else:
nu_T[x] = nu_Tstar*MOcoeff*Sb

# Step 1: Tentative velocity step
begin("Computing tentative velocity")
b1 = assemble(L1)
[bc.apply(A1) for bc in bcu]
solve(A1, u_.vector() , b1 , "gmres", "default")
end()
# Pressure correction
begin("Computing pressure correction")
b2 = assemble(L2)
solve(A2, p_.vector(), b2)
end()
#Velocity correction
begin("Computing velocity correction")
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u_.vector(), b3, "gmres","default")
end()
# Temperature
begin("Computing temperature")
b4 = assemble(L4)
solve(A4, theta.vector(), b4, "gmres", "default")
end()

#SGS model
begin("Computing turbulent energy")
b5 = assemble(L5)
solve(A5, e.vector(),b5, "gmres", "default")
end()

# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
theta_n.assign(theta)
e_n.assign(e)
t += dt
​

Community: FEniCS Project
If you were to include the plane average of a function within the variational form, each degree of freedom of the corresponding system of equations would depend on every other degree of freedom and your system matrix would be dense.  I give you an efficient method to calculate plane averages which you can incorporate into your system via fixed-point iteration below.
written 9 months ago by pf4d
How do I apply this if I need to use the averages in the variational form of the equations I'm assembling and solving?
written 9 months ago by Luz Imelda Pacheco

Using the above method, you have computed the horizontal $z$-varying average stored in $p(z)$.  Use this within your variational form.  You'll have to use fixed-point iteration to couple the horizontal averaging (HA) to your variational problem (VP):

0. give initial guess for values you average

repeat until converged :

1. solve VP
2. compute HA
written 9 months ago by pf4d

2
9 months ago by
I understand the logic behind posting your entire code for us to read, that is, more info is usually better than less.  However, you might consider reducing the size of your question to one simple piece that everyone here can easily digest.  That said, if you want to calculate a horizontal average of a field, consider the example:

which derives the average by integrating the function along a horizontal coordinate.  In this case, you'll have to integrate across both $x$ and $y$ and extrude the resulting integral into the domain after each integration.

That is, assume your 3D domain is $x \in [0, a]$, $y \in [0,b]$, and $z \in [0,c]$.  First you integrate across the domain numerically by solving a variational problem, beginning with the first fundamental theorem of calculus :
$v(x,y,z) = \int_0^x u(s,y,z) ds = U(x,y,z) - U(0,y,z) \hspace{10mm} (*)$
where $U$ is the anti-derivative of $v$ such that
$\frac{\partial U}{\partial x} = u(x,y,z).$
The integral at the start $x = 0$ will be zero --- meaning $U(0,y,z) = 0$ --- and so $(*)$ implies that $v(x,y,z) = U(x,y,z)$.  The BVP is therefore
$\frac{\partial v}{\partial x} = u(x,y,z), \hspace{5mm} v(0,y,z) = 0,$
and the associated variational problem consists of finding $v$ such that
$\int_{\Omega} \frac{\partial v}{\partial x} \phi \mathrm{d}\Omega = \int_{\Omega} u \phi \mathrm{d}\Omega$
for all test functions $\phi$.  After this is computed, the integral of $u$ in the $x$ direction will be located on the $x = a$ face of $v$.

To extrude the function into the domain --- so that $x$-average $w(y,z) = \frac{1}{a} v(a,y,z)$ using the $v$ you just computed is set within the entire horizontal coordinate --- just solve the BVP
$\frac{\partial w}{\partial x} = 0, \hspace{5mm} w(y,z) = \frac{1}{a} v(a,y,z).$

Next, repeat this procedure for the $y$ coordinate; i.e., first solve for the $y$-coordinate integral $q(y,z)$ from the BVP
$\frac{\partial q}{\partial y} = w(y,z), \hspace{5mm} w(0,z) = 0,$
resulting in the $y$-component integral of $v$ located on the $y = b$ face of $q$.  Finally, again solve for the extruded average $p(z) = \frac{1}{b} q(b,z)$ from
$\frac{\partial p}{\partial y} = 0, \hspace{5mm} p(b,z) = \frac{1}{b} q(b,z).$
The resulting function $p = p(z)$ is the $z$-varying horizontal average of the original function $u = u(x,y,z)$.

In summary, for each function you wish to horizontally average, you will solve four simple, well-behaved linear systems of equations; one integration and one extrusion for each of two horizontal coordinates.
Thanks for your responses! I guess I'm still not sure how to set up the second BVP to extrude back into the domain and repeat the process. Would I interpolate or project the solution at face x=a to a function space w using the same 3D mesh and then solve? how is the boundary condition enforced?
written 9 months ago by Luz Imelda Pacheco
from fenics import *

# MESHES TO INTERPOLATE HOMOGENEOUS AVERAGES
Lx, Ly, Lz = 1, 1, 1
rsl = 10
mesh3D = BoxMesh(Point(0,0,0),Point(Lx,Ly,Lz),rsl,rsl,rsl)

Q0 = FunctionSpace(mesh3D,'CG',1)
Q1 = FunctionSpace(mesh3D,'CG',1)
Q2 = FunctionSpace(mesh3D,'CG',1)
Q3 = FunctionSpace(mesh3D,'CG',1)

u    = interpolate(Expression('cos(x[0])', degree=2), Q0)
v    = TrialFunction(Q0)
phi0  = TestFunction(Q0)

# Define left side of Function Space for anti-derivative
def left(x, on_boundary):
return on_boundary and abs(x[0]) < 1e-14
def front(x, on_boundary):
return on_boundary and abs(x[1]) < 1e-14

# integral is zero on the left
bc0 = DirichletBC(Q0, 0.0, left)
bc1 = DirichletBC(Q2, 0.0, front)

# INTEGRATE DOMAIN IN X-DIRECTION (obtain v)
a0      = v.dx(0) * phi0 * dx
L0      = u * phi0 * dx
v      = Function(Q0)
solve(a0 == L0, v, bc0)

# v has the integral of the the domain in the x-direction
# EXTRUDE X-AVERAGE (obtain w)
w = TrialFunction(Q1)
phi1= TestFunction(Q1)
w = interpolate((1/Lx)*v,Q1)
f = 0.0

a1      = w.dx(0) * phi1 * dx
L1      = f * phi1 * dx
w      = Function(Q1)
solve(a1 == L1, w, bc0)

# INTEGRATE DOMAIN IN Y-DIRECTION (obtain q)
q = TrialFunction(Q2)
phi2 = TestFunction(Q2)
a2      = q.dx(1) * phi2 * dx
L2      = w * phi2 * dx
q = Function(Q2)
solve(a2 == L2, q, bc1)

# EXTRUDE Y-AVERAGE (obtain p)
p =TrialFunction(Q3)
phi3 = TestFunction(Q3)
p = interpolate((1/Ly)*q,Q3)

a3      = p.dx(1) * phi3 * dx
L3      = f * phi1 * dx
p      = Function(Q3)
solve(a3 == L3, p, bc1)

This is what I've tried so far, am I on the right track?

written 9 months ago by Luz Imelda Pacheco
I think I've been making some progress, I'm currently stuck on how to extract the values from the face and setup w
written 9 months ago by Luz Imelda Pacheco