Major updates. Basically, I brought in all changes from the current (hidden from public) branch. The Taylor Green vortex test now does a great job even without periodic BCs. More testing needed. Also, there's a major speed and bug fix in MappedAMRPoissonOp.

Author: EdSantilli

src/BCutil/BGScalarProfiles.H

@@ -0,0 +1,136 @@
+/*******************************************************************************
+ *  SOMAR - Stratified Ocean Model with Adaptive Refinement
+ *  Developed by Ed Santilli & Alberto Scotti
+ *  Copyright (C) 2016
+ *      Philadelphia University &
+ *      University of North Carolina at Chapel Hill
+ *
+ *  This library is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU Lesser General Public
+ *  License as published by the Free Software Foundation; either
+ *  version 2.1 of the License, or (at your option) any later version.
+ *
+ *  This library is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  Lesser General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public
+ *  License along with this library; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301
+ *  USA
+ *
+ *  For up-to-date contact information, please visit the repository homepage,
+ *  https://github.com/somarhub.
+ ******************************************************************************/
+
+#ifndef ___BGScalarProfiles_H__INCLUDED___
+#define ___BGScalarProfiles_H__INCLUDED___
+
+#include "EllipticBCInterface.H"
+
+
+// -----------------------------------------------------------------------------
+// -----------------------------------------------------------------------------
+class LinearBGScalarProfile: public EllipticBCValueClass
+{
+public:
+    LinearBGScalarProfile (Real a_Nsq)
+    : m_Nsq(a_Nsq)
+    {;}
+
+    virtual ~LinearBGScalarProfile ()
+    {;}
+
+    virtual void operator() (Real*           a_pos,
+                             int*            a_dir,
+                             Side::LoHiSide* a_side,
+                             Real*           a_value,
+                             const RealVect& a_dx,
+                             Real            a_time) const
+    {
+        a_value[0] = -m_Nsq * a_pos[CH_SPACEDIM-1];
+    }
+
+protected:
+    Real m_Nsq;
+};
+
+
+// -----------------------------------------------------------------------------
+// -----------------------------------------------------------------------------
+class QuadraticBGScalarProfile: public EllipticBCValueClass
+{
+public:
+    QuadraticBGScalarProfile (Real a_NsqBottom,
+                              Real a_NsqTop)
+    : m_NsqBottom(a_NsqBottom),
+      m_NsqTop(a_NsqTop)
+    {
+        LevelGeometry::staticDefine();
+        m_H = LevelGeometry::getDomainLength(CH_SPACEDIM-1);
+        m_zCoeff = (m_NsqBottom - m_NsqTop) / m_H;
+    }
+
+    virtual ~QuadraticBGScalarProfile ()
+    {;}
+
+    virtual void operator() (Real*           a_pos,
+                             int*            a_dir,
+                             Side::LoHiSide* a_side,
+                             Real*           a_value,
+                             const RealVect& a_dx,
+                             Real            a_time) const
+    {
+        Real z = a_pos[CH_SPACEDIM-1];
+        Real Nsq = m_NsqBottom - m_zCoeff * z;
+        a_value[0] = -Nsq * z;
+    }
+
+protected:
+    Real m_NsqBottom;
+    Real m_NsqTop;
+    Real m_H;
+    Real m_zCoeff;
+};
+
+
+// -----------------------------------------------------------------------------
+// -----------------------------------------------------------------------------
+class TanhBGScalarProfile: public EllipticBCValueClass
+{
+public:
+    TanhBGScalarProfile (Real a_z0,
+                         Real a_delta,
+                         Real a_rho0,
+                         Real a_drho)
+    : m_z0(a_z0),
+      m_delta(a_delta),
+      m_rho0(a_rho0),
+      m_drho(a_drho)
+    {;}
+
+    virtual ~TanhBGScalarProfile ()
+    {;}
+
+    virtual void operator() (Real*           a_pos,
+                             int*            a_dir,
+                             Side::LoHiSide* a_side,
+                             Real*           a_value,
+                             const RealVect& a_dx,
+                             Real            a_time) const
+    {
+        Real arg = (a_pos[CH_SPACEDIM-1] - m_z0) / m_delta;
+        a_value[0] = m_rho0 - m_drho*tanh(arg);
+    }
+
+protected:
+    Real m_z0;      // Center location
+    Real m_delta;   // Thickness
+    Real m_rho0;    // Median buoyancy value
+    Real m_drho;    // Max buoyancy difference
+};
+
+
+#endif //!___BGScalarProfiles_H__INCLUDED__
+

src/BCutil/PhysBCUtil.H

@@ -43,22 +43,48 @@
 class PhysBCUtil
 {
 public:
+    // Scalar info -------------------------------------------------------------
+    // This controls the number of scalars in the simulation, their names,
+    // the number of components for each scalar, etc...
+
+    // These can be used to identify a specific scalar.
+    // WARNING: Do not add scalars. SOMAR can't handle it yet.
+    struct ScalarIndex {
+        enum {
+            BUOYANCY_DEVIATION = 0,
+            _COUNT
+        };
+    };
+
+    // Basic accessors
+    static inline unsigned int getNumScalars ();
+    static inline const std::string& getName (unsigned int a_idx);
+    static inline int getNumComps (unsigned int a_idx);
+    static inline const IntVect& getGhostVect (unsigned int a_idx);
+
+
+    // Construction / destruction ----------------------------------------------
+
     // Constructor
     PhysBCUtil ();
 
     // Destructor
     virtual ~PhysBCUtil ();
 
-    // Derived classes should be thier own factories.
+    // Derived classes should be their own factories.
     virtual PhysBCUtil* newPhysBCUtil () const = 0;
 
-    // Reads BC data from file
+    // Overridable constructor
     virtual void define ();
 
+    // Performs static initialization (IO, etc).
+    static void staticDefine ();
+    static inline bool isStaticDefined ();
+
 
     // ICs / background fields -------------------------------------------------
     // These functions have defaults that just set the FABs to zero and should
-    // be overriden as needed.
+    // be overridden as needed.
 
     // Fills a FAB with the initial velocity.
     // Locations are in mapped space, but components are Cartesian.
@@ -82,7 +108,7 @@ public:
 
 
     // Utility functions -------------------------------------------------------
-    // Most of these functions just call thier FAB versions and should not
+    // Most of these functions just call their FAB versions and should not
     // need to be overridden.
 
     // Sets the velocity ICs over an entire level
@@ -115,6 +141,7 @@ public:
                                            const Real            a_time,
                                            const LevelGeometry&  a_levGeo) const;
 
+    // TODO: This should be moved to StratUtils.
     // Computes the Brunt–Väisälä frequency on a single grid.
     // All ins and outs must be CC.
     // a_bFAB is the background buoyancy field.
@@ -127,6 +154,7 @@ public:
                              const Box&       a_destBox,
                              const RealVect&  a_dx) const;
 
+    // TODO: This should be moved to StratUtils.
     // Computes the Brunt–Väisälä frequency over an entire level.
     // a_Nsq must be CC.
     virtual void computeNSq (LevelData<FArrayBox>& a_Nsq,
@@ -139,8 +167,8 @@ public:
     // This is in case the BC's have an effect on the timestep.
     // Pass in currently computed dt, along with the cfl and dx. If the effect
     // of the BCs requires a decreased timestep, then the newly reduced timestep
-    // is returned. In the default case, this just returns a_dt back; however,
-    // derived classes may actually have an effect.
+    // is returned. In the default case, this just returns a_dt back. Derived
+    // classes may need this.
     virtual void computeBoundaryDt (Real&                a_dt,
                                     const Real           a_cfl,
                                     const LevelGeometry& a_levGeo) const;
@@ -201,7 +229,7 @@ public:
     // Sets ghosts needed to calculate the viscous source term nu.L[vel]
     virtual Tuple<BCMethodHolder, SpaceDim> viscousSourceFuncBC () const;
 
-    // Used in single-compoent velocity TGA solves
+    // Used in single-component velocity TGA solves
     virtual BCMethodHolder viscousSolveFuncBC (int a_dir) const;
 
     // returns single-component BC for viscous refluxing solves
@@ -319,10 +347,13 @@ protected:
     // Member variables and utility functions ----------------------------------
 
     // Have we read params from the input file yet?
-    static bool s_staticMembersSet;
+    static bool m_isStaticDefined;
 
     static RealVect s_domLength;
+
     static bool s_useBackgroundScalar;
+    static int s_bgScalarProfile;
+    static RefCountedPtr<EllipticBCValueClass> s_stdProfilePtr;
 
     static RealVect s_tidalU0;
     static Real s_tidalOmega;
@@ -331,10 +362,61 @@ protected:
     static bool s_useSpongeLayer;
     static Real s_spongeWidth[CH_SPACEDIM][2];
     static Real s_spongeDtMult[CH_SPACEDIM][2];
+
+    // Scalar metadata
+    struct ScalarMetaData
+    {
+        std::string name;
+        int         numComps;
+        IntVect     ghostVect;
+    };
+    static const ScalarMetaData s_scalarMetaData[ScalarIndex::_COUNT];
 };
 
 
 
+// -----------------------------------------------------------------------------
+// Simple accessor
+// -----------------------------------------------------------------------------
+inline bool PhysBCUtil::isStaticDefined ()
+{
+    return m_isStaticDefined;
+}
+
+
+// -----------------------------------------------------------------------------
+// Simple accessor
+// -----------------------------------------------------------------------------
+inline unsigned int PhysBCUtil::getNumScalars ()
+{
+    return ScalarIndex::_COUNT;
+}
+
+
+// -----------------------------------------------------------------------------
+// Simple accessor
+// -----------------------------------------------------------------------------
+inline const std::string& PhysBCUtil::getName (unsigned int a_idx) {
+    return s_scalarMetaData[a_idx].name;
+}
+
+
+// -----------------------------------------------------------------------------
+// Simple accessor
+// -----------------------------------------------------------------------------
+inline int PhysBCUtil::getNumComps (unsigned int a_idx) {
+    return s_scalarMetaData[a_idx].numComps;
+}
+
+
+// -----------------------------------------------------------------------------
+// Simple accessor
+// -----------------------------------------------------------------------------
+inline const IntVect& PhysBCUtil::getGhostVect (unsigned int a_idx) {
+    return s_scalarMetaData[a_idx].ghostVect;
+}
+
+
 // -----------------------------------------------------------------------------
 // Simple accessor
 // -----------------------------------------------------------------------------