Merge pull request #169 from Pressio/burgers2d_outflow

Adds non-periodic ("outflow") case for 2D Burgers'

Author: fnrizzi

.github/workflows/test-all.yml

@@ -46,8 +46,11 @@ jobs:
 
       - name: Install packages
         run: |
+          echo "deb http://security.ubuntu.com/ubuntu focal-security main" | sudo tee /etc/apt/sources.list.d/focal-security.list
           sudo apt-get update
+          sudo apt-get install libssl1.1
           sudo apt-get install -y --install-suggests $APT_PACKAGES
+          sudo rm /etc/apt/sources.list.d/focal-security.list
 
       - name: Install CMake
         run: |
@@ -93,7 +96,7 @@ jobs:
             -DCMAKE_BUILD_TYPE:STRING=${{ matrix.config.mode }} \
             -DPRESSIODEMOAPPS_ENABLE_TESTS:BOOL=ON \
             -DPRESSIODEMOAPPS_ENABLE_BINDINGS:BOOL=OFF \
-            -DCMAKE_CXX_FLAGS="-Wall -Werror" \
+            -DCMAKE_CXX_FLAGS="-Wall -Wno-unused-but-set-variable -Werror" \
             -B $BUILD_DIR -S $DEMOAPPS_HOME
 
       # NOTE: -j2 caused g++ crash (out of memory)

docs/source/burgers_2d_periodic.rst

@@ -59,12 +59,12 @@ C++ synopsis
 
      const auto meshObj = pda::load_cellcentered_uniform_mesh_eigen("path-to-mesh");
 
+     const auto probId  = pda::AdvectionDiffusion2d::BurgersPeriodic;
      const auto inviscidScheme = pda::InviscidFluxReconstruction::FirstOrder; // or Weno3, Weno5
      const auto viscousScheme  = pda::ViscousFluxReconstruction::FirstOrder;  // must be FirstOrder
 
      // A. constructor for problem using default values
      {
-       const auto probId = pda::AdvectionDiffusion2d::BurgersPeriodic;
        auto problem = pda::create_problem_eigen(meshObj, probId, inviscidScheme, viscousScheme);
      }
 
@@ -76,9 +76,9 @@ C++ synopsis
        const auto D      = /* something */;
        const auto x0     = /* something */;
        const auto y0     = /* something */;
-       auto problem = pda::create_periodic_burgers_2d_problem_eigen(meshObj, inviscidScheme,
-								    viscousScheme, alpha,
-								    delta, D, x0, y0)
+       auto problem = pda::create_problem_eigen(meshObj, probId, inviscidScheme,
+						viscousScheme, alpha,
+						delta, D, x0, y0)
      }
    }
 
@@ -91,11 +91,11 @@ Python synopsis
 
    meshObj = pda.load_cellcentered_uniform_mesh("path-to-mesh")
 
+   probId  = pda.AdvectionDiffusion2d.BurgersPeriodic
    inviscidScheme = pda.InviscidFluxReconstruction.FirstOrder; # or Weno3, Weno5
    viscousScheme  = pda.ViscousFluxReconstruction.FirstOrder;  # must be FirstOrder
 
    # A. constructor for problem using default values
-   probId  = pda.AdvectionDiffusion2d.BurgersPeriodic
    problem = pda.create_problem(meshObj, probId, inviscidScheme, viscousScheme)
 
    # B. setting custom coefficients
@@ -104,9 +104,9 @@ Python synopsis
    D      = # something
    x0     = # something
    y0     = # something
-   problem = pda.create_periodic_burgers_2d_problem(meshObj, inviscidScheme,
-                                                    viscousScheme, alpha,
-                                                    delta, D, x0, y0)
+   problem = pda.create_burgers_2d_problem(meshObj, probId, inviscidScheme,
+					   viscousScheme, alpha,
+					   delta, D, x0, y0)
 
 
 Sample Plot

include/pressiodemoapps/advection_diffusion2d.hpp

@@ -53,6 +53,7 @@
 #include "./container_fncs/all.hpp"
 #include "./mesh.hpp"
 #include "./schemes_info.hpp"
+#include "./impl/custom_bc_holder.hpp"
 #include "./adapter_cpp.hpp"
 #include "./adapter_py.hpp"
 
@@ -62,7 +63,8 @@ namespace pressiodemoapps{
 // enums identifying the problems
 // ----------------------------------------------------------
 enum class AdvectionDiffusion2d{
-  BurgersPeriodic
+  BurgersPeriodic,
+  BurgersOutflow
 };
 
 }//end namespace pressiodemoapps
@@ -93,26 +95,17 @@ create_problem_eigen
 {
 
   using scalar_t = typename mesh_t::scalar_t;
-  if (problemEnum == ::pressiodemoapps::AdvectionDiffusion2d::BurgersPeriodic)
-  {
-    // default parameters
-    const auto icPulseMagnitude = static_cast<scalar_t>(0.5);
-    const auto icSpread = static_cast<scalar_t>(0.15);
-    const auto diffusion = static_cast<scalar_t>(0.00001);
-    const auto icCenterX = static_cast<scalar_t>(0.0);
-    const auto icCenterY = static_cast<scalar_t>(-0.2);
-
-    return RetType(impladvdiff2d::TagBurgersPeriodic{},
-		   meshObj,
-		   inviscidFluxRecEnum,
-		   InviscidFluxScheme::Rusanov,
-		   viscFluxRecEnum,
-		   icPulseMagnitude, icSpread, diffusion,
-		   icCenterX, icCenterY);
-  }
-  else{
-    throw std::runtime_error("advection-diffusion2d: invalid problem enum");
-  }
+  auto defaultPhsParams = impladvdiff2d::defaultPhysicalParams<scalar_t>;
+  const auto defaultICParams = (problemEnum == AdvectionDiffusion2d::BurgersPeriodic || problemEnum == AdvectionDiffusion2d::BurgersOutflow)
+    ? impladvdiff2d::defaultInitCondParams<scalar_t> : std::vector<scalar_t>();
+
+  return RetType(meshObj,
+  		 problemEnum,
+		 inviscidFluxRecEnum,
+		 InviscidFluxScheme::Rusanov,
+		 viscFluxRecEnum,
+		 defaultICParams,
+		 defaultPhsParams);
 }
 
 // ----------------------------------------------------------
@@ -126,28 +119,147 @@ template<
 RetType
 #if defined PRESSIODEMOAPPS_ENABLE_BINDINGS
 // bindings need unique naming or we get error associated with overloads
-create_periodic_burgers_2d_problem_ov1_for_py
+create_burgers_2d_problem_ov1_for_py
 #else
-create_periodic_burgers_2d_problem_eigen
+create_problem_eigen
 #endif
 (const mesh_t & meshObj,
+ AdvectionDiffusion2d problemEnum,
  InviscidFluxReconstruction inviscidFluxRecEnum,
  ViscousFluxReconstruction viscFluxRecEnum,
- typename mesh_t::scalar_t icPulseMagnitude,
- typename mesh_t::scalar_t icSpread,
- typename mesh_t::scalar_t diffusion,
- typename mesh_t::scalar_t icCenterX,
- typename mesh_t::scalar_t icCenterY)
+ const std::unordered_map<std::string, typename mesh_t::scalar_t> & userParams)
 {
 
-  return RetType(impladvdiff2d::TagBurgersPeriodic{},
-		 meshObj,
+  // only one IC for now
+  int icFlag = 1;
+
+  using scalar_t = typename mesh_t::scalar_t;
+  auto physParamsVec = impladvdiff2d::defaultPhysicalParams<scalar_t>;
+  auto icParamsVec   = impladvdiff2d::defaultInitCondParams<scalar_t>;
+  impladvdiff2d::replace_params_from_map_if_present(physParamsVec, icParamsVec, problemEnum, icFlag, userParams);
+
+  return RetType(meshObj,
+		 problemEnum,
+		 inviscidFluxRecEnum,
+		 InviscidFluxScheme::Rusanov,
+		 viscFluxRecEnum,
+		 icParamsVec,
+		 physParamsVec);
+}
+
+// repeat above with default viscous flux scheme and icFlag to fit other problems
+template<
+  class mesh_t,
+  class RetType = PublicProblemEigenMixinCpp<impladvdiff2d::EigenApp<mesh_t>>
+  >
+RetType
+create_problem_eigen(const mesh_t & meshObj,
+		     AdvectionDiffusion2d problemEnum,
+		     InviscidFluxReconstruction inviscidFluxRecEnum,
+		     int icFlag,
+		     const std::unordered_map<std::string, typename mesh_t::scalar_t> & userParams)
+{
+
+  // defaults
+  ViscousFluxReconstruction viscFluxRecEnum = ViscousFluxReconstruction::FirstOrder;
+  icFlag = 1;
+
+  using scalar_t = typename mesh_t::scalar_t;
+  auto physParamsVec = impladvdiff2d::defaultPhysicalParams<scalar_t>;
+  auto icParamsVec   = impladvdiff2d::defaultInitCondParams<scalar_t>;
+  impladvdiff2d::replace_params_from_map_if_present(physParamsVec, icParamsVec, problemEnum, icFlag, userParams);
+
+  return RetType(meshObj,
+		 problemEnum,
 		 inviscidFluxRecEnum,
 		 InviscidFluxScheme::Rusanov,
 		 viscFluxRecEnum,
-		 icPulseMagnitude, icSpread, diffusion,
-		 icCenterX, icCenterY);
+		 icParamsVec,
+		 physParamsVec);
 }
 
+//
+// custom BCs
+//
+
+#if !defined PRESSIODEMOAPPS_ENABLE_BINDINGS
+template<class mesh_t, class BCsFuncL, class BCsFuncF, class BCsFuncR, class BCsFuncB>
+auto create_problem_eigen(const mesh_t & meshObj,
+			  AdvectionDiffusion2d problemEnum,
+			  InviscidFluxReconstruction inviscidFluxRecEnum,
+			  BCsFuncL && BCsLeft,
+			  BCsFuncF && BCsFront,
+			  BCsFuncR && BCsRight,
+			  BCsFuncB && BCsBack,
+			  const int icFlag = 1)
+{
+
+  // default viscous flux order
+  // TODO: generalize
+  ViscousFluxReconstruction viscFluxRecEnum = ViscousFluxReconstruction::FirstOrder;
+
+  using BCFunctorsHolderType = impl::CustomBCsHolder<BCsFuncL, BCsFuncF, BCsFuncR, BCsFuncB>;
+  BCFunctorsHolderType bcFuncs(std::forward<BCsFuncL>(BCsLeft),
+			       std::forward<BCsFuncF>(BCsFront),
+			       std::forward<BCsFuncR>(BCsRight),
+			       std::forward<BCsFuncB>(BCsBack));
+
+  using scalar_t = typename mesh_t::scalar_t;
+  const auto defaultPhsParams = impladvdiff2d::defaultPhysicalParams<scalar_t>;
+  const auto defaultICParams = impladvdiff2d::defaultInitCondParams<scalar_t>;
+
+  using return_type = PublicProblemEigenMixinCpp<impladvdiff2d::EigenApp<mesh_t, BCFunctorsHolderType>>;
+  return return_type(meshObj,
+		     problemEnum,
+		     inviscidFluxRecEnum,
+		     InviscidFluxScheme::Rusanov,
+		     viscFluxRecEnum,
+		     std::move(bcFuncs),
+		     defaultPhsParams,
+		     defaultICParams);
+
+}
+#endif
+
+#if !defined PRESSIODEMOAPPS_ENABLE_BINDINGS
+template<class mesh_t, class BCsFuncL, class BCsFuncF, class BCsFuncR, class BCsFuncB>
+auto create_problem_eigen(const mesh_t & meshObj,
+			  AdvectionDiffusion2d problemEnum,
+			  InviscidFluxReconstruction inviscidFluxRecEnum,
+			  BCsFuncL && BCsLeft,
+			  BCsFuncF && BCsFront,
+			  BCsFuncR && BCsRight,
+			  BCsFuncB && BCsBack,
+			  const int icFlag,
+			  const std::unordered_map<std::string, typename mesh_t::scalar_t> & userParams)
+{
+
+  // default viscous flux order
+  // TODO: generalize
+  ViscousFluxReconstruction viscFluxRecEnum = ViscousFluxReconstruction::FirstOrder;
+
+  using BCFunctorsHolderType = impl::CustomBCsHolder<BCsFuncL, BCsFuncF, BCsFuncR, BCsFuncB>;
+  BCFunctorsHolderType bcFuncs(std::forward<BCsFuncL>(BCsLeft),
+			       std::forward<BCsFuncF>(BCsFront),
+			       std::forward<BCsFuncR>(BCsRight),
+			       std::forward<BCsFuncB>(BCsBack));
+
+  using scalar_t = typename mesh_t::scalar_t;
+  auto physParamsVec = impladvdiff2d::defaultPhysicalParams<scalar_t>;
+  auto icParamsVec   = impladvdiff2d::defaultInitCondParams<scalar_t>;
+  impladvdiff2d::replace_params_from_map_if_present(physParamsVec, icParamsVec, problemEnum, icFlag, userParams);
+
+  using return_type = PublicProblemEigenMixinCpp<impladvdiff2d::EigenApp<mesh_t, BCFunctorsHolderType>>;
+  return return_type(meshObj,
+		     problemEnum,
+		     inviscidFluxRecEnum,
+		     InviscidFluxScheme::Rusanov,
+		     viscFluxRecEnum,
+		     std::move(bcFuncs),
+		     icParamsVec,
+		     physParamsVec);
+}
+#endif
+
 } //end namespace pressiodemoapps
 #endif

include/pressiodemoapps/impl/advection_diffusion_2d_flux_mixin.hpp

@@ -93,13 +93,10 @@ struct ComputeDirectionalFluxValues : Parent
     switch(m_fluxEnum)
       {
       case ::pressiodemoapps::InviscidFluxScheme::Rusanov:
-	if (m_probEn == ::pressiodemoapps::AdvectionDiffusion2d::BurgersPeriodic)
-	{
-	  burgers_rusanov_flux_2d<scalar_type>(m_fluxL, uMinusHalfNeg,
-					       uMinusHalfPos, m_normal);
-	  burgers_rusanov_flux_2d<scalar_type>(m_fluxR, uPlusHalfNeg,
-					       uPlusHalfPos, m_normal);
-	}
+	burgers_rusanov_flux_2d<scalar_type>(m_fluxL, uMinusHalfNeg,
+					     uMinusHalfPos, m_normal);
+	burgers_rusanov_flux_2d<scalar_type>(m_fluxR, uPlusHalfNeg,
+					     uPlusHalfPos, m_normal);
 	break;
       }
   }
@@ -154,15 +151,12 @@ struct ComputeDirectionalFluxJacobians : Parent
       {
 
       case ::pressiodemoapps::InviscidFluxScheme::Rusanov:
-	if (m_probEn == ::pressiodemoapps::AdvectionDiffusion2d::BurgersPeriodic)
-	{
-	  burgers_rusanov_flux_jacobian_2d<scalar_type>(m_fluxJacLNeg, m_fluxJacLPos,
-							uMinusHalfNeg, uMinusHalfPos,
-							m_normal);
-	  burgers_rusanov_flux_jacobian_2d<scalar_type>(m_fluxJacRNeg, m_fluxJacRPos,
-							uPlusHalfNeg, uPlusHalfPos,
-							m_normal);
-	}
+	burgers_rusanov_flux_jacobian_2d<scalar_type>(m_fluxJacLNeg, m_fluxJacLPos,
+						      uMinusHalfNeg, uMinusHalfPos,
+						      m_normal);
+	burgers_rusanov_flux_jacobian_2d<scalar_type>(m_fluxJacRNeg, m_fluxJacRPos,
+						      uPlusHalfNeg, uPlusHalfPos,
+						      m_normal);
 	break;
       }
   }
@@ -223,18 +217,15 @@ struct ComputeDirectionalFluxValuesAndJacobians : Parent
     switch(m_fluxEnum)
       {
       case ::pressiodemoapps::InviscidFluxScheme::Rusanov:
-	if (m_probEn == ::pressiodemoapps::AdvectionDiffusion2d::BurgersPeriodic)
-        {
-	  burgers_rusanov_flux_2d<scalar_type>(m_fluxL, uMinusHalfNeg, uMinusHalfPos, m_normal);
-	  burgers_rusanov_flux_2d<scalar_type>(m_fluxR, uPlusHalfNeg, uPlusHalfPos, m_normal);
-
-	  burgers_rusanov_flux_jacobian_2d<scalar_type>(m_fluxJacLNeg, m_fluxJacLPos,
-							uMinusHalfNeg, uMinusHalfPos,
-							m_normal);
-	  burgers_rusanov_flux_jacobian_2d<scalar_type>(m_fluxJacRNeg, m_fluxJacRPos,
-							uPlusHalfNeg, uPlusHalfPos,
-							m_normal);
-	}
+	burgers_rusanov_flux_2d<scalar_type>(m_fluxL, uMinusHalfNeg, uMinusHalfPos, m_normal);
+	burgers_rusanov_flux_2d<scalar_type>(m_fluxR, uPlusHalfNeg, uPlusHalfPos, m_normal);
+
+	burgers_rusanov_flux_jacobian_2d<scalar_type>(m_fluxJacLNeg, m_fluxJacLPos,
+						      uMinusHalfNeg, uMinusHalfPos,
+						      m_normal);
+	burgers_rusanov_flux_jacobian_2d<scalar_type>(m_fluxJacRNeg, m_fluxJacRPos,
+						      uPlusHalfNeg, uPlusHalfPos,
+						      m_normal);
 	break;
       }
   }