Commit 1c582a25 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts Committed by Philippe Delandmeter
Browse files

rewrite 3d open boundaries

parent bee953fd
......@@ -212,13 +212,23 @@ class Slim3d_equations :
self._boundary_coast.append(tag)
def set_boundary_open(self, physical_tags, eta=None, u=None, v=None, salinity=None, temperature=None, flux=None, transport=False):
groups2 = self._slimSolver.groups2d
groups3 = self._slimSolver.groups3d
if (u and not v) or (v and not u):
dgpy.Msg.Fatal('Error. Either u and v or none of them must be set in open boundary conditions')
uvF = functorMerge([slim_private._load_function(u, groups3), slim_private._load_function(v, groups3)]) if u else None
uv2F = functorMerge([slim_private._load_function(u, groups2), slim_private._load_function(u, groups2)]) if u else None
etaF = slim_private._load_function(eta, groups3) if eta else None
eta2F = slim_private._load_function(eta,groups2) if eta else None
flowF = slim_private._load_function(flux, groups3) if flux else None
flow2F = slim_private._load_function(flux, groups2) if flux else None
SF = slim_private._load_function(salinity, groups3) if salinity else None
TF = slim_private._load_function(temperature, groups3) if temperature else None
if slim_private._is_string(physical_tags):
openBnd = slim_private.OpenBnd_3d(physical_tags, eta, u, v, salinity, temperature, flux, transport)
physical_tags = [physical_tags]
for tag in physical_tags:
openBnd = slim3d_private.OpenBoundary(self._slimSolver, tag, etaF, eta2F, uvF, uv2F, SF, TF, flowF, flow2F, transport)
self._boundary_open.append(openBnd)
else:
for tag in physical_tags:
openBnd = slim_private.OpenBnd_3d(tag, eta, u, v, salinity, temperature, flux, transport)
self._boundary_open.append(openBnd)
class Loop:
"""Temporal solver"""
......
......@@ -5,6 +5,108 @@ import dgpy
from dgpy.scripts import slim_private
import numpy
class OpenBoundary : # etaF can be None
def __init__(self, solver, tag, etaF, eta2F, uvF, uv2F, salinityF, temperatureF, flowF, flow2F, transport) :
self.solver = solver
self.tag = tag
self.transport = transport
self.uvF = uvF
self.uv2F = uv2F
self.etaF = etaF
self.eta2F = eta2F
self.SF = salinityF
self.TF = temperatureF
self.flowF = flowF
self.flow2F = flow2F
self.horizontal_momentum_flux_f = dgpy.functorNumpy(self.horizontal_momentum_flux)
self.sw2d_flux_f = dgpy.functorNumpy(self.sw2d_flux)
self.w_flux_f = dgpy.functorNumpy(self.w_flux)
self.s_flux_f = dgpy.functorNumpy(lambda cm : self.tracer_flux(cm, "S"))
self.t_flux_f = dgpy.functorNumpy(lambda cm : self.tracer_flux(cm, "T"))
def setup(self) :
groups2 = self.solver.groups2d
self.hF = self.solver.functions.bathFunc2d
if self.flowF :
self.section = dgpy.dgFunctionIntegratorInterface(groups2, hF).compute(self.tag)
def compute_uv_ext(self, cm, uvin, etain) :
dim = cm.elementGroup().getDimUVW()
if self.uvF and not self.transport :
return cm.get(self.uvF) if dim == 3 else cm.get(self.uv2F)
if self.etaF :
eta = cm.get(self.etaF) if dim == 3 else cm.get(self.eta2F)
else :
eta = etain
h = cm.get(self.hF)
if self.uvF and self.transport :
UV = cm.get(self.uvF) if dim == 3 else cm.get(self.uv2F)
return UV/(h+eta)
elif self.flowF :
flowF = self.flowF if dim == 3 else self.flow2F
return -cm.get(flow)/self.section * h/(h+eta) * cm.get(dgpy.function.getNormals())
else :
n = cm.get(dgpy.function.getNormals())
un = uvin[:,[0]] * n[:,[0]] + uvin[:,[1]] * n[:,[1]]
#todo : this is not correct
return un * n
def compute_eta_ext(self, cm, uvin, etain) : #only 2d
if self.etaF :
return cm.get(self.eta2F)
else :
n = cm.get(dgpy.function.getNormals())
un = uvin[:,[0]] * n[:,[0]] + uvin[:,[1]] * n[:,[1]]
#todo : this is not correct
g = dgpy.slim3dParameters.getG()
return numpy.sqrt((cm.get(self.hF)+etain)/ g) * un
def horizontal_momentum_flux(self,cm) :
cm.pushState()
cm.setDefaultToSide0(True)
eq = self.solver.getEquations().horMomEq
#solIn = cm.get(dgpy.function.getSolution(), 0)
eta_in = cm.get(eq.getEta(), 0)
#cm.set_value(dgpy.function.getSolution(), 1,self.compute_uv_ext(cm, solIn, eta_inn))
uv2d_in = cm.get(eq.getUV2d(), 0)
cm.set_value(eq.getUV2d(), 1,self.compute_uv_ext(cm, uv2d_in, eta_in))
cm.set_value(dgpy.function.getSolution(), 1, numpy.zeros((cm.nEvaluationPoint(), 2)))
f = cm.get(eq.getInterfaceTerm0(cm.interfaceGroup().physicalTag())).copy()
cm.popState()
return f
def sw2d_flux(self,cm) :
cm.pushState()
cm.setDefaultToSide0(True)
eq = self.solver.getEquations().sw2DEq
solIn = cm.get(dgpy.function.getSolution(), 0)
sol_ext = numpy.hstack([self.compute_eta_ext(cm, solIn[:,[1,2]], solIn[:,[0]]), self.compute_uv_ext(cm, solIn[:,[1,2]], solIn[:,[0]])])
cm.set_value(dgpy.function.getSolution(), 1, sol_ext)
f = cm.get(eq.getInterfaceTerm0(cm.interfaceGroup().physicalTag())).copy()
cm.popState()
return f
def w_flux(self,cm) :
cm.pushState()
cm.setDefaultToSide0(True)
eq = self.solver.getEquations().wEq
etaIn = cm.get(eq.getEta(), 0)
uvIn = cm.get(eq.getUV(), 0)
cm.set_value(eq.getEta(), 1, self.compute_eta_ext(cm, uvIn, etaIn))
cm.set_value(eq.getUV(), 1, self.compute_uv_ext(cm, uvIn, etaIn))
f = cm.get(eq.getInterfaceTerm0(cm.interfaceGroup().physicalTag())).copy()
cm.popState()
return f
def tracer_flux(self,cm, mode) :
cm.pushState()
cm.setDefaultToSide0(True)
eq = self.solver.getEquations().SEq if mode == "S" else self.solver.getEquations().TEq
etaIn = cm.get(eq.getEta(), 0)
uvIn = cm.get(eq.getUV(), 0)
if mode == "S" :
solExt = cm.get(self.SF) if self.SF else cm.get(self.solver.functions.SInitFunc)
else :
solExt = cm.get(self.TF) if self.TF else cm.get(self.solver.functions.TInitFunc)
cm.set_value(eq.getEta(), 1, self.compute_eta_ext(cm, uvIn, etaIn))
cm.set_value(eq.getUV(), 1, self.compute_uv_ext(cm, uvIn, etaIn))
cm.set_value(dgpy.function.getSolution(), 1, solExt)
f = cm.get(eq.getInterfaceTerm0(cm.interfaceGroup().physicalTag())).copy()
cm.popState()
return f
def slim3d_setup(loop):
slimSolver = loop._slimSolver
......@@ -39,7 +141,7 @@ def slim3d_setup(loop):
if eq._wind_stress[0] == 'speed':
eq._wind_speed_x = slim_private._load_function(eq._wind_stress[1], slimSolver.groups2d)
eq._wind_speed_y = slim_private._load_function(eq._wind_stress[2], slimSolver.groups2d)
f.windFunc = functorNumpy(lambda cm : mergeUVNumpy(cm, eq.wind_speed_x, eq.wind_speed_y))
f.windFunc = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq.wind_speed_x, eq.wind_speed_y))
f.windStressFunc = dgpy.windCubicStress(f.windFunc, eq._wind_stress[3])
# TODO Fatal if no wind speed and solving sediment
elif eq._wind_stress[0] == 'stress':
......@@ -154,101 +256,16 @@ def slim3d_setup(loop):
g = eq._gravity
# Open Boundaries
def uvOpenGivenEtaNumpy(cm, u, normals) :
u = cm.get(uv)
n = cm.get(normals)
un = u[:,[0]] * n[:,[0]] + u[:,[1]] * n[:,[1]]
return un * n
def uvOpenFreeNumpy(cm, bath, eta, normals) :
c = numpy.sqrt(g/(cm.get(bath)+cm.get(eta)))
return c*cm.get(eta)*cm.get(normals)
def etaOpenGivenUVNumpy(cm, u, bath, eta, normals) :
u = cm.get(uv)
n = cm.get(normals)
un = u[:,[0]] * n[:,[0]] + u[:,[1]] * n[:,[1]]
return numpy.sqrt( (cm.get(bath)+cm.get(eta)) / g) * un
def transport2UVNumpy(cm, uv, bath, eta) :
return cm.get(uv) / (cm.get(bath) + cm.get(eta))
eta = f.eta2dFunc
uv = d.uvDof.getFunction()
uvAv2d = f.uvAv2dFunc
eq._uv_open = {}
eq._uvDev_open = {}
eq._uvAv2d_open = {}
eq._uvTransport_open = {}
eq._uvAv2dTransport_open = {}
eq._u_open = {}
eq._v_open = {}
eq._uAv2d_open = {}
eq._vAv2d_open = {}
eq._flux_section = {}
eq._flux = {}
eq._flux2d = {}
eq._eta_open = {}
eq._S_open = {}
eq._T_open = {}
eq._sw2D_open = {}
def functorMerge(allf) :
return dgpy.functorNumpy(lambda cm : numpy.hstack([cm.get(f) for f in allf]))
for openBnd in eq._boundary_open:
eq._uvDev_open[openBnd] = dgpy.functionConstant([0.,0.])
if openBnd.u and openBnd.eta:
eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._v_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups3d)
eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups2d)
if openBnd.transport:
eq._uvTransport_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._u_open[openBnd], eq._v_open[openBnd]))
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : transport2UVNumpy(cm, eq._uvTransport_open[openBnd], f.bathFunc2d, eta))
eq._uvAv2dTransport_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : transport2UVNumpy(cm, eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta))
else:
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._u_open[openBnd], eq._v_open[openBnd]))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]))
eq._eta_open[openBnd] = slim_private._load_function(openBnd.eta, slimSolver.groups2d)
elif openBnd.u:
eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._v_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups3d)
eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups2d)
if openBnd.transport:
eq._uvTransport_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._u_open[openBnd], eq._v_open[openBnd]))
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : transport2UVNumpy(cm, eq._uvTransport_open[openBnd], f.bathFunc2d, eta))
eq._uvAv2dTransport_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : transport2UVNumpy(cm, eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta))
else:
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._u_open[openBnd], eq._v_open[openBnd]))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]))
eq._eta_open[openBnd] = dgpy.functorNumpy(lambda cm : etaOpenGivenUVNumpy(cm, uvAv2d, f.bathFunc2d, eta, dgpy.function.getNormals()))
elif openBnd.eta:
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : uvOpenGivenEtaNumpy(cm, uv, dgpy.function.getNormals()))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : uvOpenGivenEtaNumpy(cm, uvAv2d, dgpy.function.getNormals()))
eq._eta_open[openBnd] = slim_private._load_function(openBnd.eta, slimSolver.groups2d)
elif openBnd.flux:
integInterface = dgpy.dgFunctionIntegratorInterface(slimSolver.groups2d, f.bathFunc2d)
sectionFM = dgpy.fullMatrixDouble(1,1)
integInterface.compute(openBnd.tag, sectionFM)
eq._flux_section[openBnd] = sectionFM(0,0)
eq._flux[openBnd] = slim_private._load_function(openBnd.flux, slimSolver.groups3d)
eq._flux2d[openBnd] = slim_private._load_function(openBnd.flux, slimSolver.groups2d)
eq._uv_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups3d, eq._flux_section[openBnd], eq._flux[openBnd], f.bathFunc2d, eta)
eq._uvAv2d_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups2d, eq._flux_section[openBnd], eq._flux2d[openBnd], f.bathFunc2d, eta)
eq._eta_open[openBnd] = eta
else:
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : uvOpenFreeNumpy(cm, f.bathFunc2d, eta, dgpy.function.getNormals()))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : uvOpenFreeNumpy(cm, f.bathFunc2d, eta, dgpy.function.getNormals()))
eq._eta_open[openBnd] = eta
if openBnd.salinity:
if openBnd.salinity == 'vertical_gradient':
eq._S_open[openBnd] = f.SInitFunc
else:
eq._S_open[openBnd] = slim_private._load_function(openBnd.salinity, slimSolver.groups3d)
if openBnd.temperature:
if openBnd.temperature == 'vertical_gradient':
eq._T_open[openBnd] = f.TInitFunc
else:
eq._T_open[openBnd] = slim_private._load_function(openBnd.temperature, slimSolver.groups3d)
eq._sw2D_open[openBnd] = dgpy.functorNumpy(lambda cm : numpy.hstack([cm.get(eq._eta_open[openBnd]), cm.get(eq._uvAv2d_open[openBnd])]))
openBnd.setup()
sw2DEq = e.sw2DEq
sw2DEq.setFrom3D(True)
......@@ -259,7 +276,7 @@ def slim3d_setup(loop):
sw2DBndWall = sw2DEq.newBoundaryWall(True)
sw2DEq.addBoundaryCondition(eq._boundary_coast, sw2DBndWall)
for openBnd in eq._boundary_open:
sw2DEq.addBoundaryCondition(openBnd.tag, sw2DEq.newOutsideValueBoundary(eq._sw2D_open[openBnd]))
sw2DEq.setBoundaryNeumann(openBnd.tag, openBnd.sw2d_flux_f)
horMomEq = e.horMomEq
#horMomEq.setLaxFriedrichsFactor(eq._LFF)
......@@ -269,10 +286,8 @@ def slim3d_setup(loop):
else :
horMomEq.addBoundaryCondition(topTags, horMomEq.newBoundarySurface(f.windStressFunc)) # zero for nonconst tracers!
for openBnd in eq._boundary_open:
toReplaceHorMomOpen = dgpy.VectorFunctorConst([dgpy.function.getSolution(), uvAv2d, eta])
replaceByHorMomOpen = dgpy.VectorFunctorConst([eq._uvDev_open[openBnd], eq._uvAv2d_open[openBnd], eq._eta_open[openBnd]])
horMomBndOpen = horMomEq.newOutsideValueBoundaryGeneric("", toReplaceHorMomOpen, replaceByHorMomOpen)
horMomEq.addBoundaryCondition(openBnd.tag, horMomBndOpen)
horMomEq.setBoundaryNeumann(openBnd.tag, openBnd.horizontal_momentum_flux_f)
horMomBndWall = horMomEq.newBoundaryWall()
horMomEq.addBoundaryCondition(eq._boundary_coast, horMomBndWall)
horMomEq.addBoundaryCondition(bottomTags, horMomBndWall)
......@@ -296,10 +311,7 @@ def slim3d_setup(loop):
wEq.setBoundary0Flux(eq._boundary_coast)
wEq.setBoundary0Flux('vertical_bottom')
for openBnd in eq._boundary_open:
toReplaceWOpen = dgpy.VectorFunctorConst([uv, eta])
replaceByWOpen = dgpy.VectorFunctorConst([eq._uv_open[openBnd], eq._eta_open[openBnd]])
wBndOpen = wEq.newOutsideValueBoundaryGeneric("", toReplaceWOpen, replaceByWOpen)
wEq.addBoundaryCondition(openBnd.tag, wBndOpen)
wEq.setBoundaryNeumann(openBnd.tag, openBnd.w_flux_f)
e.newRGradEq.setBoundarySymmetry(topTags)
e.newRGradEq.setBoundarySymmetry(bottomTags)
......@@ -317,10 +329,7 @@ def slim3d_setup(loop):
e.SEq.setBoundary0Flux(eq._boundary_coast)
e.SEq.setBoundary0Flux('vertical_bottom')
for openBnd in eq._boundary_open:
toReplaceSOpen = dgpy.VectorFunctorConst([uv, eta, dgpy.function.getSolution()])
replaceBySOpen = dgpy.VectorFunctorConst([eq._uv_open[openBnd], eq._eta_open[openBnd], eq._S_open[openBnd]])
SBndOpen = e.SEq.newOutsideValueBoundaryGeneric("", toReplaceSOpen, replaceBySOpen)
e.SEq.addBoundaryCondition(openBnd.tag, SBndOpen)
e.SEq.setBoundaryNeumann(openBnd.tag, openBnd.s_flux_f)
if slimSolver.getSolveT() :
e.TEq.setLaxFriedrichsFactor(eq._LFF)
......@@ -329,10 +338,7 @@ def slim3d_setup(loop):
e.TEq.setBoundary0Flux(eq._boundary_coast)
e.TEq.setBoundary0Flux('vertical_bottom')
for openBnd in eq._boundary_open:
toReplaceTOpen = dgpy.VectorFunctorConst([uv, eta, dgpy.function.getSolution()])
replaceByTOpen = dgpy.VectorFunctorConst([eq._uv_open[openBnd], eq._eta_open[openBnd], eq._T_open[openBnd]])
TBndOpen = e.TEq.newOutsideValueBoundaryGeneric("", toReplaceTOpen, replaceByTOpen)
e.TEq.addBoundaryCondition(openBnd.tag, TBndOpen)
e.TEq.setBoundaryNeumann(openBnd.tag, openBnd.t_flux_f)
if slimSolver.getAdvectTurbulentEnergy() :
for turbEq in [e.tkeAdvEq, e.epsAdvEq] :
......
......@@ -247,19 +247,6 @@ def _findTopAndBottomTags(mesh_file_name):
fin.close()
return (topTags, bottomTags)
class OpenBnd_3d:
def __init__(self, tag, eta, u, v, salinity, temperature, flux, transport):
self.tag = tag
self.u = u
self.v = v
if (u and not v) or (v and not u):
dgpy.Msg.Fatal('Error. Either u and v or none of them must be set in open boundary conditions')
self.eta = eta
self.salinity = salinity
self.temperature = temperature
self.flux = flux
self.transport = transport
class pre_evaluator :
""" base class for objects that can be evaluated one point at a time (in mesh generation)
or on a whole region (in preprocessing)
......
......@@ -61,6 +61,7 @@ public :
void pushState();
void popState();
fullMatrix<double> &set(const functor &f, int side = -1);
void set_value(const functor &f, int side, const fullMatrix<double> &v) {set(f,side).copy(v);}
bool sideDefined() const;
int side() const;
void clear();
......
......@@ -38,9 +38,9 @@ def InitialCondition(cm):
def Velocity(cm):
v = numpy.ndarray((cm.nEvaluationPoint(), 3))
v[:, 0] = Vx
v[:, 1] = Vy
v[:, 2] = Vz
v[:,0] = Vx
v[:,1] = Vy
v[:,2] = Vz
return v
print('*** Loading the Mesh ***')
......
......@@ -24,8 +24,10 @@ class dgConservationLawSW3dContinuity : public dgConservationLawFunction {
void setup();
/**slip wall boundary */
dgBoundaryCondition *newBoundaryWall();
inline void setUV(const functor *f) { _uv = f; };
inline void setEta(const functor *f) { _eta = f; };
inline void setUV(const functor *f) { _uv = f; }
inline void setEta(const functor *f) { _eta = f; }
const functor *getEta()const {return _eta;}
const functor *getUV()const {return _uv;}
inline void setWOld(const functor *f) { _wOld = f; };
inline void setBathymetryMinCG(const functor *bathymetryMinCG) { _bathymetryMinCG = bathymetryMinCG;}
inline void setLaxFriedrichsFactor(double lff){ _laxFriedrichsFactor = lff;};
......
......@@ -31,11 +31,12 @@ class dgConservationLawSW3dMomentum : public dgConservationLawFunction {
void setUseConservativeALE(bool b) { _useConservativeALE = b; };
/**set eta function [eta] */
inline void setEta(const functor *eta){ _eta = eta;}
const functor * getEta()const{return _eta;}
/**set function for horizontal gradient of the baroclinic head [drdx drdy] */
inline void setRGrad(const functor *f){ _rGrad = f;}
/**set depth averaged uv function [uInt, vInt] */
inline void setUV2d(const functor *uv2d){ _uv2d = uv2d;}
inline const functor* getUV2d(){ return _uv2d;}
const functor *getUV2d()const {return _uv2d;}
/**set depth averaged uvGrad function [dudx, dudy, dudz, dvdx, dvdy, dvdz] */
inline void setUV2dGrad(const functor *uv2dGrad){ _uv2dGrad = uv2dGrad;}
/**set vertical velocity function [w] */
......
......@@ -50,6 +50,8 @@ class dgConservationLawSW3dTracer : public dgConservationLawFunction {
inline void setWMeshDz(const functor *dwMdz) {_dwMdz = dwMdz;};
inline void setIsoDiff(const functor *S, const functor *Kedd) {_slope = S; _Kedd = Kedd;};
inline void setGMVel(const functor *uGM) {_uGM = uGM; };
const functor *getUV()const {return _uv;}
const functor *getEta()const {return _eta;}
~dgConservationLawSW3dTracer();
/** Boundary condition that sets the in-flowing tracer concentration to a prescribed value */
void boundaryInFlux(functorCache &cm, fullMatrix<double> &val, const functor *outValue) const;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment