slim3d: boundary conditions for IMEX equations

dgpy/scripts/slim3d_private.py

@@ -67,7 +67,10 @@ class OpenBoundary : # etaF can be None
         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_in))
+        uv_ext = self.compute_uv_ext(cm, solIn, eta_in)
+        cm.set_value(dgpy.function.getSolution(), 1, uv_ext)
+        if eq.isImplicitVertical():
+            cm.set_value(eq.getExplicitSolution(), 1, uv_ext)
         uv2d_in = cm.get(eq.getUV2d(), 0)
         cm.set_value(eq.getUV2d(), 1,self.compute_uv_ext(cm, uv2d_in, eta_in))
         f = cm.get(eq.getInterfaceTerm0("")).copy()
@@ -107,6 +110,8 @@ class OpenBoundary : # etaF can be None
         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)
+        if eq.isImplicitVertical():
+            cm.set_value(eq.getExplicitSolution(), 1, solExt)
         f = cm.get(eq.getInterfaceTerm0("")).copy()
         cm.popState()
         return f

modules/slim3d/dgConservationLawSW3dTracer.h

@@ -15,9 +15,11 @@ class dgConservationLawSW3dIMEX : public dgConservationLawFunction {
   inline void setW(const functor *f) {_w = f;};
   const functor *getUV()const {return _uv;}
   const functor *getEta()const {return _eta;}
+  const functor *getExplicitSolution() const {return _explicitSol;}
   const functor *getW()const {return _w;}
   void setFullExplicit();
   void setImplicitVertical(const functor *explicitSol, const functor *explicitSolGrad);
+  bool isImplicitVertical() {return (_explicitSol != NULL);};
 };
 
 class dgConservationLawSW3dTracer : public dgConservationLawSW3dIMEX {

modules/slim3d/slim3dEquations.cpp

@@ -690,9 +690,9 @@ slim3dSolverEquations::~slim3dSolverEquations()
   delete TEq;
   delete tkeAdvEq;
   delete epsAdvEq;
-  delete depthIntegratorRho;
+  //delete depthIntegratorRho;
   delete sw2DEq;
-  delete uvHorDivEq;
+  //delete uvHorDivEq;
   if ( _solver->getSolveUVGrad())
     delete uvGradEq;
   if ( _solver->getSolveGMVel() ) {

modules/slim3d/slim3dEquations.h

@@ -109,7 +109,7 @@ public:
   dgConservationLawSW3dMomentum *horMomEq;
   /** 2D shallow water equations */
   dgConservationLawShallowWater2d *sw2DEq;
-  dgUVHorDivergence *uvHorDivEq;
+  //dgUVHorDivergence *uvHorDivEq;
   /** 3d continuity equation for computing w (vertical velocity) */
   dgConservationLawSW3dContinuity *wEq;
 #ifdef HAVE_GOTM
@@ -127,7 +127,7 @@ public:
   /** 3d transport equation for turbulence dissipation rate */
   dgConservationLawSW3dTracer *epsAdvEq;
   /** Vertical integrator for integrating r (baroclinic head) */
-  dgSW3dDepthIntegrate *depthIntegratorRho;
+  //dgSW3dDepthIntegrate *depthIntegratorRho;
   /** new way to compute rGRad */
   dgConservationLawSW3dPressure *newRGradEq;