diff --git a/.github/workflows/test-all.yml b/.github/workflows/test-all.yml index f81048e2..7e6ef5dd 100644 --- a/.github/workflows/test-all.yml +++ b/.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) diff --git a/docs/source/burgers_2d_periodic.rst b/docs/source/burgers_2d_periodic.rst index 241bc55d..f1764c14 100644 --- a/docs/source/burgers_2d_periodic.rst +++ b/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 diff --git a/include/pressiodemoapps/advection_diffusion2d.hpp b/include/pressiodemoapps/advection_diffusion2d.hpp index c7481f3e..3892cacb 100644 --- a/include/pressiodemoapps/advection_diffusion2d.hpp +++ b/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(0.5); - const auto icSpread = static_cast(0.15); - const auto diffusion = static_cast(0.00001); - const auto icCenterX = static_cast(0.0); - const auto icCenterY = static_cast(-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; + const auto defaultICParams = (problemEnum == AdvectionDiffusion2d::BurgersPeriodic || problemEnum == AdvectionDiffusion2d::BurgersOutflow) + ? impladvdiff2d::defaultInitCondParams : std::vector(); + + 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 & 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; + auto icParamsVec = impladvdiff2d::defaultInitCondParams; + 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> + > +RetType +create_problem_eigen(const mesh_t & meshObj, + AdvectionDiffusion2d problemEnum, + InviscidFluxReconstruction inviscidFluxRecEnum, + int icFlag, + const std::unordered_map & userParams) +{ + + // defaults + ViscousFluxReconstruction viscFluxRecEnum = ViscousFluxReconstruction::FirstOrder; + icFlag = 1; + + using scalar_t = typename mesh_t::scalar_t; + auto physParamsVec = impladvdiff2d::defaultPhysicalParams; + auto icParamsVec = impladvdiff2d::defaultInitCondParams; + 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 +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; + BCFunctorsHolderType bcFuncs(std::forward(BCsLeft), + std::forward(BCsFront), + std::forward(BCsRight), + std::forward(BCsBack)); + + using scalar_t = typename mesh_t::scalar_t; + const auto defaultPhsParams = impladvdiff2d::defaultPhysicalParams; + const auto defaultICParams = impladvdiff2d::defaultInitCondParams; + + using return_type = PublicProblemEigenMixinCpp>; + return return_type(meshObj, + problemEnum, + inviscidFluxRecEnum, + InviscidFluxScheme::Rusanov, + viscFluxRecEnum, + std::move(bcFuncs), + defaultPhsParams, + defaultICParams); + +} +#endif + +#if !defined PRESSIODEMOAPPS_ENABLE_BINDINGS +template +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 & userParams) +{ + + // default viscous flux order + // TODO: generalize + ViscousFluxReconstruction viscFluxRecEnum = ViscousFluxReconstruction::FirstOrder; + + using BCFunctorsHolderType = impl::CustomBCsHolder; + BCFunctorsHolderType bcFuncs(std::forward(BCsLeft), + std::forward(BCsFront), + std::forward(BCsRight), + std::forward(BCsBack)); + + using scalar_t = typename mesh_t::scalar_t; + auto physParamsVec = impladvdiff2d::defaultPhysicalParams; + auto icParamsVec = impladvdiff2d::defaultInitCondParams; + impladvdiff2d::replace_params_from_map_if_present(physParamsVec, icParamsVec, problemEnum, icFlag, userParams); + + using return_type = PublicProblemEigenMixinCpp>; + return return_type(meshObj, + problemEnum, + inviscidFluxRecEnum, + InviscidFluxScheme::Rusanov, + viscFluxRecEnum, + std::move(bcFuncs), + icParamsVec, + physParamsVec); +} +#endif + } //end namespace pressiodemoapps #endif diff --git a/include/pressiodemoapps/impl/advection_diffusion_2d_flux_mixin.hpp b/include/pressiodemoapps/impl/advection_diffusion_2d_flux_mixin.hpp index 1ad90b4c..7ce34e05 100644 --- a/include/pressiodemoapps/impl/advection_diffusion_2d_flux_mixin.hpp +++ b/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(m_fluxL, uMinusHalfNeg, - uMinusHalfPos, m_normal); - burgers_rusanov_flux_2d(m_fluxR, uPlusHalfNeg, - uPlusHalfPos, m_normal); - } + burgers_rusanov_flux_2d(m_fluxL, uMinusHalfNeg, + uMinusHalfPos, m_normal); + burgers_rusanov_flux_2d(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(m_fluxJacLNeg, m_fluxJacLPos, - uMinusHalfNeg, uMinusHalfPos, - m_normal); - burgers_rusanov_flux_jacobian_2d(m_fluxJacRNeg, m_fluxJacRPos, - uPlusHalfNeg, uPlusHalfPos, - m_normal); - } + burgers_rusanov_flux_jacobian_2d(m_fluxJacLNeg, m_fluxJacLPos, + uMinusHalfNeg, uMinusHalfPos, + m_normal); + burgers_rusanov_flux_jacobian_2d(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(m_fluxL, uMinusHalfNeg, uMinusHalfPos, m_normal); - burgers_rusanov_flux_2d(m_fluxR, uPlusHalfNeg, uPlusHalfPos, m_normal); - - burgers_rusanov_flux_jacobian_2d(m_fluxJacLNeg, m_fluxJacLPos, - uMinusHalfNeg, uMinusHalfPos, - m_normal); - burgers_rusanov_flux_jacobian_2d(m_fluxJacRNeg, m_fluxJacRPos, - uPlusHalfNeg, uPlusHalfPos, - m_normal); - } + burgers_rusanov_flux_2d(m_fluxL, uMinusHalfNeg, uMinusHalfPos, m_normal); + burgers_rusanov_flux_2d(m_fluxR, uPlusHalfNeg, uPlusHalfPos, m_normal); + + burgers_rusanov_flux_jacobian_2d(m_fluxJacLNeg, m_fluxJacLPos, + uMinusHalfNeg, uMinusHalfPos, + m_normal); + burgers_rusanov_flux_jacobian_2d(m_fluxJacRNeg, m_fluxJacRPos, + uPlusHalfNeg, uPlusHalfPos, + m_normal); break; } } diff --git a/include/pressiodemoapps/impl/advection_diffusion_2d_ghost_filler_outflow.hpp b/include/pressiodemoapps/impl/advection_diffusion_2d_ghost_filler_outflow.hpp new file mode 100644 index 00000000..5d0d0c53 --- /dev/null +++ b/include/pressiodemoapps/impl/advection_diffusion_2d_ghost_filler_outflow.hpp @@ -0,0 +1,246 @@ +/* +//@HEADER +// ************************************************************************ +// +// advection_diffusion_2d_ghost_filler_outflow.hpp +// Pressio +// Copyright 2019 +// National Technology & Engineering Solutions of Sandia, LLC (NTESS) +// +// Under the terms of Contract DE-NA0003525 with NTESS, the +// U.S. Government retains certain rights in this software. +// +// Pressio is licensed under BSD-3-Clause terms of use: +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions +// are met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +// COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING +// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Francesco Rizzi (fnrizzi@sandia.gov) +// +// ************************************************************************ +//@HEADER +*/ + +#ifndef PRESSIODEMOAPPS_GHOST_FILLER_ADVECTION_DIFFUSION_OUTFLOW_HPP_ +#define PRESSIODEMOAPPS_GHOST_FILLER_ADVECTION_DIFFUSION_OUTFLOW_HPP_ + +namespace pressiodemoapps{ +namespace impladvdiff2d{ + +template +class Ghost2dOutflowFiller +{ + +public: + Ghost2dOutflowFiller() = delete; + Ghost2dOutflowFiller(const int stencilSize, + int numDofPerCell, + const state_t & stateIn, + const mesh_t & meshIn, + ghost_t & ghostLeft, + ghost_t & ghostFront, + ghost_t & ghostRight, + ghost_t & ghostBack) + : m_stencilSize(stencilSize), + m_numDofPerCell(numDofPerCell), + m_state(stateIn), + m_meshObj(meshIn), + m_ghostLeft(ghostLeft), + m_ghostFront(ghostFront), + m_ghostRight(ghostRight), + m_ghostBack(ghostBack) + {} + + template + void operator()(index_t smPt, int gRow) + { + if (m_stencilSize == 3){ + stencilThreeImpl(smPt, gRow); + } + + else if (m_stencilSize == 5){ + stencilFiveImpl(smPt, gRow); + } + + else if (m_stencilSize == 7){ + stencilSevenImpl(smPt, gRow); + } + + else{ + throw std::runtime_error("advection diffusion neumann ghost filler: invalid stencil size"); + } + } + +private: + + template + void stencilThreeImpl(index_t smPt, int gRow) + { + const auto & graph = m_meshObj.graph(); + const auto cellGID = graph(smPt, 0); + const auto uIndex = cellGID*m_numDofPerCell; + + const auto left0 = graph(smPt, 1); + const auto front0 = graph(smPt, 2); + const auto right0 = graph(smPt, 3); + const auto back0 = graph(smPt, 4); + + if (left0 == -1) + { + m_ghostLeft(gRow, 0) = 0.0; + m_ghostLeft(gRow, 1) = 0.0; + } + + if (front0 == -1) + { + m_ghostFront(gRow, 0) = m_state(uIndex); + m_ghostFront(gRow, 1) = m_state(uIndex+1); + } + + if (right0 == -1) + { + m_ghostRight(gRow, 0) = m_state(uIndex); + m_ghostRight(gRow, 1) = m_state(uIndex+1); + } + + if (back0 == -1) + { + m_ghostBack(gRow, 0) = 0.0; + m_ghostBack(gRow, 1) = 0.0; + } + } + + template + void stencilFiveImpl(index_t smPt, int gRow) + { + const auto & graph = m_meshObj.graph(); + const auto cellGID = graph(smPt, 0); + const auto uIndex = cellGID*m_numDofPerCell; + + stencilThreeImpl(smPt, gRow); + const auto left0 = graph(smPt, 1); + const auto front0 = graph(smPt, 2); + const auto right0 = graph(smPt, 3); + const auto back0 = graph(smPt, 4); + const auto left1 = graph(smPt, 5); + const auto front1 = graph(smPt, 6); + const auto right1 = graph(smPt, 7); + const auto back1 = graph(smPt, 8); + + if (left1 == -1){ + m_ghostLeft(gRow, 2) = 0.0; + m_ghostLeft(gRow, 3) = 0.0; + } + + if (front1 == -1){ + auto ind = uIndex; + if (front0==-1){ ind = back0*m_numDofPerCell; } + else { ind = front0*m_numDofPerCell; } + + m_ghostFront(gRow, 2) = m_state(ind); + m_ghostFront(gRow, 3) = m_state(ind+1); + } + + if (right1 == -1){ + auto ind = uIndex; + if (right0==-1){ ind = left0*m_numDofPerCell; } + else { ind = right0*m_numDofPerCell; } + + m_ghostRight(gRow, 2) = m_state(ind); + m_ghostRight(gRow, 3) = m_state(ind+1); + } + + if (back1 == -1){ + m_ghostBack(gRow, 2) = 0.0; + m_ghostBack(gRow, 3) = 0.0; + } + } + + template + void stencilSevenImpl(index_t smPt, int gRow) + { + const auto & graph = m_meshObj.graph(); + const auto cellGID = graph(smPt, 0); + const auto uIndex = cellGID*m_numDofPerCell; + + stencilFiveImpl(smPt, gRow); + const auto front0 = graph(smPt, 2); + const auto right0 = graph(smPt, 3); + const auto left1 = graph(smPt, 5); + const auto front1 = graph(smPt, 6); + const auto right1 = graph(smPt, 7); + const auto back1 = graph(smPt, 8); + const auto left2 = graph(smPt, 9); + const auto front2 = graph(smPt, 10); + const auto right2 = graph(smPt, 11); + const auto back2 = graph(smPt, 12); + + if (left2 == -1){ + m_ghostLeft(gRow, 4) = 0.0; + m_ghostLeft(gRow, 5) = 0.0; + } + + if (front2 == -1){ + auto ind = uIndex; ; + if (front1!=-1 && front0!=-1){ ind = front1*m_numDofPerCell; } + if (front1==-1 && front0!=-1){ ind = uIndex; } + if (front1==-1 && front0==-1){ ind = back1*m_numDofPerCell; } + + m_ghostFront(gRow, 4) = m_state(ind); + m_ghostFront(gRow, 5) = m_state(ind+1); + } + + if (right2 == -1){ + auto ind = uIndex; ; + if (right1!=-1 && right0!=-1){ ind = right1*m_numDofPerCell; } + if (right1==-1 && right0!=-1){ ind = uIndex; } + if (right1==-1 && right0==-1){ ind = left1*m_numDofPerCell; } + + m_ghostRight(gRow, 4) = m_state(ind); + m_ghostRight(gRow, 5) = m_state(ind+1); + } + + if (back2 == -1){ + m_ghostBack(gRow, 4) = 0.0; + m_ghostBack(gRow, 5) = 0.0; + } + } + +private: + int m_stencilSize; + int m_numDofPerCell; + const state_t & m_state; + const mesh_t & m_meshObj; + ghost_t & m_ghostLeft; + ghost_t & m_ghostFront; + ghost_t & m_ghostRight; + ghost_t & m_ghostBack; +}; + +}} +#endif \ No newline at end of file diff --git a/include/pressiodemoapps/impl/advection_diffusion_2d_parametrization_helpers.hpp b/include/pressiodemoapps/impl/advection_diffusion_2d_parametrization_helpers.hpp new file mode 100644 index 00000000..b7132c52 --- /dev/null +++ b/include/pressiodemoapps/impl/advection_diffusion_2d_parametrization_helpers.hpp @@ -0,0 +1,141 @@ +/* +//@HEADER +// ************************************************************************ +// +// advection_diffusion_2d_parametrization_helpers.hpp +// Pressio +// Copyright 2019 +// National Technology & Engineering Solutions of Sandia, LLC (NTESS) +// +// Under the terms of Contract DE-NA0003525 with NTESS, the +// U.S. Government retains certain rights in this software. +// +// Pressio is licensed under BSD-3-Clause terms of use: +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions +// are met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +// COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING +// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Francesco Rizzi (fnrizzi@sandia.gov) +// +// ************************************************************************ +//@HEADER +*/ + +#ifndef ADVECTION_DIFFUSION_2D_PARAMETRIZATION_HELPERS_HPP_ +#define ADVECTION_DIFFUSION_2D_PARAMETRIZATION_HELPERS_HPP_ + +namespace pressiodemoapps{ +namespace impladvdiff2d{ + +constexpr int invalidIndex = std::numeric_limits::max(); + +// NOTE: physical and IC params must be associated with +// contiguous but separate indexes because they are stored in separate vectors +constexpr int diffusion_i = 0; + +constexpr int ic1_pulseMag_i = 0; +constexpr int ic1_pulseSpread_i = 1; +constexpr int ic1_pulseX_i = 2; +constexpr int ic1_pulseY_i = 3; + +// IMPORTANT: we need to fill the vectors with default vaues +// with ALL default paramters **consistently** with the indexing defined above +template +const std::vector defaultPhysicalParams +({ + static_cast(0.00001), // diffusion +}); + +template +const std::vector defaultInitCondParams +({ + static_cast(0.5), // pulseMag + static_cast(0.15), // pulseSpread + static_cast(0.0), // pulseX + static_cast(-0.2) // pulseY +}); + +const std::vector paramNames({ + "diffusion", + "pulseMagnitude", "pulseSpread", + "pulseX", "pulseY" + }); + +template +bool is_physical(const std::string & s){ + return (s == "diffusion"); +} + +template +int phys_param_string_to_index(const std::string & s){ + if (s == "diffusion") { return diffusion_i; } + else { return invalidIndex; } +} + +template +int ic_param_string_to_index(const int icFlag, const std::string & s){ + switch(icFlag){ + case 1: + if (s == "pulseMagnitude") { return ic1_pulseMag_i; } + else if (s == "pulseSpread") { return ic1_pulseSpread_i; } + else if (s == "pulseX") { return ic1_pulseX_i; } + else if (s == "pulseY") { return ic1_pulseY_i; } + else { return invalidIndex; } + + default: return invalidIndex; + } +} + +template +bool valid_parameter_name(const std::string & pname){ + auto match = [&](const std::string & testS){ return pname == testS; }; + return std::any_of(paramNames.cbegin(), paramNames.cend(), match); +} + +template +bool contains_valid_parameter_names(const std::unordered_map & map){ + return std::all_of(map.cbegin(), map.cend(), [](const auto & it){ return valid_parameter_name(it.first); }); +} + +template +void replace_params_from_map_if_present(std::vector & physVec, + std::vector & icVec, + const AdvectionDiffusion2d probEn, + const int icFlag, + const std::unordered_map & map) +{ + + auto action = [&](auto it){ + if (is_physical(it.first)){ physVec[phys_param_string_to_index<>(it.first)] = it.second; } + else { icVec[ic_param_string_to_index<>(icFlag, it.first)] = it.second;} + }; + std::for_each(map.cbegin(), map.cend(), action); +} + +}}//end namespace +#endif diff --git a/include/pressiodemoapps/impl/advection_diffusion_2d_prob_class.hpp b/include/pressiodemoapps/impl/advection_diffusion_2d_prob_class.hpp index ea7c5dc1..39580320 100644 --- a/include/pressiodemoapps/impl/advection_diffusion_2d_prob_class.hpp +++ b/include/pressiodemoapps/impl/advection_diffusion_2d_prob_class.hpp @@ -51,6 +51,8 @@ #include "advection_diffusion_2d_flux_functions.hpp" #include "advection_diffusion_2d_initial_condition.hpp" +#include "advection_diffusion_2d_ghost_filler_outflow.hpp" +#include "advection_diffusion_2d_parametrization_helpers.hpp" #include "functor_fill_stencil.hpp" #include "functor_reconstruct_from_stencil.hpp" #include "functor_reconstruct_from_state.hpp" @@ -58,6 +60,9 @@ #include "mixin_directional_flux_balance.hpp" #include "mixin_directional_flux_balance_jacobian.hpp" #include "Eigen/Sparse" +#include "noop.hpp" +#include "custom_bcs_functions.hpp" +#include "ghost_relative_locations.hpp" #ifdef PRESSIODEMOAPPS_ENABLE_OPENMP #include @@ -65,6 +70,16 @@ namespace pressiodemoapps{ namespace impladvdiff2d{ +template +bool valid_ic_flag(const AdvectionDiffusion2d probEn, const int flag) +{ + switch(probEn){ + case AdvectionDiffusion2d::BurgersPeriodic: return flag == 1; + case AdvectionDiffusion2d::BurgersOutflow: return flag == 1; + default: return false; + } +} + // tags are used inside he public create function: create_problem_...() // in the file ../advection_diffusion.hpp // to dispatch to the proper problem @@ -75,7 +90,10 @@ struct TagBurgersPeriodic{}; ///////////////////////// // eigen class ///////////////////////// -template +template< + class MeshType, + class BCFunctorsHolderType = impl::NoOperation + > class EigenApp { @@ -88,84 +106,116 @@ class EigenApp using jacobian_type = Eigen::SparseMatrix; using mesh_connectivity_graph_type = typename MeshType::graph_t; -private: static constexpr int dimensionality{2}; + static constexpr int numDofPerCell{2}; + +private: using ghost_container_type = Eigen::Matrix; using stencil_container_type = Eigen::Matrix; - using flux_type = Eigen::Matrix; - using edge_rec_type = Eigen::Matrix; - using flux_jac_type = Eigen::Matrix; + using flux_type = Eigen::Matrix; + using edge_rec_type = Eigen::Matrix; + using flux_jac_type = Eigen::Matrix; using reconstruction_gradient_t = Eigen::Matrix; public: EigenApp() = delete; // - // constructor for Burgers2d periodic + // constructor for Burgers2d // - EigenApp(TagBurgersPeriodic /*tag*/, - const MeshType & meshObj, + EigenApp(const MeshType & meshObj, + ::pressiodemoapps::AdvectionDiffusion2d probEnum, ::pressiodemoapps::InviscidFluxReconstruction inviscidFluxRecEn, ::pressiodemoapps::InviscidFluxScheme invFluxSchemeEn, ::pressiodemoapps::ViscousFluxReconstruction visFluxRecEn, - scalar_type icPulseMagnitude, - scalar_type icSpread, - scalar_type diffusionCoeff, - scalar_type x0, - scalar_type y0) - : m_numDofPerCell(2), - m_probEn(::pressiodemoapps::AdvectionDiffusion2d::BurgersPeriodic), + const std::vector & icParameters, + const std::vector & physParameters) + : m_probEn(probEnum), m_inviscidFluxRecEn(inviscidFluxRecEn), m_inviscidFluxSchemeEn(invFluxSchemeEn), m_viscousFluxRecEn(visFluxRecEn), m_meshObj(meshObj), - m_burgers2d_icPulse(icPulseMagnitude), - m_burgers2d_icSpread(icSpread), - m_burgers2d_diffusion(diffusionCoeff), - m_burgers2d_x0(x0), - m_burgers2d_y0(y0) + m_ic_parameters(icParameters), + m_phys_parameters(physParameters) { - m_numDofStencilMesh = m_meshObj.get().stencilMeshSize() * m_numDofPerCell; - m_numDofSampleMesh = m_meshObj.get().sampleMeshSize() * m_numDofPerCell; - // don't need to allocate ghosts because it is periodic + m_numDofStencilMesh = m_meshObj.get().stencilMeshSize() * numDofPerCell; + m_numDofSampleMesh = m_meshObj.get().sampleMeshSize() * numDofPerCell; + + if (m_probEn == pressiodemoapps::AdvectionDiffusion2d::BurgersOutflow) { + allocateGhosts(); + } + } - state_type initialCondition() const +#if !defined PRESSIODEMOAPPS_ENABLE_BINDINGS + EigenApp(const MeshType & meshObj, + ::pressiodemoapps::AdvectionDiffusion2d probEnum, + ::pressiodemoapps::InviscidFluxReconstruction inviscidFluxRecEn, + ::pressiodemoapps::InviscidFluxScheme invFluxSchemeEn, + ::pressiodemoapps::ViscousFluxReconstruction visFluxRecEn, + BCFunctorsHolderType && bcHolder, + const std::vector & icParameters, + const std::vector & physParameters) + : m_probEn(probEnum), + m_inviscidFluxRecEn(inviscidFluxRecEn), + m_inviscidFluxSchemeEn(invFluxSchemeEn), + m_viscousFluxRecEn(visFluxRecEn), + m_meshObj(meshObj), + m_ic_parameters(icParameters), + m_phys_parameters(physParameters), + m_bcFuncsHolder(std::move(bcHolder)) { + m_numDofStencilMesh = m_meshObj.get().stencilMeshSize() * numDofPerCell; + m_numDofSampleMesh = m_meshObj.get().sampleMeshSize() * numDofPerCell; + + if (m_probEn == pressiodemoapps::AdvectionDiffusion2d::BurgersOutflow) { + allocateGhosts(); + } + + } +#endif + + state_type initialCondition() const{ state_type initialState(m_numDofStencilMesh); - if (m_probEn == ::pressiodemoapps::AdvectionDiffusion2d::BurgersPeriodic) + if (m_probEn == ::pressiodemoapps::AdvectionDiffusion2d::BurgersPeriodic || + m_probEn == ::pressiodemoapps::AdvectionDiffusion2d::BurgersOutflow) { burgers2d_gaussian(initialState, m_meshObj.get(), - m_burgers2d_icPulse, - m_burgers2d_icSpread, - m_burgers2d_x0, - m_burgers2d_y0); + m_ic_parameters[ic1_pulseMag_i], + m_ic_parameters[ic1_pulseSpread_i], + m_ic_parameters[ic1_pulseX_i], + m_ic_parameters[ic1_pulseY_i]); } return initialState; } +public: + template + void setBCPointer(::pressiodemoapps::impl::GhostRelativeLocation rloc, T* ptr) { + m_bcFuncsHolder.setInternalPointer(rloc, ptr); + } + protected: int numDofPerCellImpl() const { - return m_numDofPerCell; + return numDofPerCell; } void initializeJacobian(jacobian_type & J) { J.resize(m_numDofSampleMesh, m_numDofStencilMesh); - using Tr = Eigen::Triplet; std::vector trList; - // only need this since this is for now periodic + initializeJacobianForNearBoundaryCells(trList); initializeJacobianForInnerCells(trList); - J.setFromTriplets(trList.begin(), trList.end()); + J.setFromTriplets(trList.begin(), trList.end()); // compress to make it Csr if (!J.isCompressed()){ J.makeCompressed(); @@ -188,15 +238,11 @@ class EigenApp // for omp, these are private variables for each thread // edge reconstructions - edge_rec_type uMinusHalfNeg(m_numDofPerCell); - edge_rec_type uMinusHalfPos(m_numDofPerCell); - edge_rec_type uPlusHalfNeg(m_numDofPerCell); - edge_rec_type uPlusHalfPos(m_numDofPerCell); + edge_rec_type uMinusHalfNeg, uMinusHalfPos; + edge_rec_type uPlusHalfNeg, uPlusHalfPos; // fluxes - flux_type fluxL(m_numDofPerCell); - flux_type fluxF(m_numDofPerCell);; - flux_type fluxR(m_numDofPerCell); - flux_type fluxB(m_numDofPerCell); + flux_type fluxL, fluxF; + flux_type fluxR, fluxB; #ifdef PRESSIODEMOAPPS_ENABLE_OPENMP ::pressiodemoapps::set_zero_omp(V); @@ -212,13 +258,28 @@ class EigenApp #endif } + if constexpr(std::is_same_v>){ + fillGhosts(U, currentTime); + } + else { + fillGhostsUseCustomFunctors(U, currentTime, m_meshObj, m_bcFuncsHolder, + m_ghostLeft, m_ghostFront, + m_ghostRight, m_ghostBack, numDofPerCell); + } + if (J){ - velocityAndJacobianPeriodicImpl(U, currentTime, V, *J, - fluxL, fluxF, fluxR, fluxB, - uMinusHalfNeg, uMinusHalfPos, - uPlusHalfNeg, uPlusHalfPos); + velocityAndJacobianImpl(U, currentTime, V, *J, + fluxL, fluxF, fluxR, fluxB, + uMinusHalfNeg, uMinusHalfPos, + uPlusHalfNeg, uPlusHalfPos); } else{ + if (!m_meshObj.get().isFullyPeriodic()) { + velocityOnlyNearBdCellsImpl(U, currentTime, V, + fluxL, fluxF, fluxR, fluxB, + uMinusHalfNeg, uMinusHalfPos, + uPlusHalfNeg, uPlusHalfPos); + } velocityOnlyInnerCellsImpl(U, currentTime, V, fluxL, fluxF, fluxR, fluxB, @@ -233,6 +294,28 @@ class EigenApp } private: + template + void fillGhosts(const U_t & U, const scalar_type currentTime) const + { + const auto stencilSize = reconstructionTypeToStencilSize(m_inviscidFluxRecEn); + if (m_probEn == pressiodemoapps::AdvectionDiffusion2d::BurgersOutflow) + { + using ghost_filler_t = Ghost2dOutflowFiller< + U_t, MeshType, ghost_container_type>; + ghost_filler_t ghF(stencilSize, numDofPerCell, + U, m_meshObj.get(), + m_ghostLeft, m_ghostFront, + m_ghostRight, m_ghostBack); + const auto & rowsBd = m_meshObj.get().graphRowsOfCellsNearBd(); +#ifdef PRESSIODEMOAPPS_ENABLE_OPENMP +#pragma omp for schedule(static) +#endif + for (decltype(rowsBd.size()) it=0; it void initializeJacobianForInnerCells(std::vector & trList) { @@ -248,12 +331,12 @@ class EigenApp const auto smPt = targetGraphRows[it]; // find out which row in the jacobian we are dealing with - const auto jacRowOfCurrCellFirstDof = smPt*m_numDofPerCell; + const auto jacRowOfCurrCellFirstDof = smPt*numDofPerCell; // initialize jacobian block entries wrt current cell's dofs - const auto jacColOfCurrCellRho = graph(smPt, 0)*m_numDofPerCell; - for (int k=0; k 0); for (int i=1; i<=numNeighbors; ++i){ - const auto colInd = graph(smPt, i)*m_numDofPerCell; - for (int k=0; k + void initializeJacobianForNearBoundaryCells(std::vector & trList) + { + const scalar_type zero = 0; + const auto & graph = m_meshObj.get().graph(); + const auto & targetGraphRows = m_meshObj.get().graphRowsOfCellsNearBd(); + for (std::size_t it=0; it - void velocityAndJacobianPeriodicImpl(const U_t & U, - const scalar_type currentTime, - V_t & V, - jacobian_type & J, - flux_type & fluxL, - flux_type & fluxF, - flux_type & fluxR, - flux_type & fluxB, - edge_rec_type & uMinusHalfNeg, - edge_rec_type & uMinusHalfPos, - edge_rec_type & uPlusHalfNeg, - edge_rec_type & uPlusHalfPos) const + void velocityAndJacobianImpl(const U_t & U, + const scalar_type currentTime, + V_t & V, + jacobian_type & J, + flux_type & fluxL, + flux_type & fluxF, + flux_type & fluxR, + flux_type & fluxB, + edge_rec_type & uMinusHalfNeg, + edge_rec_type & uMinusHalfPos, + edge_rec_type & uPlusHalfNeg, + edge_rec_type & uPlusHalfPos) const { // flux jacobians - flux_jac_type fluxJacLNeg(m_numDofPerCell, m_numDofPerCell); - flux_jac_type fluxJacLPos(m_numDofPerCell, m_numDofPerCell); - flux_jac_type fluxJacFNeg(m_numDofPerCell, m_numDofPerCell); - flux_jac_type fluxJacFPos(m_numDofPerCell, m_numDofPerCell); - flux_jac_type fluxJacRNeg(m_numDofPerCell, m_numDofPerCell); - flux_jac_type fluxJacRPos(m_numDofPerCell, m_numDofPerCell); - flux_jac_type fluxJacBNeg(m_numDofPerCell, m_numDofPerCell); - flux_jac_type fluxJacBPos(m_numDofPerCell, m_numDofPerCell); + flux_jac_type fluxJacLNeg; + flux_jac_type fluxJacLPos; + flux_jac_type fluxJacFNeg; + flux_jac_type fluxJacFPos; + flux_jac_type fluxJacRNeg; + flux_jac_type fluxJacRPos; + flux_jac_type fluxJacBNeg; + flux_jac_type fluxJacBPos; int nonZerosCountBeforeComputing = J.nonZeros(); - // this is for periodic so we don't have boundary cells, - // all cells count as "inner cells" + if (!m_meshObj.get().isFullyPeriodic()){ + if (m_inviscidFluxRecEn == InviscidFluxReconstruction::FirstOrder){ + velocityAndJacNearBDCellsImplFirstOrder(U, currentTime, V, J, + fluxL, fluxF, fluxR, fluxB, + fluxJacLNeg, fluxJacLPos, + fluxJacFNeg, fluxJacFPos, + fluxJacRNeg, fluxJacRPos, + fluxJacBNeg, fluxJacBPos, + uMinusHalfNeg, uMinusHalfPos, + uPlusHalfNeg, uPlusHalfPos); + } + else{ + velocityAndJacNearBDCellsImplDifferentSchemeNotFirstOrderInviscid( + U, currentTime, V, J, + fluxL, fluxF, fluxR, fluxB, + fluxJacLNeg, fluxJacLPos, + fluxJacFNeg, fluxJacFPos, + fluxJacRNeg, fluxJacRPos, + fluxJacBNeg, fluxJacBPos, + uMinusHalfNeg, uMinusHalfPos, + uPlusHalfNeg, uPlusHalfPos); + } + } + velocityAndJacInnerCellsImpl(U, currentTime, V, J, fluxL, fluxF, fluxR, fluxB, fluxJacLNeg, fluxJacLPos, @@ -363,14 +504,14 @@ class EigenApp const auto stencilSize = reconstructionTypeToStencilSize(m_inviscidFluxRecEn); assert(stencilSize >= reconstructionTypeToStencilSize(m_viscousFluxRecEn)); - reconstruction_gradient_t gradLNeg(m_numDofPerCell, stencilSize-1); - reconstruction_gradient_t gradLPos(m_numDofPerCell, stencilSize-1); - reconstruction_gradient_t gradFNeg(m_numDofPerCell, stencilSize-1); - reconstruction_gradient_t gradFPos(m_numDofPerCell, stencilSize-1); - reconstruction_gradient_t gradRNeg(m_numDofPerCell, stencilSize-1); - reconstruction_gradient_t gradRPos(m_numDofPerCell, stencilSize-1); - reconstruction_gradient_t gradBNeg(m_numDofPerCell, stencilSize-1); - reconstruction_gradient_t gradBPos(m_numDofPerCell, stencilSize-1); + reconstruction_gradient_t gradLNeg(numDofPerCell, stencilSize-1); + reconstruction_gradient_t gradLPos(numDofPerCell, stencilSize-1); + reconstruction_gradient_t gradFNeg(numDofPerCell, stencilSize-1); + reconstruction_gradient_t gradFPos(numDofPerCell, stencilSize-1); + reconstruction_gradient_t gradRNeg(numDofPerCell, stencilSize-1); + reconstruction_gradient_t gradRPos(numDofPerCell, stencilSize-1); + reconstruction_gradient_t gradBNeg(numDofPerCell, stencilSize-1); + reconstruction_gradient_t gradBPos(numDofPerCell, stencilSize-1); using functor_type = pda::impl::ComputeDirectionalFluxBalance< @@ -416,15 +557,515 @@ class EigenApp for (decltype(graphRows.size()) it=0; it + void velocityAndJacNearBDCellsImplDifferentSchemeNotFirstOrderInviscid( + const state_t & state, + const scalar_type currentTime, + V_t & V, + jacobian_type & J, + flux_type & fluxL, + flux_type & fluxF, + flux_type & fluxR, + flux_type & fluxB, + flux_jac_type & fluxJacLNeg, + flux_jac_type & fluxJacLPos, + flux_jac_type & fluxJacFNeg, + flux_jac_type & fluxJacFPos, + flux_jac_type & fluxJacRNeg, + flux_jac_type & fluxJacRPos, + flux_jac_type & fluxJacBNeg, + flux_jac_type & fluxJacBPos, + edge_rec_type & uMinusHalfNeg, + edge_rec_type & uMinusHalfPos, + edge_rec_type & uPlusHalfNeg, + edge_rec_type & uPlusHalfPos) const + { + namespace pda = ::pressiodemoapps; + constexpr int xAxis = 1; + constexpr int yAxis = 2; + + // if here, then the velocity must be computed with Weno, + /// while the jacobian must be computed with first order + + using stencil_filler_t = pda::impl::StencilFiller; + + // ***************************** + // *** functors for velocity *** + // ***************************** + const auto stencilSizeForV = reconstructionTypeToStencilSize(m_inviscidFluxRecEn); + // this method should only be called for stencilsize 5 or 7, because + // the case =3 is handled specifically in another method + assert(stencilSizeForV == 5 or stencilSizeForV == 7); + + stencil_container_type stencilValsForV(numDofPerCell*stencilSizeForV); + + stencil_filler_t FillStencilVeloX(reconstructionTypeToStencilSize(m_inviscidFluxRecEn), + state, m_meshObj.get(), m_ghostLeft, m_ghostRight, + stencilValsForV, xAxis); + stencil_filler_t FillStencilVeloY(reconstructionTypeToStencilSize(m_inviscidFluxRecEn), + state, m_meshObj.get(), m_ghostBack, m_ghostFront, + stencilValsForV, yAxis); + + using velo_functor_type = + pda::impl::ComputeDirectionalFluxBalance< + pda::impladvdiff2d::ComputeDirectionalFluxValues< + pda::impl::ReconstructorFromStencil< + edge_rec_type, stencil_container_type>, + scalar_type, flux_type>, + V_t, scalar_type + >; + + velo_functor_type funcVeloX(V, m_meshObj.get().dxInv(), + /* end args for velo */ + m_probEn, m_inviscidFluxSchemeEn, normalX_, fluxL, fluxR, + /* end args for flux */ + toReconstructionScheme(m_inviscidFluxRecEn), stencilValsForV, + uMinusHalfNeg, uMinusHalfPos, uPlusHalfNeg, uPlusHalfPos + /* end args for reconstructor */ + ); + + velo_functor_type funcVeloY(V, m_meshObj.get().dyInv(), + /* end args for velo */ + m_probEn, m_inviscidFluxSchemeEn, normalY_, fluxB, fluxF, + /* end args for flux */ + toReconstructionScheme(m_inviscidFluxRecEn), stencilValsForV, + uMinusHalfNeg, uMinusHalfPos, uPlusHalfNeg, uPlusHalfPos + /* end args for reconstructor */ + ); + + // ***************************** + // *** functors for jacobian *** + // ***************************** + const auto firstOrderRec = pda::InviscidFluxReconstruction::FirstOrder; + const auto stencilSizeForJ = reconstructionTypeToStencilSize(firstOrderRec); + stencil_container_type stencilValsForJ(numDofPerCell*stencilSizeForJ); + stencil_filler_t FillStencilJacX(stencilSizeForJ, + state, m_meshObj.get(), m_ghostLeft, m_ghostRight, + stencilValsForJ, xAxis); + stencil_filler_t FillStencilJacY(stencilSizeForJ, + state, m_meshObj.get(), m_ghostBack, m_ghostFront, + stencilValsForJ, yAxis); + + using jac_functor_type = + pda::impl::ComputeDirectionalFluxBalanceFirstOrderJacobianOnBoundaryCell< + pda::impladvdiff2d::ComputeDirectionalFluxJacobians< + pda::impl::ReconstructorFromStencil< + edge_rec_type, stencil_container_type>, + scalar_type, flux_jac_type>, + dimensionality, MeshType, jacobian_type + >; + + jac_functor_type funcJacX(J, xAxis, m_meshObj.get(), + /* end args for jac */ + m_probEn, m_inviscidFluxSchemeEn, normalX_, + fluxJacLNeg, fluxJacLPos, fluxJacRNeg, fluxJacRPos, + /* end args for flux */ + toReconstructionScheme(firstOrderRec), stencilValsForJ, + uMinusHalfNeg, uMinusHalfPos, uPlusHalfNeg, uPlusHalfPos + /* end args for reconstructor */ + ); + + jac_functor_type funcJacY(J, yAxis, m_meshObj.get(), + /* end args for jac */ + m_probEn, m_inviscidFluxSchemeEn, normalY_, + fluxJacBNeg, fluxJacBPos, fluxJacFNeg, fluxJacFPos, + /* end args for flux */ + toReconstructionScheme(firstOrderRec), stencilValsForJ, + uMinusHalfNeg, uMinusHalfPos, uPlusHalfNeg, uPlusHalfPos + /* end args for reconstructor */ + ); + + // ************ + // loop + // ************ + const auto & graph = m_meshObj.get().graph(); + const auto & graphRows = m_meshObj.get().graphRowsOfCellsNearBd(); + const auto dxInvSq = m_meshObj.get().dxInv()*m_meshObj.get().dxInv(); + const auto dyInvSq = m_meshObj.get().dyInv()*m_meshObj.get().dyInv(); + const auto diffDxInvSq = m_phys_parameters[diffusion_i]*dxInvSq; + const auto diffDyInvSq = m_phys_parameters[diffusion_i]*dyInvSq; + constexpr auto two = static_cast(2); + +#ifdef PRESSIODEMOAPPS_ENABLE_OPENMP +#pragma omp for schedule(static) +#endif + for (decltype(graphRows.size()) it=0; it>){ + fillJacFactorsForCellBd(smPt, xAxis); + } + else{ + fillJacFactorsCustomBCs(smPt, xAxis, m_meshObj, m_bcFuncsHolder, + m_bcCellJacFactors, numDofPerCell); + } + funcJacX(smPt, numDofPerCell, m_bcCellJacFactors); + // diffusion contribution to velocity + if (stencilSizeForV == 5){ + V(vIndex) += diffDxInvSq*( stencilValsForV(6) - two*stencilValsForV(4) + stencilValsForV(2) ); + V(vIndex+1) += diffDxInvSq*( stencilValsForV(7) - two*stencilValsForV(5) + stencilValsForV(3) ); + } + else if (stencilSizeForV == 7){ + V(vIndex) += diffDxInvSq*( stencilValsForV(8) - two*stencilValsForV(6) + stencilValsForV(4) ); + V(vIndex+1) += diffDxInvSq*( stencilValsForV(9) - two*stencilValsForV(7) + stencilValsForV(5) ); + } + + // y-direction + FillStencilVeloY(smPt, it, numDofPerCell); + funcVeloY(smPt, numDofPerCell); + FillStencilJacY(smPt, it, numDofPerCell); + if constexpr(std::is_same_v>){ + fillJacFactorsForCellBd(smPt, yAxis); + } + else{ + fillJacFactorsCustomBCs(smPt, yAxis, m_meshObj, m_bcFuncsHolder, + m_bcCellJacFactors, numDofPerCell); + } + funcJacY(smPt, numDofPerCell, m_bcCellJacFactors); + // diffusion contribution to velocity + if (stencilSizeForV == 5){ + V(vIndex) += diffDyInvSq*( stencilValsForV(6) - two*stencilValsForV(4) + stencilValsForV(2) ); + V(vIndex+1) += diffDyInvSq*( stencilValsForV(7) - two*stencilValsForV(5) + stencilValsForV(3) ); + } + else if (stencilSizeForV == 7){ + V(vIndex) += diffDyInvSq*( stencilValsForV(8) - two*stencilValsForV(6) + stencilValsForV(4) ); + V(vIndex+1) += diffDyInvSq*( stencilValsForV(9) - two*stencilValsForV(7) + stencilValsForV(5) ); + } - if (m_probEn == ::pressiodemoapps::AdvectionDiffusion2d::BurgersPeriodic){ - addBurgersDiffusionAndSourceToVelocityAndJacobianInnerCells(U, V, J, smPt); + // diffusion contribution to Jacobian + J.coeffRef(vIndex, uIndex) += -two*diffDxInvSq - two*diffDyInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -two*diffDxInvSq - two*diffDyInvSq; + + if (uIndexLeft > -1){ + J.coeffRef(vIndex, uIndexLeft) += diffDxInvSq; + J.coeffRef(vIndex+1, uIndexLeft+1) += diffDxInvSq; + } + else{ + J.coeffRef(vIndex, uIndex) += -diffDxInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -diffDxInvSq; + } + + if (uIndexFront > -1){ + J.coeffRef(vIndex, uIndexFront) += diffDyInvSq; + J.coeffRef(vIndex+1, uIndexFront+1) += diffDyInvSq; + } + else{ + J.coeffRef(vIndex, uIndex) += -diffDyInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -diffDyInvSq; + } + + if (uIndexRight > -1){ + J.coeffRef(vIndex, uIndexRight) += diffDxInvSq; + J.coeffRef(vIndex+1, uIndexRight+1) += diffDxInvSq; + } + else{ + J.coeffRef(vIndex, uIndex) += -diffDxInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -diffDxInvSq; + } + + if (uIndexBack > -1){ + J.coeffRef(vIndex, uIndexBack) += diffDyInvSq; + J.coeffRef(vIndex+1, uIndexBack+1) += diffDyInvSq; + } + else{ + J.coeffRef(vIndex, uIndex) += -diffDyInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -diffDyInvSq; } } } + template + void velocityAndJacNearBDCellsImplFirstOrder(const state_t & state, + const scalar_type currentTime, + V_t & V, + jacobian_type & J, + flux_type & fluxL, + flux_type & fluxF, + flux_type & fluxR, + flux_type & fluxB, + flux_jac_type & fluxJacLNeg, + flux_jac_type & fluxJacLPos, + flux_jac_type & fluxJacFNeg, + flux_jac_type & fluxJacFPos, + flux_jac_type & fluxJacRNeg, + flux_jac_type & fluxJacRPos, + flux_jac_type & fluxJacBNeg, + flux_jac_type & fluxJacBPos, + edge_rec_type & uMinusHalfNeg, + edge_rec_type & uMinusHalfPos, + edge_rec_type & uPlusHalfNeg, + edge_rec_type & uPlusHalfPos) const + { + namespace pda = ::pressiodemoapps; + constexpr int xAxis = 1; + constexpr int yAxis = 2; + assert(m_inviscidFluxRecEn == InviscidFluxReconstruction::FirstOrder); + + // if here, then the scheme for velocity matches + // the one for Jacobian so we can use same functors + + const auto stencilSize = reconstructionTypeToStencilSize(m_inviscidFluxRecEn); + stencil_container_type stencilVals(numDofPerCell*stencilSize); + + using stencil_filler_t = pda::impl::StencilFiller< + dimensionality, stencil_container_type, state_t, MeshType, ghost_container_type>; + stencil_filler_t FillStencilX(reconstructionTypeToStencilSize(m_inviscidFluxRecEn), + state, m_meshObj.get(), m_ghostLeft, m_ghostRight, + stencilVals, xAxis); + + stencil_filler_t FillStencilY(reconstructionTypeToStencilSize(m_inviscidFluxRecEn), + state, m_meshObj.get(), m_ghostBack, m_ghostFront, + stencilVals, yAxis); + + using functor_type = + pda::impl::ComputeDirectionalFluxBalance< + pda::impl::ComputeDirectionalFluxBalanceFirstOrderJacobianOnBoundaryCell< + pda::impladvdiff2d::ComputeDirectionalFluxValuesAndJacobians< + pda::impl::ReconstructorFromStencil< + edge_rec_type, stencil_container_type>, + scalar_type, flux_type, flux_jac_type>, + dimensionality, MeshType, jacobian_type>, + V_t, scalar_type + >; + + functor_type funcx(V, m_meshObj.get().dxInv(), + /* end args for velo */ + J, xAxis, m_meshObj.get(), + /* end args for jac */ + m_probEn, m_inviscidFluxSchemeEn, normalX_, fluxL, fluxR, + fluxJacLNeg, fluxJacLPos, fluxJacRNeg, fluxJacRPos, + /* end args for flux */ + toReconstructionScheme(m_inviscidFluxRecEn), stencilVals, + uMinusHalfNeg, uMinusHalfPos, uPlusHalfNeg, uPlusHalfPos + /* end args for reconstructor */ + ); + + functor_type funcy(V, m_meshObj.get().dyInv(), + /* end args for velo */ + J, yAxis, m_meshObj.get(), + /* end args for jac */ + m_probEn, m_inviscidFluxSchemeEn, normalY_, fluxB, fluxF, + fluxJacBNeg, fluxJacBPos, fluxJacFNeg, fluxJacFPos, + /* end args for flux */ + toReconstructionScheme(m_inviscidFluxRecEn), stencilVals, + uMinusHalfNeg, uMinusHalfPos, uPlusHalfNeg, uPlusHalfPos + /* end args for reconstructor */ + ); + + const auto & graph = m_meshObj.get().graph(); + constexpr auto two = static_cast(2); + const auto & graphRows = m_meshObj.get().graphRowsOfCellsNearBd(); + const auto dxInvSq = m_meshObj.get().dxInv()*m_meshObj.get().dxInv(); + const auto dyInvSq = m_meshObj.get().dyInv()*m_meshObj.get().dyInv(); + const auto diffDxInvSq = m_phys_parameters[diffusion_i]*dxInvSq; + const auto diffDyInvSq = m_phys_parameters[diffusion_i]*dyInvSq; + +#ifdef PRESSIODEMOAPPS_ENABLE_OPENMP +#pragma omp for schedule(static) +#endif + for (decltype(graphRows.size()) it=0; it>){ + fillJacFactorsForCellBd(smPt, xAxis); + } + else { + fillJacFactorsCustomBCs(smPt, xAxis, m_meshObj, m_bcFuncsHolder, + m_bcCellJacFactors, numDofPerCell); + } + funcx(smPt, numDofPerCell, m_bcCellJacFactors); + // x-direction velocity contributions + V(vIndex) += diffDxInvSq*( stencilVals(4) - two*stencilVals(2) + stencilVals(0) ); + V(vIndex+1) += diffDxInvSq*( stencilVals(5) - two*stencilVals(3) + stencilVals(1) ); + + // y-direction inviscid contributions + FillStencilY(smPt, it, numDofPerCell); + if constexpr(std::is_same_v>){ + fillJacFactorsForCellBd(smPt, yAxis); + } + else { + fillJacFactorsCustomBCs(smPt, yAxis, m_meshObj, m_bcFuncsHolder, + m_bcCellJacFactors, numDofPerCell); + } + funcy(smPt, numDofPerCell, m_bcCellJacFactors); + // y-direction velocity contributions + V(vIndex) += diffDyInvSq*( stencilVals(4) - two*stencilVals(2) + stencilVals(0) ); + V(vIndex+1) += diffDyInvSq*( stencilVals(5) - two*stencilVals(3) + stencilVals(1) ); + + // diffusion Jacobian contributions + J.coeffRef(vIndex, uIndex) += -two*diffDxInvSq - two*diffDyInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -two*diffDxInvSq - two*diffDyInvSq; + + if (uIndexLeft > -1){ + J.coeffRef(vIndex, uIndexLeft) += diffDxInvSq; + J.coeffRef(vIndex+1, uIndexLeft+1) += diffDxInvSq; + } + else{ + J.coeffRef(vIndex, uIndex) += -diffDxInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -diffDxInvSq; + } + + if (uIndexFront > -1){ + J.coeffRef(vIndex, uIndexFront) += diffDyInvSq; + J.coeffRef(vIndex+1, uIndexFront+1) += diffDyInvSq; + } + else{ + J.coeffRef(vIndex, uIndex) += -diffDyInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -diffDyInvSq; + } + + if (uIndexRight > -1){ + J.coeffRef(vIndex, uIndexRight) += diffDxInvSq; + J.coeffRef(vIndex+1, uIndexRight+1) += diffDxInvSq; + } + else{ + J.coeffRef(vIndex, uIndex) += -diffDxInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -diffDxInvSq; + } + + if (uIndexBack > -1){ + J.coeffRef(vIndex, uIndexBack) += diffDyInvSq; + J.coeffRef(vIndex+1, uIndexBack+1) += diffDyInvSq; + } + else{ + J.coeffRef(vIndex, uIndex) += -diffDyInvSq; + J.coeffRef(vIndex+1, uIndex+1) += -diffDyInvSq; + } + } + } + + template + void velocityOnlyNearBdCellsImpl(const state_t & state, + const scalar_type currentTime, + V_t & V, + flux_type & fluxL, + flux_type & fluxF, + flux_type & fluxR, + flux_type & fluxB, + edge_rec_type & uMinusHalfNeg, + edge_rec_type & uMinusHalfPos, + edge_rec_type & uPlusHalfNeg, + edge_rec_type & uPlusHalfPos) const + { + namespace pda = ::pressiodemoapps; + constexpr int xAxis = 1; + constexpr int yAxis = 2; + + const auto stencilSize = reconstructionTypeToStencilSize(m_inviscidFluxRecEn); + stencil_container_type stencilVals(numDofPerCell*stencilSize); + + // stencil filler needed because we are doing cells near boundaries + using sfiller_t = ::pressiodemoapps::impl::StencilFiller< + dimensionality, stencil_container_type, state_t, MeshType, ghost_container_type>; + + sfiller_t StencilFillerX(stencilSize, state, m_meshObj.get(), m_ghostLeft, m_ghostRight, stencilVals, xAxis); + sfiller_t StencilFillerY(stencilSize, state, m_meshObj.get(), m_ghostBack, m_ghostFront, stencilVals, yAxis); + + using functor_type = + pda::impl::ComputeDirectionalFluxBalance< + pda::impladvdiff2d::ComputeDirectionalFluxValues< + pda::impl::ReconstructorFromStencil< + edge_rec_type, stencil_container_type>, + scalar_type, flux_type>, + V_t, scalar_type + >; + + functor_type Fx(V, m_meshObj.get().dxInv(), + /* end args for velo */ + m_probEn, m_inviscidFluxSchemeEn, normalX_, fluxL, fluxR, + /* end args for flux */ + toReconstructionScheme(m_inviscidFluxRecEn), stencilVals, + uMinusHalfNeg, uMinusHalfPos, uPlusHalfNeg, uPlusHalfPos + /* end args for reconstructor */ + ); + + functor_type Fy(V, m_meshObj.get().dyInv(), + /* end args for velo */ + m_probEn, m_inviscidFluxSchemeEn, normalY_, fluxB, fluxF, + /* end args for flux */ + toReconstructionScheme(m_inviscidFluxRecEn), stencilVals, + uMinusHalfNeg, uMinusHalfPos, uPlusHalfNeg, uPlusHalfPos + /* end args for reconstructor */ + ); + + constexpr auto two = static_cast(2); + const auto dxInvSq = m_meshObj.get().dxInv()*m_meshObj.get().dxInv(); + const auto dyInvSq = m_meshObj.get().dyInv()*m_meshObj.get().dyInv(); + const auto diffDxInvSq = m_phys_parameters[diffusion_i]*dxInvSq; + const auto diffDyInvSq = m_phys_parameters[diffusion_i]*dyInvSq; + + const auto & rows = m_meshObj.get().graphRowsOfCellsNearBd(); +#if defined PRESSIODEMOAPPS_ENABLE_OPENMP && !defined PRESSIODEMOAPPS_ENABLE_BINDINGS +#pragma omp for schedule(static) +#endif + for (std::size_t it=0; it void velocityOnlyInnerCellsImpl(const U_t & U, const scalar_type /*currentTime*/, @@ -475,86 +1116,116 @@ class EigenApp #endif for (decltype(graphRows.size()) it=0; it - void addBurgersDiffusionToVelocityInnerCells(const U_t & U, - V_t & V, - index_t smPt) const + void fillJacFactorsForCellBd(index_t graphRow, int axis) const { - constexpr auto two = static_cast(2); - const auto dxInvSq = m_meshObj.get().dxInv()*m_meshObj.get().dxInv(); - const auto dyInvSq = m_meshObj.get().dyInv()*m_meshObj.get().dyInv(); - const auto diffDxInvSq = m_burgers2d_diffusion*dxInvSq; - const auto diffDyInvSq = m_burgers2d_diffusion*dyInvSq; + assert(axis <= 2); - const auto & graph = m_meshObj.get().graph(); - const auto vIndex = smPt*m_numDofPerCell; - const auto uIndex = graph(smPt, 0)*m_numDofPerCell; - const auto uIndexLeft = graph(smPt, 1)*m_numDofPerCell; - const auto uIndexFront = graph(smPt, 2)*m_numDofPerCell; - const auto uIndexRight = graph(smPt, 3)*m_numDofPerCell; - const auto uIndexBack = graph(smPt, 4)*m_numDofPerCell; - V(vIndex) += diffDxInvSq*( U(uIndexRight) - two*U(uIndex) + U(uIndexLeft) ); - V(vIndex) += diffDyInvSq*( U(uIndexFront) - two*U(uIndex) + U(uIndexBack) ); + if (m_probEn == ::pressiodemoapps::AdvectionDiffusion2d::BurgersOutflow) { + if (axis == 1 && m_meshObj.get().hasBdLeft2d(graphRow)){ + // homogeneous dirichlet + m_bcCellJacFactors = {0., 0.}; + return; + } - V(vIndex+1) += diffDxInvSq*( U(uIndexRight+1) - two*U(uIndex+1) + U(uIndexLeft+1) ); - V(vIndex+1) += diffDyInvSq*( U(uIndexFront+1) - two*U(uIndex+1) + U(uIndexBack+1) ); + if (axis == 1 && m_meshObj.get().hasBdRight2d(graphRow)){ + // homogeneous neumann + m_bcCellJacFactors = {1., 1.}; + return; + } + + if (axis == 2 && m_meshObj.get().hasBdBack2d(graphRow)) + { + // homogeneous dirichlet + m_bcCellJacFactors = {0., 0.}; + return; + } + + if (axis == 2 && m_meshObj.get().hasBdFront2d(graphRow)) + { + // homogeneous neumann + m_bcCellJacFactors = {1., 1.}; + return; + } + } + else if (m_probEn == ::pressiodemoapps::AdvectionDiffusion2d::BurgersPeriodic) { + throw std::runtime_error("Should not be getting Jacobian factors for fully periodic problem."); + } + else { + throw std::runtime_error("Invalid problem enumeration."); + } } template - void addBurgersDiffusionAndSourceToVelocityAndJacobianInnerCells(const U_t & U, + void addBurgersDiffusionToVelocityAndOptionalJacobianInnerCells(const U_t & U, V_t & V, - jacobian_type & J, + jacobian_type * J, index_t smPt) const { constexpr auto two = static_cast(2); const auto dxInvSq = m_meshObj.get().dxInv()*m_meshObj.get().dxInv(); const auto dyInvSq = m_meshObj.get().dyInv()*m_meshObj.get().dyInv(); - const auto diffDxInvSq = m_burgers2d_diffusion*dxInvSq; - const auto diffDyInvSq = m_burgers2d_diffusion*dyInvSq; + const auto diffDxInvSq = m_phys_parameters[diffusion_i]*dxInvSq; + const auto diffDyInvSq = m_phys_parameters[diffusion_i]*dyInvSq; const auto & graph = m_meshObj.get().graph(); - const auto vIndex = smPt*m_numDofPerCell; - const auto uIndex = graph(smPt, 0)*m_numDofPerCell; - const auto uIndexLeft = graph(smPt, 1)*m_numDofPerCell; - const auto uIndexFront = graph(smPt, 2)*m_numDofPerCell; - const auto uIndexRight = graph(smPt, 3)*m_numDofPerCell; - const auto uIndexBack = graph(smPt, 4)*m_numDofPerCell; + const auto vIndex = smPt*numDofPerCell; + const auto uIndex = graph(smPt, 0)*numDofPerCell; + const auto uIndexLeft = graph(smPt, 1)*numDofPerCell; + const auto uIndexFront = graph(smPt, 2)*numDofPerCell; + const auto uIndexRight = graph(smPt, 3)*numDofPerCell; + const auto uIndexBack = graph(smPt, 4)*numDofPerCell; + V(vIndex) += diffDxInvSq*( U(uIndexRight) - two*U(uIndex) + U(uIndexLeft) ); V(vIndex) += diffDyInvSq*( U(uIndexFront) - two*U(uIndex) + U(uIndexBack) ); V(vIndex+1) += diffDxInvSq*( U(uIndexRight+1) - two*U(uIndex+1) + U(uIndexLeft+1) ); V(vIndex+1) += diffDyInvSq*( U(uIndexFront+1) - two*U(uIndex+1) + U(uIndexBack+1) ); - J.coeffRef(vIndex, uIndex) += -two*diffDxInvSq - two*diffDyInvSq; - J.coeffRef(vIndex, uIndexLeft) += diffDxInvSq; - J.coeffRef(vIndex, uIndexFront) += diffDyInvSq; - J.coeffRef(vIndex, uIndexRight) += diffDxInvSq; - J.coeffRef(vIndex, uIndexBack) += diffDyInvSq; - J.coeffRef(vIndex+1, uIndex+1) += -two*diffDxInvSq - two*diffDyInvSq; - J.coeffRef(vIndex+1, uIndexLeft+1) += diffDxInvSq; - J.coeffRef(vIndex+1, uIndexFront+1) += diffDyInvSq; - J.coeffRef(vIndex+1, uIndexRight+1) += diffDxInvSq; - J.coeffRef(vIndex+1, uIndexBack+1) += diffDyInvSq; + if (J){ + J->coeffRef(vIndex, uIndex) += -two*diffDxInvSq - two*diffDyInvSq; + J->coeffRef(vIndex, uIndexLeft) += diffDxInvSq; + J->coeffRef(vIndex, uIndexFront) += diffDyInvSq; + J->coeffRef(vIndex, uIndexRight) += diffDxInvSq; + J->coeffRef(vIndex, uIndexBack) += diffDyInvSq; + J->coeffRef(vIndex+1, uIndex+1) += -two*diffDxInvSq - two*diffDyInvSq; + J->coeffRef(vIndex+1, uIndexLeft+1) += diffDxInvSq; + J->coeffRef(vIndex+1, uIndexFront+1) += diffDyInvSq; + J->coeffRef(vIndex+1, uIndexRight+1) += diffDxInvSq; + J->coeffRef(vIndex+1, uIndexBack+1) += diffDyInvSq; + } + } + + void allocateGhosts() + { + const auto stencilSize = reconstructionTypeToStencilSize(m_inviscidFluxRecEn); + const auto numGhostValues = numDofPerCell*((stencilSize-1)/2); + + const index_t s1 = m_meshObj.get().numCellsNearBd(); + ::pressiodemoapps::resize(m_ghostLeft, s1, numGhostValues); + ::pressiodemoapps::resize(m_ghostFront,s1, numGhostValues); + ::pressiodemoapps::resize(m_ghostRight,s1, numGhostValues); + ::pressiodemoapps::resize(m_ghostBack, s1, numGhostValues); } protected: // common to all problems - int m_numDofPerCell = {}; ::pressiodemoapps::AdvectionDiffusion2d m_probEn; ::pressiodemoapps::InviscidFluxReconstruction m_inviscidFluxRecEn; ::pressiodemoapps::InviscidFluxScheme m_inviscidFluxSchemeEn; ::pressiodemoapps::ViscousFluxReconstruction m_viscousFluxRecEn; - std::reference_wrapper m_meshObj; + std::vector m_ic_parameters; + std::vector m_phys_parameters; + BCFunctorsHolderType m_bcFuncsHolder = {}; + index_t m_numDofStencilMesh = {}; index_t m_numDofSampleMesh = {}; @@ -566,16 +1237,11 @@ class EigenApp std::array normalX_{1, 0}; std::array normalY_{0, 1}; - // parameters specific to problems - // will need to handle this better later - scalar_type m_burgers2d_icPulse = {}; - scalar_type m_burgers2d_icSpread = {}; - scalar_type m_burgers2d_diffusion = {}; - scalar_type m_burgers2d_x0 = {}; - scalar_type m_burgers2d_y0 = {}; + mutable std::array m_bcCellJacFactors; }; -template constexpr int EigenApp::dimensionality; +template constexpr int EigenApp::numDofPerCell; +template constexpr int EigenApp::dimensionality; }}//end namespace #endif diff --git a/include/pressiodemoapps/impl/swe_2d_prob_class.hpp b/include/pressiodemoapps/impl/swe_2d_prob_class.hpp index 92566884..d21b4df5 100644 --- a/include/pressiodemoapps/impl/swe_2d_prob_class.hpp +++ b/include/pressiodemoapps/impl/swe_2d_prob_class.hpp @@ -97,10 +97,11 @@ class EigenApp using jacobian_type = Eigen::SparseMatrix; using mesh_connectivity_graph_type = typename MeshType::graph_t; -private: static constexpr int dimensionality{2}; static constexpr int numDofPerCell{3}; +private: + using ghost_container_type = Eigen::Matrix; using stencil_container_type = Eigen::Matrix; using flux_type = Eigen::Matrix; diff --git a/pressiodemoapps/__init__.py b/pressiodemoapps/__init__.py index fdf2f489..f4886d8e 100644 --- a/pressiodemoapps/__init__.py +++ b/pressiodemoapps/__init__.py @@ -75,7 +75,7 @@ import create_problem, \ create_linear_advection_1d_problem, \ create_diffusion_reaction_1d_problem_A,\ - create_periodic_burgers_2d_problem,\ + create_burgers_2d_problem,\ create_diffusion_reaction_2d_problem_A,\ create_adv_diff_reac_2d_problem_A, \ create_gray_scott_2d_problem,\ diff --git a/src_py/main_binder.cc b/src_py/main_binder.cc index 95111eb8..0c3eec9b 100644 --- a/src_py/main_binder.cc +++ b/src_py/main_binder.cc @@ -164,7 +164,9 @@ void bind2dProblemEnums(pybind11::module & mParent) pybind11::enum_(mParent, "AdvectionDiffusion2d") .value("BurgersPeriodic", - pda::AdvectionDiffusion2d::BurgersPeriodic); + pda::AdvectionDiffusion2d::BurgersPeriodic) + .value("BurgersOutflow", + pda::AdvectionDiffusion2d::BurgersOutflow); pybind11::enum_(mParent, "DiffusionReaction2d") .value("ProblemA", @@ -389,13 +391,12 @@ void bindAdvectionDiffusion2d(pybind11::module & mParent) pybind11::arg().noconvert(), pybind11::arg().noconvert(), pybind11::arg().noconvert()); - mParent.def("create_periodic_burgers_2d_problem", - &pda::create_periodic_burgers_2d_problem_ov1_for_py, + mParent.def("create_burgers_2d_problem", + &pda::create_burgers_2d_problem_ov1_for_py, pybind11::return_value_policy::take_ownership, pybind11::arg().noconvert(), pybind11::arg().noconvert(), pybind11::arg().noconvert(), pybind11::arg().noconvert(), - pybind11::arg().noconvert(), pybind11::arg().noconvert(), - pybind11::arg().noconvert(), pybind11::arg().noconvert()); + pybind11::arg().noconvert()); } // ----------------------- diff --git a/tests_cpp/CMakeLists.txt b/tests_cpp/CMakeLists.txt index 616fa26e..01e77141 100644 --- a/tests_cpp/CMakeLists.txt +++ b/tests_cpp/CMakeLists.txt @@ -87,6 +87,7 @@ add_subdirectory(eigen_2d_gray_scott_explicit) add_subdirectory(eigen_2d_gray_scott_sample_mesh_test) add_subdirectory(eigen_2d_burgers_periodic_explicit) add_subdirectory(eigen_2d_burgers_periodic_sample_mesh_test) +add_subdirectory(eigen_2d_burgers_outflow_implicit) add_subdirectory(eigen_2d_advdiffreac_probA_explicit) add_subdirectory(eigen_2d_advdiffreac_probA_sample_mesh_test) diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/CMakeLists.txt b/tests_cpp/eigen_2d_burgers_outflow_implicit/CMakeLists.txt new file mode 100644 index 00000000..3de92773 --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/CMakeLists.txt @@ -0,0 +1,4 @@ + +add_subdirectory(firstorder) +add_subdirectory(weno3) +add_subdirectory(weno5) diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/compare.py b/tests_cpp/eigen_2d_burgers_outflow_implicit/compare.py new file mode 100644 index 00000000..0f364f1e --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/compare.py @@ -0,0 +1,20 @@ + +import numpy as np +import sys, os + +if __name__== "__main__": + nx=20 + ny=20 + fomTotDofs = nx*ny*2 + + D = np.fromfile("burgers2d_solution.bin") + nt = int(np.size(D)/fomTotDofs) + D = np.reshape(D, (nt, fomTotDofs)) + D = D[-1, :] + np.savetxt("field.txt", D) + + goldD = np.loadtxt("gold.txt") + assert(np.allclose(D.shape, goldD.shape)) + assert(np.isnan(D).all() == False) + assert(np.isnan(goldD).all() == False) + assert(np.allclose(D, goldD,rtol=1e-10, atol=1e-12)) diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/firstorder/CMakeLists.txt b/tests_cpp/eigen_2d_burgers_outflow_implicit/firstorder/CMakeLists.txt new file mode 100644 index 00000000..2e7e1cbb --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/firstorder/CMakeLists.txt @@ -0,0 +1,18 @@ + +set(testname eigen_2d_burgers_outflow_firstorder_implicit) +set(exename ${testname}_exe) + +configure_file(../compare.py compare.py COPYONLY) +configure_file(gold.txt gold.txt COPYONLY) +configure_file(../plot.py plot.py COPYONLY) + +add_executable(${exename} ${CMAKE_CURRENT_SOURCE_DIR}/../main.cc) + +add_test(NAME ${testname} +COMMAND ${CMAKE_COMMAND} +-DMESHDRIVER=${MESHSRC}/create_full_mesh.py +-DOUTDIR=${CMAKE_CURRENT_BINARY_DIR} +-DEXENAME=$ +-DSTENCILVAL=3 +-P ${CMAKE_CURRENT_SOURCE_DIR}/../test.cmake +) diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/firstorder/gold.txt b/tests_cpp/eigen_2d_burgers_outflow_implicit/firstorder/gold.txt new file mode 100644 index 00000000..29ea6a61 --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/firstorder/gold.txt @@ -0,0 +1,800 @@ +2.873982525357717896e-05 +2.873982525357717896e-05 +9.518451723826902263e-05 +9.518451723826902263e-05 +2.748100603588286730e-04 +2.748100603588286730e-04 +6.898420427176647730e-04 +6.898420427176647730e-04 +1.500367128108482653e-03 +1.500367128108482653e-03 +2.823443253766142781e-03 +2.823443253766142781e-03 +4.613711001676902915e-03 +4.613711001676902915e-03 +6.606348232252464391e-03 +6.606348232252464391e-03 +8.386056372392937039e-03 +8.386056372392937039e-03 +9.526208910376336367e-03 +9.526208910376336367e-03 +9.710417656489738272e-03 +9.710417656489738272e-03 +8.842713975724484157e-03 +8.842713975724484157e-03 +7.122917846142971104e-03 +7.122917846142971104e-03 +5.006099646306888182e-03 +5.006099646306888182e-03 +3.035700161556099345e-03 +3.035700161556099345e-03 +1.583060375335185741e-03 +1.583060375335185741e-03 +7.132564019247373101e-04 +7.132564019247373101e-04 +2.797044747132363621e-04 +2.797044747132363621e-04 +9.595270342935998645e-05 +9.595270342935998645e-05 +2.888939192729877558e-05 +2.888939192729877558e-05 +7.299101790638417263e-05 +7.299101790638417263e-05 +2.410940502829049753e-04 +2.410940502829049753e-04 +6.914248354664422839e-04 +6.914248354664422839e-04 +1.712203068150484203e-03 +1.712203068150484203e-03 +3.640781883844851526e-03 +3.640781883844851526e-03 +6.648767095784709029e-03 +6.648767095784709029e-03 +1.052927963725216326e-02 +1.052927963725216326e-02 +1.470750584158317040e-02 +1.470750584158317040e-02 +1.844214827200072945e-02 +1.844214827200072945e-02 +2.101590791367392852e-02 +2.101590791367392852e-02 +2.182563182749713143e-02 +2.182563182749713143e-02 +2.051116346444867289e-02 +2.051116346444867289e-02 +1.719890276001106305e-02 +1.719890276001106305e-02 +1.258888919332892126e-02 +1.258888919332892126e-02 +7.866895093151702115e-03 +7.866895093151702115e-03 +4.151060628453971842e-03 +4.151060628453971842e-03 +1.862882414570770240e-03 +1.862882414570770240e-03 +7.233736877822215033e-04 +7.233736877822215033e-04 +2.461065532856013672e-04 +2.461065532856013672e-04 +7.373597191880513287e-05 +7.373597191880513287e-05 +1.619162163677340181e-04 +1.619162163677340181e-04 +5.320530437926035460e-04 +5.320530437926035460e-04 +1.506710618129291105e-03 +1.506710618129291105e-03 +3.640889328820760568e-03 +3.640889328820760568e-03 +7.454093597361336788e-03 +7.454093597361336788e-03 +1.299954772666421163e-02 +1.299954772666421163e-02 +1.969835095915133463e-02 +1.969835095915133810e-02 +2.662108719676254884e-02 +2.662108719676254884e-02 +3.281567818967875277e-02 +3.281567818967875277e-02 +3.742839416225328397e-02 +3.742839416225328397e-02 +3.966388214772299109e-02 +3.966388214772299109e-02 +3.876716768981489786e-02 +3.876716768981489786e-02 +3.440119465550674910e-02 +3.440119465550674910e-02 +2.693769294484439983e-02 +2.693769294484439983e-02 +1.791481737015239439e-02 +1.791481737015239439e-02 +9.790495122709425707e-03 +9.790495122709425707e-03 +4.388411340120931467e-03 +4.388411340120931467e-03 +1.669356529052715577e-03 +1.669356529052715360e-03 +5.574497953119461916e-04 +5.574497953119461916e-04 +1.652431740165662697e-04 +1.652431740165662697e-04 +3.138536123155618663e-04 +3.138536123155618663e-04 +1.023190008095733712e-03 +1.023190008095733712e-03 +2.844572417833858297e-03 +2.844572417833858297e-03 +6.649406098495009958e-03 +6.649406098495009958e-03 +1.299957336727236221e-02 +1.299957336727236221e-02 +2.157754374143133064e-02 +2.157754374143133064e-02 +3.135919392792172605e-02 +3.135919392792173299e-02 +4.119232010067142014e-02 +4.119232010067142707e-02 +5.007530624758425170e-02 +5.007530624758425863e-02 +5.714958938500840080e-02 +5.714958938500840080e-02 +6.157953383164312122e-02 +6.157953383164312122e-02 +6.233964633366768232e-02 +6.233964633366768232e-02 +5.847518933583577966e-02 +5.847518933583577966e-02 +4.943524681113868613e-02 +4.943524681113868613e-02 +3.591738957422349610e-02 +3.591738957422349610e-02 +2.107569397401294620e-02 +2.107569397401294620e-02 +9.586743056381845374e-03 +9.586743056381845374e-03 +3.510201536893694119e-03 +3.510201536893694119e-03 +1.124933936521681005e-03 +1.124933936521681005e-03 +3.261780234227465371e-04 +3.261780234227465371e-04 +5.317318670383966136e-04 +5.317318670383966136e-04 +1.716532924291051029e-03 +1.716532924291051029e-03 +4.669521818085370542e-03 +4.669521818085370542e-03 +1.053194531115930928e-02 +1.053194531115930928e-02 +1.969851698454277295e-02 +1.969851698454277295e-02 +3.135920652902691314e-02 +3.135920652902691314e-02 +4.415435180969642243e-02 +4.415435180969642243e-02 +5.685269791553066582e-02 +5.685269791553066582e-02 +6.847982145022933620e-02 +6.847982145022933620e-02 +7.822171232738263369e-02 +7.822171232738263369e-02 +8.529230625399511490e-02 +8.529230625399511490e-02 +8.864435728193642561e-02 +8.864435728193642561e-02 +8.691404701734552207e-02 +8.691404701734552207e-02 +7.865043528669900252e-02 +7.865043528669900252e-02 +6.282338211036585240e-02 +6.282338211036585240e-02 +4.100389163469370102e-02 +4.100389163469370102e-02 +1.983329039346224190e-02 +1.983329039346224190e-02 +6.953628767390820052e-03 +6.953628767390820052e-03 +2.049105401645301503e-03 +2.049105401645301503e-03 +5.685272794434416800e-04 +5.685272794434416800e-04 +7.882906212555158658e-04 +7.882906212555158658e-04 +2.520635574158126228e-03 +2.520635574158126228e-03 +6.722063494166427740e-03 +6.722063494166427740e-03 +1.471528456981613386e-02 +1.471528456981613386e-02 +2.662175952444097507e-02 +2.662175952444097507e-02 +4.119239098414981953e-02 +4.119239098414981953e-02 +5.685270629024663513e-02 +5.685270629024663513e-02 +7.234215468871942056e-02 +7.234215468871942056e-02 +8.671936101895622029e-02 +8.671936101895622029e-02 +9.921806956495843322e-02 +9.921806956495844709e-02 +1.091247022483874674e-01 +1.091247022483874535e-01 +1.155101194134125625e-01 +1.155101194134125625e-01 +1.169597355069472294e-01 +1.169597355069472294e-01 +1.115021413786817794e-01 +1.115021413786817794e-01 +9.661176996321826915e-02 +9.661176996321826915e-02 +7.071741894936955930e-02 +7.071741894936955930e-02 +3.843771075688261207e-02 +3.843771075688261207e-02 +1.351733930811111449e-02 +1.351733930811111449e-02 +3.457411149549197323e-03 +3.457411149549197323e-03 +8.768393408823488494e-04 +8.768393408823489579e-04 +1.024469824696957500e-03 +1.024469824696957500e-03 +3.255410691692508639e-03 +3.255410691692508639e-03 +8.575647610993098607e-03 +8.575647610993098607e-03 +1.845832510335571053e-02 +1.845832510335571053e-02 +3.281744540194538040e-02 +3.281744540194538040e-02 +5.007554183929467662e-02 +5.007554183929467662e-02 +6.847985739685481799e-02 +6.847985739685481799e-02 +8.671936682720306211e-02 +8.671936682720306211e-02 +1.038597212620602367e-01 +1.038597212620602367e-01 +1.191671172007175566e-01 +1.191671172007175566e-01 +1.319911836759235513e-01 +1.319911836759235513e-01 +1.415454961931899014e-01 +1.415454961931899014e-01 +1.466090182984615131e-01 +1.466090182984615131e-01 +1.451411366134933156e-01 +1.451411366134933156e-01 +1.338509701225856496e-01 +1.338509701225856496e-01 +1.081913713009177574e-01 +1.081913713009177436e-01 +6.742790903549861214e-02 +6.742790903549861214e-02 +2.596140075687665763e-02 +2.596140075687665763e-02 +5.683638013680356732e-03 +5.683638013680355865e-03 +1.200603655052420249e-03 +1.200603655052420249e-03 +1.169140521971979464e-03 +1.169140521971979464e-03 +3.712762116879750329e-03 +3.712762116879750329e-03 +9.773841340687235493e-03 +9.773841340687235493e-03 +2.104034378075574432e-02 +2.104034378075574432e-02 +3.743151131542591520e-02 +3.743151131542591520e-02 +5.715007938343958510e-02 +5.715007938343958510e-02 +7.822180149083812761e-02 +7.822180149083812761e-02 +9.921808737422396773e-02 +9.921808737422396773e-02 +1.191671207917373992e-01 +1.191671207917373715e-01 +1.373432027709245995e-01 +1.373432027709245995e-01 +1.531393636030465855e-01 +1.531393636030465855e-01 +1.658892505617783619e-01 +1.658892505617783619e-01 +1.746140670466522771e-01 +1.746140670466522771e-01 +1.775815225144266873e-01 +1.775815225144266873e-01 +1.714380953504801852e-01 +1.714380953504801852e-01 +1.497746779069446810e-01 +1.497746779069446810e-01 +1.055294432533241678e-01 +1.055294432533241678e-01 +4.700474608257837472e-02 +4.700474608257836778e-02 +9.613831405004649375e-03 +9.613831405004647640e-03 +1.477208741864330203e-03 +1.477208741864330203e-03 +1.172259265432595322e-03 +1.172259265432595322e-03 +3.742885157004168396e-03 +3.742885157004168396e-03 +9.968328507257488008e-03 +9.968328507257488008e-03 +2.185252257262234807e-02 +2.185252257262234807e-02 +3.966761266000910086e-02 +3.966761266000910779e-02 +6.158018828446369664e-02 +6.158018828446369664e-02 +8.529244174323030725e-02 +8.529244174323030725e-02 +1.091247335052700496e-01 +1.091247335052700496e-01 +1.319911913584337626e-01 +1.319911913584337626e-01 +1.531393654013509242e-01 +1.531393654013509242e-01 +1.719852103320660408e-01 +1.719852103320660408e-01 +1.879519514253438217e-01 +1.879519514253438217e-01 +2.002635423745368903e-01 +2.002635423745368903e-01 +2.076131370365280893e-01 +2.076131370365280893e-01 +2.071944695371971679e-01 +2.071944695371971679e-01 +1.920415858822571709e-01 +1.920415858822571709e-01 +1.493450568351995522e-01 +1.493450568351995522e-01 +7.683890521167344168e-02 +7.683890521167345555e-02 +1.650787558751332096e-02 +1.650787558751333137e-02 +1.689095293933452497e-03 +1.689095293933452497e-03 +1.031877449084909288e-03 +1.031877449084909288e-03 +3.327893491351383796e-03 +3.327893491351383796e-03 +9.055090787826160226e-03 +9.055090787826160226e-03 +2.053243058508896768e-02 +2.053243058508896768e-02 +3.877015473854902167e-02 +3.877015473854902861e-02 +6.234020409327012402e-02 +6.234020409327012402e-02 +8.864448491766446780e-02 +8.864448491766445393e-02 +1.155101527094306579e-01 +1.155101527094306579e-01 +1.415455055840548204e-01 +1.415455055840548204e-01 +1.658892532952233601e-01 +1.658892532952233601e-01 +1.879519521263174897e-01 +1.879519521263174897e-01 +2.072232870572120189e-01 +2.072232870572120189e-01 +2.230787598045296249e-01 +2.230787598045296249e-01 +2.345739925080181643e-01 +2.345739925080181643e-01 +2.398226555245457081e-01 +2.398226555245457081e-01 +2.324573289639490503e-01 +2.324573289639490503e-01 +1.950686791263761544e-01 +1.950686791263761544e-01 +1.121649280712701258e-01 +1.121649280712701258e-01 +2.624056464278549244e-02 +2.624056464278549244e-02 +1.855874179394101766e-03 +1.855874179394101766e-03 +7.960110845126047111e-04 +7.960110845126047111e-04 +2.598075333385254770e-03 +2.598075333385254770e-03 +7.259337858595617601e-03 +7.259337858595617601e-03 +1.721074465909781109e-02 +1.721074465909781109e-02 +3.440276795771838453e-02 +3.440276795771837759e-02 +5.847549194464776723e-02 +5.847549194464775335e-02 +8.691412270456867761e-02 +8.691412270456866374e-02 +1.169597578284046480e-01 +1.169597578284046480e-01 +1.466090255394078279e-01 +1.466090255394078279e-01 +1.746140695049198110e-01 +1.746140695049198110e-01 +2.002635431968111113e-01 +2.002635431968111113e-01 +2.230787600269563686e-01 +2.230787600269563409e-01 +2.425774737217365196e-01 +2.425774737217364641e-01 +2.581085646600990158e-01 +2.581085646600990158e-01 +2.685117438136923051e-01 +2.685117438136923051e-01 +2.691171569138732234e-01 +2.691171569138732234e-01 +2.393864497956572357e-01 +2.393864497956572635e-01 +1.468697860583690062e-01 +1.468697860583690340e-01 +3.486747978138997339e-02 +3.486747978138998033e-02 +1.810406215364018367e-03 +1.810406215364019017e-03 +5.370364921721321096e-04 +5.370364921721321096e-04 +1.771341744497385894e-03 +1.771341744497385894e-03 +5.073234873003343477e-03 +5.073234873003343477e-03 +1.259336380411173878e-02 +1.259336380411173878e-02 +2.693821692782434260e-02 +2.693821692782434260e-02 +4.943534986289626482e-02 +4.943534986289626482e-02 +7.865046417341728380e-02 +7.865046417341726992e-02 +1.115021512762264522e-01 +1.115021512762264522e-01 +1.451411403732800476e-01 +1.451411403732800476e-01 +1.775815240122357030e-01 +1.775815240122357030e-01 +2.076131376353408109e-01 +2.076131376353408109e-01 +2.345739927302898098e-01 +2.345739927302897543e-01 +2.581085647221513790e-01 +2.581085647221513235e-01 +2.775039786917873896e-01 +2.775039786917873341e-01 +2.917030980105513871e-01 +2.917030980105513871e-01 +3.013404987164187698e-01 +3.013404987164187698e-01 +2.808980624867751241e-01 +2.808980624867751796e-01 +1.692898996791745991e-01 +1.692898996791746269e-01 +3.374917818627542337e-02 +3.374917818627543031e-02 +1.199353212326397791e-03 +1.199353212326398225e-03 +3.164636245147167508e-04 +3.164636245147167508e-04 +1.050897393823448317e-03 +1.050897393823448317e-03 +3.060774537530709227e-03 +3.060774537530709227e-03 +7.867997666154346800e-03 +7.867997666154346800e-03 +1.791492126434729090e-02 +1.791492126434729090e-02 +3.591741099167507223e-02 +3.591741099167507223e-02 +6.282338933216155552e-02 +6.282338933216155552e-02 +9.661177300704749948e-02 +9.661177300704749948e-02 +1.338509715260565724e-01 +1.338509715260565724e-01 +1.714380960172517399e-01 +1.714380960172517399e-01 +2.071944698524973705e-01 +2.071944698524973705e-01 +2.398226556671150533e-01 +2.398226556671150256e-01 +2.685117438704012760e-01 +2.685117438704012760e-01 +2.917030980265991613e-01 +2.917030980265991058e-01 +3.117617136625963070e-01 +3.117617136625963070e-01 +3.376590428723901249e-01 +3.376590428723901249e-01 +3.102915466361798602e-01 +3.102915466361799712e-01 +1.507316568240995580e-01 +1.507316568240996690e-01 +1.759753657089992152e-02 +1.759753657089993886e-02 +4.405142112716267438e-04 +4.405142112716269064e-04 +1.628637770638636021e-04 +1.628637770638636021e-04 +5.423098632020054675e-04 +5.423098632020054675e-04 +1.590292070246753387e-03 +1.590292070246753387e-03 +4.151235251214945195e-03 +4.151235251214945195e-03 +9.790506451386415260e-03 +9.790506451386415260e-03 +2.107569652885251527e-02 +2.107569652885251527e-02 +4.100389280779209028e-02 +4.100389280779209722e-02 +7.071741961453627268e-02 +7.071741961453627268e-02 +1.081913716925064078e-01 +1.081913716925064078e-01 +1.497746781348633094e-01 +1.497746781348633094e-01 +1.920415860109393469e-01 +1.920415860109393469e-01 +2.324573290329800257e-01 +2.324573290329800257e-01 +2.691171569475998004e-01 +2.691171569475997449e-01 +3.013404987304033056e-01 +3.013404987304033056e-01 +3.376590428766200191e-01 +3.376590428766200747e-01 +3.613740464733033453e-01 +3.613740464733033453e-01 +2.645751509534173262e-01 +2.645751509534173818e-01 +7.065945601997783432e-02 +7.065945601997788983e-02 +2.637729311442317797e-03 +2.637729311442319098e-03 +1.680293702632584352e-04 +1.680293702632584352e-04 +7.324869501053931438e-05 +7.324869501053931438e-05 +2.439136107724817787e-04 +2.439136107724817787e-04 +7.149736777567683920e-04 +7.149736777567683920e-04 +1.862901414515245619e-03 +1.862901414515245619e-03 +4.388411973358616845e-03 +4.388411973358616845e-03 +9.586743208275457645e-03 +9.586743208275457645e-03 +1.983329050383641351e-02 +1.983329050383641351e-02 +3.843771085360912776e-02 +3.843771085360912776e-02 +6.742790911508014795e-02 +6.742790911508014795e-02 +1.055294433128078080e-01 +1.055294433128078080e-01 +1.493450568761576780e-01 +1.493450568761576780e-01 +1.950686791534719522e-01 +1.950686791534719522e-01 +2.393864498139401942e-01 +2.393864498139401942e-01 +2.808980624958546390e-01 +2.808980624958546390e-01 +3.102915466374367437e-01 +3.102915466374367437e-01 +2.645751509529743473e-01 +2.645751509529743473e-01 +1.057330141696317949e-01 +1.057330141696318226e-01 +7.209461018410067215e-03 +7.209461018410069817e-03 +2.625311877190467796e-04 +2.625311877190467796e-04 +7.400926428664477666e-05 +7.400926428664477666e-05 +2.881444435594164553e-05 +2.881444435594164553e-05 +9.585299182492693299e-05 +9.585299182492693299e-05 +2.800828484928022432e-04 +2.800828484928022432e-04 +7.233754436172108422e-04 +7.233754436172108422e-04 +1.669356548098330922e-03 +1.669356548098330922e-03 +3.510201538606926974e-03 +3.510201538606926974e-03 +6.953628765219156162e-03 +6.953628765219156162e-03 +1.351733929608290967e-02 +1.351733929608290967e-02 +2.596140072472026550e-02 +2.596140072472026203e-02 +4.700474602104129374e-02 +4.700474602104128680e-02 +7.683890508789568974e-02 +7.683890508789568974e-02 +1.121649278074332984e-01 +1.121649278074332984e-01 +1.468697856600742202e-01 +1.468697856600741924e-01 +1.692898993480654313e-01 +1.692898993480654035e-01 +1.507316566991467599e-01 +1.507316566991467321e-01 +7.065945599816894629e-02 +7.065945599816893241e-02 +7.209461017210366286e-03 +7.209461017210368021e-03 +3.085949365568110879e-04 +3.085949365568110879e-04 +9.665582052955351132e-05 +9.665582052955351132e-05 +2.896535004338501613e-05 +2.896535004338501613e-05 +9.919720688948248041e-06 +9.919720688948248041e-06 +3.296537694468057461e-05 +3.296537694468057461e-05 +9.604030039669556431e-05 +9.604030039669556431e-05 +2.461066162508814925e-04 +2.461066162508814925e-04 +5.574494660385840347e-04 +5.574494660385840347e-04 +1.124933161035150545e-03 +1.124933161035150545e-03 +2.049103856703332008e-03 +2.049103856703332008e-03 +3.457408707238642138e-03 +3.457408707238641705e-03 +5.683634998539646710e-03 +5.683634998539646710e-03 +9.613828068273325839e-03 +9.613828068273325839e-03 +1.650787081785519417e-02 +1.650787081785519070e-02 +2.624055620554653634e-02 +2.624055620554653634e-02 +3.486746856699922403e-02 +3.486746856699922403e-02 +3.374916999859873395e-02 +3.374916999859873395e-02 +1.759753391667905759e-02 +1.759753391667905412e-02 +2.637728870528548326e-03 +2.637728870528547893e-03 +2.625310582352186604e-04 +2.625310582352186604e-04 +9.665577487229198405e-05 +9.665577487229198405e-05 +3.305799867477848613e-05 +3.305799867477848613e-05 +9.950518956072408193e-06 +9.950518956072408193e-06 +2.989301098584977337e-06 +2.989301098584977337e-06 +9.927680931472465884e-06 +9.927680931472465884e-06 +2.886967493763391347e-05 +2.886967493763391347e-05 +7.363351099190124940e-05 +7.363351099190124940e-05 +1.650288558657967786e-04 +1.650288558657967786e-04 +3.258067070379667388e-04 +3.258067070379667388e-04 +5.680091865251434913e-04 +5.680091865251434913e-04 +8.762743685381023991e-04 +8.762743685381023991e-04 +1.200125112387124583e-03 +1.200125112387124583e-03 +1.476848071725970655e-03 +1.476848071725970655e-03 +1.688731901793414638e-03 +1.688731901793414638e-03 +1.855375866019437622e-03 +1.855375866019437622e-03 +1.809797494126035369e-03 +1.809797494126035369e-03 +1.198801557605039389e-03 +1.198801557605039389e-03 +4.401349738458046812e-04 +4.401349738458046269e-04 +1.678131055669875526e-04 +1.678131055669875526e-04 +7.390620292520634919e-05 +7.390620292520634919e-05 +2.892370277705373573e-05 +2.892370277705373573e-05 +9.936027232859333486e-06 +9.936027232859333486e-06 +2.996263814266988420e-06 +2.996263814266988420e-06 +7.885662595402906146e-07 +7.885662595402906146e-07 +2.618043009336446800e-06 +2.618043009336446800e-06 +7.606402243416762565e-06 +7.606402243416762565e-06 +1.935744263419941931e-05 +1.935744263419941931e-05 +4.317732465290848265e-05 +4.317732465290848265e-05 +8.447239664758985484e-05 +8.447239664758985484e-05 +1.450138302030300902e-04 +1.450138302030300902e-04 +2.183485535690893212e-04 +2.183485535690893212e-04 +2.878685314670370891e-04 +2.878685314670370891e-04 +3.314117564275045522e-04 +3.314117564275045522e-04 +3.323712941417572087e-04 +3.323712941417572087e-04 +2.901662479352816936e-04 +2.901662479352816936e-04 +2.205138907228784923e-04 +2.205138907228784923e-04 +1.459423275069515779e-04 +1.459423275069515779e-04 +8.470455514567227801e-05 +8.470455514567227801e-05 +4.324678265665701853e-05 +4.324678265665701853e-05 +1.937579710520107211e-05 +1.937579710520107211e-05 +7.610128104894790865e-06 +7.610128104894790865e-06 +2.618623909600831834e-06 +2.618623909600831834e-06 +7.902092518265551476e-07 +7.902092518265551476e-07 +1.823829432661276714e-07 +1.823829432661276714e-07 +6.054324388772379224e-07 +6.054324388772379224e-07 +1.758367701267488359e-06 +1.758367701267488359e-06 +4.470901254621804705e-06 +4.470901254621804705e-06 +9.954087888303038127e-06 +9.954087888303038127e-06 +1.940949665430117800e-05 +1.940949665430117800e-05 +3.314874007852850520e-05 +3.314874007852850520e-05 +4.957547991972452450e-05 +4.957547991972452450e-05 +6.488706005392352156e-05 +6.488706005392352156e-05 +7.426542406889497091e-05 +7.426542406889497091e-05 +7.427975061048183562e-05 +7.427975061048183562e-05 +6.492016138602213558e-05 +6.492016138602213558e-05 +4.960840950866481525e-05 +4.960840950866481525e-05 +3.317028165434255158e-05 +3.317028165434255158e-05 +1.941971858674511545e-05 +1.941971858674511545e-05 +9.957715942325221677e-06 +9.957715942325221677e-06 +4.471875145462200742e-06 +4.471875145462200742e-06 +1.758566644689487472e-06 +1.758566644689487472e-06 +6.054638239274369353e-07 +6.054638239274369353e-07 +1.827505789116900402e-07 +1.827505789116900402e-07 diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/main.cc b/tests_cpp/eigen_2d_burgers_outflow_implicit/main.cc new file mode 100644 index 00000000..c6239496 --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/main.cc @@ -0,0 +1,45 @@ + +#include "pressio/ode_steppers_implicit.hpp" +#include "pressio/ode_advancers.hpp" +#include "pressiodemoapps/advection_diffusion2d.hpp" +#include "../observer.hpp" + +int main() +{ + + namespace pda = pressiodemoapps; + const auto meshObj = pda::load_cellcentered_uniform_mesh_eigen("."); +#ifdef USE_WENO5 + const auto scheme = pda::InviscidFluxReconstruction::Weno5; +#elif defined USE_WENO3 + const auto scheme = pda::InviscidFluxReconstruction::Weno3; +#else + const auto scheme = pda::InviscidFluxReconstruction::FirstOrder; +#endif + + const auto schemeVisc = pda::ViscousFluxReconstruction::FirstOrder; + + const auto probId = pda::AdvectionDiffusion2d::BurgersOutflow; + auto appObj = pda::create_problem_eigen(meshObj, probId, scheme, schemeVisc); + + using app_t = decltype(appObj); + using state_t = typename app_t::state_type; + using jacob_t = typename app_t::jacobian_type; + state_t state = appObj.initialCondition(); + + auto stepperObj = pressio::ode::create_implicit_stepper( + pressio::ode::StepScheme::CrankNicolson, appObj); + + using lin_solver_t = pressio::linearsolvers::Solver< + pressio::linearsolvers::iterative::Bicgstab, jacob_t>; + lin_solver_t linSolverObj; + auto NonLinSolver=pressio::create_newton_solver(stepperObj, linSolverObj); + NonLinSolver.setStopTolerance(1e-5); + + const auto dt = 0.01; + const auto Nsteps = pressio::ode::StepCount(2./dt); + FomObserver Obs("burgers2d_solution.bin", 50); + pressio::ode::advance_n_steps(stepperObj, state, 0., dt, Nsteps, Obs, NonLinSolver); + + return 0; +} diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/plot.py b/tests_cpp/eigen_2d_burgers_outflow_implicit/plot.py new file mode 100644 index 00000000..231de8f2 --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/plot.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python3 + +import matplotlib.pyplot as plt +from matplotlib import cm +import numpy as np +from numpy import linalg as LA +import re + +def extractN(ns): + reg = re.compile(r''+ns+'.+') + file1 = open('info.dat', 'r') + strings = re.search(reg, file1.read()) + file1.close() + assert(strings) + return int(strings.group().split()[1]) + +########################## +if __name__== "__main__": +########################## + nx = extractN('nx') + ny = extractN('ny') + print(nx, ny) + fomTotDofs = nx*ny*2 + + fomCoords = np.loadtxt('coordinates.dat', dtype=float) + x_fom, y_fom = fomCoords[:,1], fomCoords[:,2] + x_fom = np.reshape(x_fom, (ny,nx)) + y_fom = np.reshape(y_fom, (ny,nx)) + + fomTestD = np.fromfile("burgers2d_solution.bin") + nt = int(np.size(fomTestD)/fomTotDofs) + print("fomTest: nt = ", nt) + fomTestD = np.reshape(fomTestD, (nt, fomTotDofs)) + + fig = plt.figure(1) + for i in range(0, nt): + fomS = fomTestD[i,:] + fomS = np.reshape(fomS, (ny*nx, 2)) + fomS = np.reshape(fomS[:,0], (nx,ny)) + plt.clf() + ax = plt.gca() + h = plt.contourf(x_fom, y_fom, fomS) + ax.set_aspect(aspect=1.) + plt.colorbar() + plt.pause(0.001) + # plt.show() diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/test.cmake b/tests_cpp/eigen_2d_burgers_outflow_implicit/test.cmake new file mode 100644 index 00000000..1275f0bb --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/test.cmake @@ -0,0 +1,24 @@ +include(FindUnixCommands) + +set(CMD "python3 ${MESHDRIVER} -n 20 20 --outDir ${OUTDIR} -s ${STENCILVAL} --bounds -1.0 1.0 -1.0 1.0") +execute_process(COMMAND ${BASH} -c ${CMD} RESULT_VARIABLE RES) +if(RES) + message(FATAL_ERROR "Mesh generation failed") +else() + message("Mesh generation succeeded!") +endif() + +execute_process(COMMAND ${EXENAME} RESULT_VARIABLE CMD_RESULT) +if(RES) + message(FATAL_ERROR "run failed") +else() + message("run succeeded!") +endif() + +set(CMD "python3 compare.py") +execute_process(COMMAND ${BASH} -c ${CMD} RESULT_VARIABLE RES) +if(RES) + message(FATAL_ERROR "comparison failed") +else() + message("comparison succeeded!") +endif() diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/weno3/CMakeLists.txt b/tests_cpp/eigen_2d_burgers_outflow_implicit/weno3/CMakeLists.txt new file mode 100644 index 00000000..b9a114b4 --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/weno3/CMakeLists.txt @@ -0,0 +1,19 @@ + +set(testname eigen_2d_burgers_outflow_weno3_implicit) +set(exename ${testname}_exe) + +configure_file(../compare.py compare.py COPYONLY) +configure_file(gold.txt gold.txt COPYONLY) +configure_file(../plot.py plot.py COPYONLY) + +add_executable(${exename} ${CMAKE_CURRENT_SOURCE_DIR}/../main.cc) +target_compile_definitions(${exename} PUBLIC -DUSE_WENO3) + +add_test(NAME ${testname} +COMMAND ${CMAKE_COMMAND} +-DMESHDRIVER=${MESHSRC}/create_full_mesh.py +-DOUTDIR=${CMAKE_CURRENT_BINARY_DIR} +-DEXENAME=$ +-DSTENCILVAL=5 +-P ${CMAKE_CURRENT_SOURCE_DIR}/../test.cmake +) diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/weno3/gold.txt b/tests_cpp/eigen_2d_burgers_outflow_implicit/weno3/gold.txt new file mode 100644 index 00000000..b8cd3325 --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/weno3/gold.txt @@ -0,0 +1,800 @@ +2.872340186914654436e-05 +2.872340186914654436e-05 +9.503371622176118934e-05 +9.503371622176118934e-05 +2.736638066969191698e-04 +2.736638066969191698e-04 +6.836853189277268495e-04 +6.836853189277268495e-04 +1.475077367394037053e-03 +1.475077367394037053e-03 +2.741859063866173771e-03 +2.741859063866173771e-03 +4.415224292534288171e-03 +4.415224292534288171e-03 +6.235438338296738110e-03 +6.235438338296738110e-03 +7.834301295433171095e-03 +7.834301295433171095e-03 +8.859773213057376903e-03 +8.859773213057376903e-03 +8.976340063147696460e-03 +8.976340063147696460e-03 +8.255025052703977503e-03 +8.255025052703977503e-03 +6.734701183216371817e-03 +6.734701183216371817e-03 +4.777821600789483938e-03 +4.777821600789483938e-03 +2.924430612614365484e-03 +2.924430612614365484e-03 +1.540881771891522509e-03 +1.540881771891522509e-03 +7.012066306563262148e-04 +7.012066306563262148e-04 +2.770679887837833484e-04 +2.770679887837833484e-04 +9.549390979845627373e-05 +9.549390979845627373e-05 +2.882619817514248288e-05 +2.882619817514248288e-05 +7.289957889827519027e-05 +7.289957889827519027e-05 +2.403025822925897062e-04 +2.403025822925897062e-04 +6.870280598577883317e-04 +6.870280598577883317e-04 +1.695313684748896415e-03 +1.695313684748896415e-03 +3.577784524506906309e-03 +3.577784524506906309e-03 +6.470040839552547961e-03 +6.470040839552547961e-03 +1.016489592116763020e-02 +1.016489592116763020e-02 +1.414496187640165714e-02 +1.414496187640165714e-02 +1.768525800639152851e-02 +1.768525800639152851e-02 +2.037165950048302684e-02 +2.037165950048302684e-02 +2.081283471153486109e-02 +2.081283471153486109e-02 +1.975316774982575907e-02 +1.975316774982575907e-02 +1.675911163151006661e-02 +1.675911163151006661e-02 +1.219811506063076499e-02 +1.219811506063076499e-02 +7.541861092850325916e-03 +7.541861092850325916e-03 +3.978597375643455093e-03 +3.978597375643455093e-03 +1.802586494639276922e-03 +1.802586494639276922e-03 +7.086543341844778014e-04 +7.086543341844778014e-04 +2.435167909453695475e-04 +2.435167909453695475e-04 +7.335442776350058148e-05 +7.335442776350058148e-05 +1.614679488491355380e-04 +1.614679488491355380e-04 +5.292050462580282638e-04 +5.292050462580282638e-04 +1.492733316017893922e-03 +1.492733316017893922e-03 +3.569766212889014401e-03 +3.569766212889014401e-03 +7.191474570551842534e-03 +7.191474570551842534e-03 +1.235101271229600145e-02 +1.235101271229600145e-02 +1.854669225628418835e-02 +1.854669225628418835e-02 +2.504274286320257170e-02 +2.504274286320257170e-02 +3.067992302346804132e-02 +3.067992302346804825e-02 +3.573254471794228698e-02 +3.573254471794228004e-02 +3.706381318431191396e-02 +3.706381318431190008e-02 +3.617128089950789188e-02 +3.617128089950788494e-02 +3.258648585821055654e-02 +3.258648585821055654e-02 +2.550232720483260399e-02 +2.550232720483260745e-02 +1.653244771208326955e-02 +1.653244771208326608e-02 +8.917088106113713400e-03 +8.917088106113713400e-03 +4.064143259353041743e-03 +4.064143259353041743e-03 +1.592824384352219602e-03 +1.592824384352219602e-03 +5.448914977682216450e-04 +5.448914977682216450e-04 +1.635042311972274582e-04 +1.635042311972274582e-04 +3.123394228607603157e-04 +3.123394228607603157e-04 +1.016921740875942901e-03 +1.016921740875942901e-03 +2.803529591411414550e-03 +2.803529591411414550e-03 +6.443360276280795523e-03 +6.443360276280795523e-03 +1.234416541207333448e-02 +1.234416541207333448e-02 +2.021012090543007622e-02 +2.021012090543007622e-02 +2.918599008091407598e-02 +2.918599008091407598e-02 +3.868055418122373634e-02 +3.868055418122374328e-02 +4.625904106773792712e-02 +4.625904106773792712e-02 +5.475787338323014758e-02 +5.475787338323013370e-02 +5.797591427638688205e-02 +5.797591427638688205e-02 +5.746198027024760485e-02 +5.746198027024759791e-02 +5.431668995068502509e-02 +5.431668995068501815e-02 +4.621692130227948725e-02 +4.621692130227948725e-02 +3.257392879010995862e-02 +3.257392879010995862e-02 +1.802217647656975674e-02 +1.802217647656975674e-02 +8.197246263766671745e-03 +8.197246263766671745e-03 +3.172430345147924085e-03 +3.172430345147924085e-03 +1.072736377465606430e-03 +1.072736377465606430e-03 +3.196814615666886150e-04 +3.196814615666886150e-04 +5.283585913230372331e-04 +5.283585913230372331e-04 +1.704190545738113946e-03 +1.704190545738113946e-03 +4.575595642600381770e-03 +4.575595642600381770e-03 +1.010518685414016002e-02 +1.010518685414016002e-02 +1.852341063020004203e-02 +1.852341063020004203e-02 +2.918481544901341956e-02 +2.918481544901341956e-02 +4.102256220956653626e-02 +4.102256220956653626e-02 +5.291231795079547051e-02 +5.291231795079547051e-02 +6.357633875948665025e-02 +6.357633875948666413e-02 +7.301938538218026575e-02 +7.301938538218025188e-02 +8.138510364208428882e-02 +8.138510364208427494e-02 +8.211178655807753468e-02 +8.211178655807754856e-02 +8.030995790746689456e-02 +8.030995790746689456e-02 +7.316450966861615113e-02 +7.316450966861613725e-02 +5.752641158159096241e-02 +5.752641158159096241e-02 +3.394955501565876527e-02 +3.394955501565876527e-02 +1.507043765835094223e-02 +1.507043765835094223e-02 +5.647837388019499513e-03 +5.647837388019499513e-03 +1.864381373183205065e-03 +1.864381373183205065e-03 +5.479127906684417518e-04 +5.479127906684417518e-04 +7.824930372726533788e-04 +7.824930372726533788e-04 +2.499974633517328300e-03 +2.499974633517328300e-03 +6.561233787683196962e-03 +6.561233787683196962e-03 +1.404503008089030662e-02 +1.404503008089030662e-02 +2.498560768770746038e-02 +2.498560768770746038e-02 +3.867438382954749287e-02 +3.867438382954749287e-02 +5.291376646699646330e-02 +5.291376646699646330e-02 +7.054355449308523374e-02 +7.054355449308521986e-02 +8.117468553241330431e-02 +8.117468553241327656e-02 +9.589001407612221528e-02 +9.589001407612217365e-02 +1.046903251896006137e-01 +1.046903251896005999e-01 +1.088895771666096163e-01 +1.088895771666096302e-01 +1.083166999404693021e-01 +1.083166999404693159e-01 +1.037135625516276000e-01 +1.037135625516276000e-01 +9.000031139722490525e-02 +9.000031139722491913e-02 +6.117882292872835193e-02 +6.117882292872835887e-02 +2.625461141643805132e-02 +2.625461141643805479e-02 +9.098494898999050923e-03 +9.098494898999052657e-03 +2.866167596076472229e-03 +2.866167596076472229e-03 +8.240523508639703955e-04 +8.240523508639703955e-04 +1.016131835596077553e-03 +1.016131835596077553e-03 +3.228742920749622970e-03 +3.228742920749622970e-03 +8.362044354872467250e-03 +8.362044354872467250e-03 +1.755234799977231988e-02 +1.755234799977231988e-02 +3.057816832186041958e-02 +3.057816832186041958e-02 +4.623837516144367610e-02 +4.623837516144367610e-02 +6.358020909984013480e-02 +6.358020909984013480e-02 +8.117728315017914975e-02 +8.117728315017914975e-02 +9.944627069482142590e-02 +9.944627069482142590e-02 +1.122497666433840829e-01 +1.122497666433841107e-01 +1.336512873137510049e-01 +1.336512873137510604e-01 +1.391489457596560009e-01 +1.391489457596560286e-01 +1.393349594226837296e-01 +1.393349594226837018e-01 +1.371299607166312773e-01 +1.371299607166313050e-01 +1.272764558610018637e-01 +1.272764558610018915e-01 +9.994683941416608597e-02 +9.994683941416608597e-02 +4.634733427213014179e-02 +4.634733427213014179e-02 +1.347461527214672064e-02 +1.347461527214672064e-02 +3.894671158227016716e-03 +3.894671158227016716e-03 +1.086709291616683214e-03 +1.086709291616683214e-03 +1.159192693890569647e-03 +1.159192693890569647e-03 +3.683881710216504058e-03 +3.683881710216504058e-03 +9.560006779147573949e-03 +9.560006779147573949e-03 +2.022132362507709363e-02 +2.022132362507709363e-02 +3.557297582774773814e-02 +3.557297582774773814e-02 +5.471528105718446305e-02 +5.471528105718444918e-02 +7.301447549445044816e-02 +7.301447549445044816e-02 +9.589554024669337540e-02 +9.589554024669341703e-02 +1.122512110042423045e-01 +1.122512110042422490e-01 +1.297310611919354484e-01 +1.297310611919353929e-01 +1.459469380120116677e-01 +1.459469380120116955e-01 +1.681296805529282345e-01 +1.681296805529282623e-01 +1.731640601003850422e-01 +1.731640601003850422e-01 +1.714094141591715470e-01 +1.714094141591716025e-01 +1.654129004608445863e-01 +1.654129004608445863e-01 +1.446096742386915601e-01 +1.446096742386915601e-01 +8.359383107241911970e-02 +8.359383107241911970e-02 +1.950658379612308507e-02 +1.950658379612308507e-02 +4.666935839643961657e-03 +4.666935839643961657e-03 +1.252935227969043108e-03 +1.252935227969043108e-03 +1.162611199309747552e-03 +1.162611199309747552e-03 +3.714620558941753199e-03 +3.714620558941753199e-03 +9.696195700540831408e-03 +9.696195700540831408e-03 +2.066136134607747532e-02 +2.066136134607747532e-02 +3.688932084783497906e-02 +3.688932084783497906e-02 +5.791215714533537451e-02 +5.791215714533537451e-02 +8.137786236140506779e-02 +8.137786236140506779e-02 +1.046956640730065530e-01 +1.046956640730065807e-01 +1.336552947297852401e-01 +1.336552947297851568e-01 +1.459483683661985942e-01 +1.459483683661985942e-01 +1.657300745854751733e-01 +1.657300745854751456e-01 +1.812640573982127190e-01 +1.812640573982127190e-01 +1.991976086242382427e-01 +1.991976086242382427e-01 +2.034579852073410133e-01 +2.034579852073409578e-01 +2.013956011472494700e-01 +2.013956011472494423e-01 +1.891826196622281575e-01 +1.891826196622281298e-01 +1.363307463481452930e-01 +1.363307463481452930e-01 +3.289729207532702399e-02 +3.289729207532703786e-02 +4.963230958038794896e-03 +4.963230958038794896e-03 +1.258509508911925291e-03 +1.258509508911925291e-03 +1.024161483898132280e-03 +1.024161483898132280e-03 +3.307066622378775252e-03 +3.307066622378775252e-03 +8.869110361388031841e-03 +8.869110361388031841e-03 +1.960616629666059318e-02 +1.960616629666059318e-02 +3.600987233648447056e-02 +3.600987233648447056e-02 +5.740196372463032826e-02 +5.740196372463032826e-02 +8.210499134320792258e-02 +8.210499134320792258e-02 +1.088955199657752754e-01 +1.088955199657753031e-01 +1.391548352957413248e-01 +1.391548352957413526e-01 +1.681322833674999406e-01 +1.681322833674999129e-01 +1.812651307284526003e-01 +1.812651307284525448e-01 +2.101874607502330350e-01 +2.101874607502330072e-01 +2.198458183210821870e-01 +2.198458183210821593e-01 +2.312728747144571495e-01 +2.312728747144571217e-01 +2.316869680866294190e-01 +2.316869680866294190e-01 +2.262458749754666099e-01 +2.262458749754666099e-01 +1.915902697952102662e-01 +1.915902697952102662e-01 +6.844912372541800405e-02 +6.844912372541800405e-02 +5.450133809269877040e-03 +5.450133809269877040e-03 +1.101519818175843504e-03 +1.101519818175843504e-03 +7.905923703832496872e-04 +7.905923703832496872e-04 +2.582783009528992566e-03 +2.582783009528992566e-03 +7.139419895782943629e-03 +7.139419895782943629e-03 +1.663335771800429977e-02 +1.663335771800429977e-02 +3.246728758825490752e-02 +3.246728758825490752e-02 +5.427285346291520901e-02 +5.427285346291520901e-02 +8.030592146248492769e-02 +8.030592146248492769e-02 +1.083224740571800543e-01 +1.083224740571800820e-01 +1.393405002642600643e-01 +1.393405002642600365e-01 +1.731667034361413227e-01 +1.731667034361414059e-01 +1.991989828909614957e-01 +1.991989828909615234e-01 +2.198463231210124880e-01 +2.198463231210125157e-01 +2.583319367746053086e-01 +2.583319367746053086e-01 +2.584689605965893588e-01 +2.584689605965893588e-01 +2.607107326281520310e-01 +2.607107326281520865e-01 +2.617064935877325138e-01 +2.617064935877325138e-01 +2.436381960401965974e-01 +2.436381960401965974e-01 +1.256651850996303443e-01 +1.256651850996303443e-01 +7.943945888664259367e-03 +7.943945888664259367e-03 +8.463624427084775418e-04 +8.463624427084775418e-04 +5.336031539895435627e-04 +5.336031539895435627e-04 +1.759009512168862189e-03 +1.759009512168862189e-03 +4.975850299303851212e-03 +4.975850299303851212e-03 +1.211641910507626693e-02 +1.211641910507626693e-02 +2.544556776053265301e-02 +2.544556776053265301e-02 +4.619883544530085928e-02 +4.619883544530085928e-02 +7.316497938650552701e-02 +7.316497938650552701e-02 +1.037187201716345819e-01 +1.037187201716345958e-01 +1.371345262271489729e-01 +1.371345262271489729e-01 +1.714117809110042301e-01 +1.714117809110042301e-01 +2.034588453072114778e-01 +2.034588453072114778e-01 +2.312733087937207110e-01 +2.312733087937207943e-01 +2.584690747928307664e-01 +2.584690747928308219e-01 +3.007113899750168406e-01 +3.007113899750168406e-01 +3.064169061568519647e-01 +3.064169061568519647e-01 +3.004478869615822001e-01 +3.004478869615822001e-01 +2.868257442179072414e-01 +2.868257442179072414e-01 +2.002480599454643795e-01 +2.002480599454644072e-01 +1.779862631734403441e-02 +1.779862631734404482e-02 +5.870204880368896016e-04 +5.870204880368896016e-04 +3.147277081680782121e-04 +3.147277081680782121e-04 +1.042404972244733970e-03 +1.042404972244733970e-03 +2.997498901140243071e-03 +2.997498901140243071e-03 +7.505878653663883120e-03 +7.505878653663883120e-03 +1.651808004830011833e-02 +1.651808004830011833e-02 +3.257294957993620010e-02 +3.257294957993620010e-02 +5.752893176817616694e-02 +5.752893176817617388e-02 +9.000350198956109937e-02 +9.000350198956111325e-02 +1.272791893161914489e-01 +1.272791893161914767e-01 +1.654145340150749433e-01 +1.654145340150749433e-01 +2.013963449001563433e-01 +2.013963449001563710e-01 +2.316873292598280176e-01 +2.316873292598280176e-01 +2.607108635040313316e-01 +2.607108635040312206e-01 +3.064168540401578378e-01 +3.064168540401576712e-01 +3.210995861791196559e-01 +3.210995861791196559e-01 +3.304923994516762997e-01 +3.304923994516763552e-01 +3.500557044886096203e-01 +3.500557044886097313e-01 +2.449307137395454192e-01 +2.449307137395454470e-01 +1.635581696571078089e-02 +1.635581696571078089e-02 +3.411754664322945202e-04 +3.411754664322945202e-04 +1.622648989579520405e-04 +1.622648989579520405e-04 +5.378249410532217906e-04 +5.378249410532217906e-04 +1.560946321250956514e-03 +1.560946321250956514e-03 +3.968186647472144758e-03 +3.968186647472144758e-03 +8.914822858613755813e-03 +8.914822858613755813e-03 +1.802258619466899983e-02 +1.802258619466899983e-02 +3.395023196160790874e-02 +3.395023196160790180e-02 +6.117929025123992198e-02 +6.117929025123992892e-02 +9.994747039964534119e-02 +9.994747039964535507e-02 +1.446103117989511744e-01 +1.446103117989511744e-01 +1.891830526279451863e-01 +1.891830526279451863e-01 +2.262461173979958551e-01 +2.262461173979958551e-01 +2.617065646220276820e-01 +2.617065646220276820e-01 +3.004478576219950514e-01 +3.004478576219950514e-01 +3.304923615909618473e-01 +3.304923615909617918e-01 +3.907843126682680190e-01 +3.907843126682680746e-01 +3.937718916138377434e-01 +3.937718916138377434e-01 +1.944147639351959389e-01 +1.944147639351959389e-01 +4.349197479748864725e-03 +4.349197479748864725e-03 +1.656410457378680373e-04 +1.656410457378680373e-04 +7.309473690716465828e-05 +7.309473690716465828e-05 +2.424378594455041272e-04 +2.424378594455041272e-04 +7.047285033990149158e-04 +7.047285033990149158e-04 +1.800960849943555962e-03 +1.800960849943555962e-03 +4.063968234616795010e-03 +4.063968234616795010e-03 +8.197284670643779220e-03 +8.197284670643779220e-03 +1.507046669001023236e-02 +1.507046669001023236e-02 +2.625453117718182586e-02 +2.625453117718182239e-02 +4.634706705478338490e-02 +4.634706705478339184e-02 +8.359368316757800899e-02 +8.359368316757800899e-02 +1.363307891829103891e-01 +1.363307891829103613e-01 +1.915903510981476487e-01 +1.915903510981476765e-01 +2.436382181181556139e-01 +2.436382181181556694e-01 +2.868257299815103467e-01 +2.868257299815103467e-01 +3.500556506753876307e-01 +3.500556506753876862e-01 +3.937718690698299717e-01 +3.937718690698299717e-01 +3.243689173803241399e-01 +3.243689173803241399e-01 +3.406837207512914550e-02 +3.406837207512914550e-02 +2.918265366534085817e-04 +2.918265366534085817e-04 +7.356081891974136579e-05 +7.356081891974136579e-05 +2.878183745031362198e-05 +2.878183745031362198e-05 +9.552659700789000309e-05 +9.552659700789000309e-05 +2.775808717671619124e-04 +2.775808717671619124e-04 +7.085096333400628648e-04 +7.085096333400628648e-04 +1.592818717373823724e-03 +1.592818717373823724e-03 +3.172434316021055696e-03 +3.172434316021055696e-03 +5.647850893493710708e-03 +5.647850893493710708e-03 +9.098528695521385101e-03 +9.098528695521385101e-03 +1.347467414500176132e-02 +1.347467414500176305e-02 +1.950665026562079310e-02 +1.950665026562079310e-02 +3.289729591650134011e-02 +3.289729591650133317e-02 +6.844910433844786368e-02 +6.844910433844786368e-02 +1.256651709004969475e-01 +1.256651709004969752e-01 +2.002479903120474303e-01 +2.002479903120474303e-01 +2.449306453863263622e-01 +2.449306453863264454e-01 +1.944147576326029359e-01 +1.944147576326029636e-01 +3.406837078005506281e-02 +3.406837078005505587e-02 +3.705685450209890539e-04 +3.705685450209890539e-04 +9.600415944829897644e-05 +9.600415944829897644e-05 +2.888477147711751710e-05 +2.888477147711751710e-05 +9.914279787327621058e-06 +9.914279787327621058e-06 +3.290932911471119142e-05 +3.290932911471119142e-05 +9.559268071419554404e-05 +9.559268071419554404e-05 +2.434973769309800974e-04 +2.434973769309800974e-04 +5.448585062943892866e-04 +5.448585062943892866e-04 +1.072601844357821358e-03 +1.072601844357821358e-03 +1.864021586280514581e-03 +1.864021586280514581e-03 +2.865592414997476627e-03 +2.865592414997476627e-03 +3.893909957995662918e-03 +3.893909957995662918e-03 +4.666072565865747764e-03 +4.666072565865747764e-03 +4.962432131613385249e-03 +4.962432131613385249e-03 +5.449547418337070612e-03 +5.449547418337070612e-03 +7.943436500071312384e-03 +7.943436500071314119e-03 +1.779799915920607031e-02 +1.779799915920607725e-02 +1.635461733984533775e-02 +1.635461733984535510e-02 +4.347776236650361011e-03 +4.347776236650361878e-03 +2.918135486684698488e-04 +2.918135486684698488e-04 +9.600349278731042416e-05 +9.600349278731042416e-05 +3.296242638936149548e-05 +3.296242638936149548e-05 +9.939424378931682237e-06 +9.939424378931682237e-06 +2.988601987316248402e-06 +2.988601987316248402e-06 +9.920327280683448418e-06 +9.920327280683448418e-06 +2.880928885667598295e-05 +2.880928885667598295e-05 +7.326318781081628967e-05 +7.326318781081628967e-05 +1.633686847550297765e-04 +1.633686847550297765e-04 +3.197441030623565075e-04 +3.197441030623565075e-04 +5.490586746553867328e-04 +5.490586746553867328e-04 +8.274262064250647881e-04 +8.274262064250647881e-04 +1.093002018242471805e-03 +1.093002018242471805e-03 +1.261447233189359925e-03 +1.261447233189359925e-03 +1.267033307655745589e-03 +1.267033307655745589e-03 +1.107745381565444429e-03 +1.107745381565444429e-03 +8.494663299082708422e-04 +8.494663299082708422e-04 +5.877255640199390311e-04 +5.877255640199390311e-04 +3.410895532096535287e-04 +3.410895532096535287e-04 +1.655468338569301230e-04 +1.655468338569301230e-04 +7.347044708653211582e-05 +7.347044708653211582e-05 +2.884491700633570014e-05 +2.884491700633570014e-05 +9.925154114133244695e-06 +9.925154114133244695e-06 +2.995052844888061068e-06 +2.995052844888061068e-06 +7.884975930434087479e-07 +7.884975930434087479e-07 +2.617310956375428986e-06 +2.617310956375428986e-06 +7.600327080210691689e-06 +7.600327080210691689e-06 +1.931890898934767326e-05 +1.931890898934767326e-05 +4.299312990030456310e-05 +4.299312990030456310e-05 +8.382367323520377802e-05 +8.382367323520377802e-05 +1.432687851615411416e-04 +1.432687851615411416e-04 +2.145244519072102827e-04 +2.145244519072102827e-04 +2.811140138506892645e-04 +2.811140138506892645e-04 +3.220817850439932403e-04 +3.220817850439932403e-04 +3.223741720986069586e-04 +3.223741720986069586e-04 +2.817937563431330522e-04 +2.817937563431330522e-04 +2.152047465755699827e-04 +2.152047465755699827e-04 +1.437098539980276790e-04 +1.437098539980276790e-04 +8.401730087581817553e-05 +8.401730087581817553e-05 +4.305063488548597009e-05 +4.305063488548597009e-05 +1.933264489797818320e-05 +1.933264489797818320e-05 +7.602819359828572604e-06 +7.602819359828572604e-06 +2.617648832003952086e-06 +2.617648832003952086e-06 +7.901049517367656961e-07 +7.901049517367656961e-07 +1.823777696231993586e-07 +1.823777696231993586e-07 +6.053765949213519006e-07 +6.053765949213519006e-07 +1.757901891818919833e-06 +1.757901891818919833e-06 +4.467923017837591905e-06 +4.467923017837591905e-06 +9.939495171188612887e-06 +9.939495171188612887e-06 +1.935495142054928958e-05 +1.935495142054928958e-05 +3.299435709305326364e-05 +3.299435709305326364e-05 +4.924497836989638117e-05 +4.924497836989638117e-05 +6.434454308033667517e-05 +6.434454308033667517e-05 +7.357174501377471620e-05 +7.357174501377471620e-05 +7.358631144388814516e-05 +7.358631144388814516e-05 +6.437793731422292501e-05 +6.437793731422292501e-05 +4.927754965352772610e-05 +4.927754965352772610e-05 +3.301491077152221852e-05 +3.301491077152221852e-05 +1.936419606626931577e-05 +1.936419606626931577e-05 +9.942539684234432941e-06 +9.942539684234432941e-06 +4.468662792935623921e-06 +4.468662792935623921e-06 +1.758035637746577950e-06 +1.758035637746577950e-06 +6.053950015225363516e-07 +6.053950015225363516e-07 +1.827435046263828092e-07 +1.827435046263828092e-07 diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/weno5/CMakeLists.txt b/tests_cpp/eigen_2d_burgers_outflow_implicit/weno5/CMakeLists.txt new file mode 100644 index 00000000..9a5716c0 --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/weno5/CMakeLists.txt @@ -0,0 +1,19 @@ + +set(testname eigen_2d_burgers_outflow_weno5_implicit) +set(exename ${testname}_exe) + +configure_file(../compare.py compare.py COPYONLY) +configure_file(gold.txt gold.txt COPYONLY) +configure_file(../plot.py plot.py COPYONLY) + +add_executable(${exename} ${CMAKE_CURRENT_SOURCE_DIR}/../main.cc) +target_compile_definitions(${exename} PUBLIC -DUSE_WENO5) + +add_test(NAME ${testname} +COMMAND ${CMAKE_COMMAND} +-DMESHDRIVER=${MESHSRC}/create_full_mesh.py +-DOUTDIR=${CMAKE_CURRENT_BINARY_DIR} +-DEXENAME=$ +-DSTENCILVAL=7 +-P ${CMAKE_CURRENT_SOURCE_DIR}/../test.cmake +) diff --git a/tests_cpp/eigen_2d_burgers_outflow_implicit/weno5/gold.txt b/tests_cpp/eigen_2d_burgers_outflow_implicit/weno5/gold.txt new file mode 100644 index 00000000..3de59eec --- /dev/null +++ b/tests_cpp/eigen_2d_burgers_outflow_implicit/weno5/gold.txt @@ -0,0 +1,800 @@ +2.872555083269606567e-05 +2.872555083269606567e-05 +9.505115307558087900e-05 +9.505115307558087900e-05 +2.737642781910634639e-04 +2.737642781910634639e-04 +6.837668209684796524e-04 +6.837668209684796524e-04 +1.473508491452706114e-03 +1.473508491452706114e-03 +2.736899927773278011e-03 +2.736899927773278011e-03 +4.408370636399741929e-03 +4.408370636399741929e-03 +6.225264348165724868e-03 +6.225264348165724868e-03 +7.811161536498911762e-03 +7.811161536498911762e-03 +8.806162388682703401e-03 +8.806162388682703401e-03 +8.966227232250816026e-03 +8.966227232250816026e-03 +8.206641259407938468e-03 +8.206641259407938468e-03 +6.674324709850210438e-03 +6.674324709850210438e-03 +4.746886933102251356e-03 +4.746886933102252223e-03 +2.914379696683458247e-03 +2.914379696683458247e-03 +1.538689839356357624e-03 +1.538689839356357624e-03 +7.008477900491153151e-04 +7.008477900491153151e-04 +2.770081926227910643e-04 +2.770081926227910643e-04 +9.550145182502520810e-05 +9.550145182502520810e-05 +2.882952016266132430e-05 +2.882952016266132430e-05 +7.291539594699776605e-05 +7.291539594699776605e-05 +2.404216705044922003e-04 +2.404216705044922003e-04 +6.865876080061136680e-04 +6.865876080061136680e-04 +1.687962981378506120e-03 +1.687962981378506120e-03 +3.552082567853952471e-03 +3.552082567853952471e-03 +6.415064511833978682e-03 +6.415064511833978682e-03 +1.007004310182199602e-02 +1.007004310182199602e-02 +1.398956842399657653e-02 +1.398956842399657653e-02 +1.749359236389133251e-02 +1.749359236389133251e-02 +1.990219199734364133e-02 +1.990219199734364133e-02 +2.068696049310670901e-02 +2.068696049310670901e-02 +1.948321198774304736e-02 +1.948321198774304736e-02 +1.635124655045732539e-02 +1.635124655045732192e-02 +1.195803159937665686e-02 +1.195803159937665686e-02 +7.464228093995859542e-03 +7.464228093995859542e-03 +3.953459616124052264e-03 +3.953459616124052264e-03 +1.794919274403298549e-03 +1.794919274403298549e-03 +7.072455585601387972e-04 +7.072455585601387972e-04 +2.433398996651103092e-04 +2.433398996651103092e-04 +7.335911591675769224e-05 +7.335911591675769224e-05 +1.615408936559489122e-04 +1.615408936559489122e-04 +5.290583392372023505e-04 +5.290583392372023505e-04 +1.486204959082870498e-03 +1.486204959082870498e-03 +3.541750875007921347e-03 +3.541750875007921347e-03 +7.129550458351711743e-03 +7.129550458351711743e-03 +1.226022654941474541e-02 +1.226022654941474541e-02 +1.842448734906254951e-02 +1.842448734906254951e-02 +2.481164114701116252e-02 +2.481164114701116252e-02 +3.055525468494655131e-02 +3.055525468494655131e-02 +3.480387358944584536e-02 +3.480387358944584536e-02 +3.682734971967805554e-02 +3.682734971967805554e-02 +3.595587819511916750e-02 +3.595587819511916056e-02 +3.175007994279734763e-02 +3.175007994279734763e-02 +2.462811694471289670e-02 +2.462811694471289323e-02 +1.619311778476283450e-02 +1.619311778476283450e-02 +8.824823771124550367e-03 +8.824823771124550367e-03 +4.030691208769150713e-03 +4.030691208769150713e-03 +1.584460935055327429e-03 +1.584460935055327429e-03 +5.434715869394736009e-04 +5.434715869394736009e-04 +1.633500208839404235e-04 +1.633500208839404235e-04 +3.124231930288363370e-04 +3.124231930288363370e-04 +1.014229417011464377e-03 +1.014229417011464377e-03 +2.782790587058836702e-03 +2.782790587058836702e-03 +6.385676632275385546e-03 +6.385676632275385546e-03 +1.226841351863383264e-02 +1.226841351863383264e-02 +2.015517576216856524e-02 +2.015517576216856524e-02 +2.919575138534359973e-02 +2.919575138534359973e-02 +3.836352187511051642e-02 +3.836352187511051642e-02 +4.670414872968226688e-02 +4.670414872968225994e-02 +5.330163686131859230e-02 +5.330163686131859230e-02 +5.730880536242663470e-02 +5.730880536242660694e-02 +5.786056151044519785e-02 +5.786056151044519785e-02 +5.392778562613095850e-02 +5.392778562613095156e-02 +4.482119583310888861e-02 +4.482119583310887473e-02 +3.154595337816208656e-02 +3.154595337816208656e-02 +1.784952796324977206e-02 +1.784952796324977206e-02 +8.116030766075957650e-03 +8.116030766075957650e-03 +3.141706431651832086e-03 +3.141706431651832086e-03 +1.067587814296756376e-03 +1.067587814296756376e-03 +3.188561142043909100e-04 +3.188561142043909100e-04 +5.277218241356750527e-04 +5.277218241356750527e-04 +1.696153606730550507e-03 +1.696153606730550507e-03 +4.528570305369526532e-03 +4.528570305369526532e-03 +1.000498311277565021e-02 +1.000498311277565021e-02 +1.843760703748644619e-02 +1.843760703748644619e-02 +2.920441989991957171e-02 +2.920441989991957865e-02 +4.119031129426077575e-02 +4.119031129426077575e-02 +5.324584904006132019e-02 +5.324584904006134795e-02 +6.439469413751547355e-02 +6.439469413751547355e-02 +7.372264801176460491e-02 +7.372264801176457716e-02 +8.034606923716346161e-02 +8.034606923716346161e-02 +8.328222769155511906e-02 +8.328222769155511906e-02 +8.131228179763733133e-02 +8.131228179763734520e-02 +7.252673975919805838e-02 +7.252673975919803062e-02 +5.544457121022797197e-02 +5.544457121022797197e-02 +3.324902942402468281e-02 +3.324902942402468975e-02 +1.493549989980793094e-02 +1.493549989980793094e-02 +5.551398601227701228e-03 +5.551398601227701228e-03 +1.848715578244682512e-03 +1.848715578244682512e-03 +5.459317698005885502e-04 +5.459317698005885502e-04 +7.798587395375970949e-04 +7.798587395375970949e-04 +2.483694405156190114e-03 +2.483694405156190114e-03 +6.472119596743276457e-03 +6.472119596743276457e-03 +1.387627712150927423e-02 +1.387627712150927423e-02 +2.482185544261909793e-02 +2.482185544261909793e-02 +3.838191996981184934e-02 +3.838191996981184934e-02 +5.324942241977104890e-02 +5.324942241977104890e-02 +6.819143565729945955e-02 +6.819143565729947343e-02 +8.222774606242612416e-02 +8.222774606242615192e-02 +9.448104949937510821e-02 +9.448104949937503882e-02 +1.041027923342827260e-01 +1.041027923342826705e-01 +1.101072418407885306e-01 +1.101072418407884890e-01 +1.112021368713179137e-01 +1.112021368713178998e-01 +1.054410871890366641e-01 +1.054410871890367057e-01 +8.851242059599379042e-02 +8.851242059599381817e-02 +5.795706578794054820e-02 +5.795706578794054820e-02 +2.588771508951703479e-02 +2.588771508951703479e-02 +8.807368670131928562e-03 +8.807368670131928562e-03 +2.825682360804956158e-03 +2.825682360804956158e-03 +8.199376101029054453e-04 +8.199376101029054453e-04 +1.010781603622840823e-03 +1.010781603622840823e-03 +3.201528214919670962e-03 +3.201528214919670962e-03 +8.216844828091313460e-03 +8.216844828091313460e-03 +1.732990021939252723e-02 +1.732990021939252723e-02 +3.055135006264700359e-02 +3.055135006264700359e-02 +4.673397670130106235e-02 +4.673397670130104847e-02 +6.440414548044542942e-02 +6.440414548044542942e-02 +8.222773527488005940e-02 +8.222773527488005940e-02 +9.919418254012495051e-02 +9.919418254012496439e-02 +1.144620413978425233e-01 +1.144620413978424817e-01 +1.272367102730990462e-01 +1.272367102730989907e-01 +1.366592838948322841e-01 +1.366592838948322841e-01 +1.414550125761958321e-01 +1.414550125761958599e-01 +1.398468980764651648e-01 +1.398468980764651648e-01 +1.282684871807271698e-01 +1.282684871807271421e-01 +9.493363543532806759e-02 +9.493363543532809534e-02 +4.381181251915330227e-02 +4.381181251915330921e-02 +1.266267884510490112e-02 +1.266267884510490112e-02 +3.804384817367174008e-03 +3.804384817367174008e-03 +1.079239925500768103e-03 +1.079239925500768103e-03 +1.151788756815191406e-03 +1.151788756815191406e-03 +3.649076801161624251e-03 +3.649076801161624251e-03 +9.343043434831506400e-03 +9.343043434831506400e-03 +1.969695573773456998e-02 +1.969695573773456998e-02 +3.477897414550479988e-02 +3.477897414550479988e-02 +5.334159901200899445e-02 +5.334159901200900139e-02 +7.374164792823792447e-02 +7.374164792823789671e-02 +9.448301153092594773e-02 +9.448301153092596161e-02 +1.144607814626998205e-01 +1.144607814626998760e-01 +1.328318448493940040e-01 +1.328318448493940318e-01 +1.488578910437412528e-01 +1.488578910437412528e-01 +1.617806827816759274e-01 +1.617806827816758997e-01 +1.705654906427084283e-01 +1.705654906427084005e-01 +1.735827065086567589e-01 +1.735827065086567866e-01 +1.681150460606422037e-01 +1.681150460606422592e-01 +1.441848757230303757e-01 +1.441848757230304312e-01 +7.440935418782608501e-02 +7.440935418782607114e-02 +1.691376571571403392e-02 +1.691376571571403392e-02 +4.483522018030907698e-03 +4.483522018030907698e-03 +1.242306277728105583e-03 +1.242306277728105583e-03 +1.155075910490895117e-03 +1.155075910490895117e-03 +3.682227897346315032e-03 +3.682227897346315032e-03 +9.539938258905924021e-03 +9.539938258905924021e-03 +2.045570225273929052e-02 +2.045570225273929052e-02 +3.678642748143022556e-02 +3.678642748143021862e-02 +5.735572078103282068e-02 +5.735572078103283455e-02 +8.037513799857867736e-02 +8.037513799857867736e-02 +1.041090043690522382e-01 +1.041090043690522382e-01 +1.272354419511117218e-01 +1.272354419511117218e-01 +1.488567077380867909e-01 +1.488567077380867909e-01 +1.682123678707450054e-01 +1.682123678707450332e-01 +1.846949522738964844e-01 +1.846949522738964844e-01 +1.975146050885957794e-01 +1.975146050885957238e-01 +2.052186682430302833e-01 +2.052186682430302278e-01 +2.057040271891829586e-01 +2.057040271891829586e-01 +1.937191168454417411e-01 +1.937191168454417689e-01 +1.272688177103888651e-01 +1.272688177103888374e-01 +2.410592957258068589e-02 +2.410592957258066854e-02 +4.599641659053638719e-03 +4.599641659053638719e-03 +1.247368427422860778e-03 +1.247368427422860778e-03 +1.018493102078286234e-03 +1.018493102078286234e-03 +3.279550995383244145e-03 +3.279550995383244145e-03 +8.692584150039374749e-03 +8.692584150039376484e-03 +1.925058132490464410e-02 +1.925058132490464063e-02 +3.591745109181514184e-02 +3.591745109181512796e-02 +5.791286916015237624e-02 +5.791286916015239011e-02 +8.331796145210679239e-02 +8.331796145210680626e-02 +1.101176173302625960e-01 +1.101176173302625544e-01 +1.366589912141640828e-01 +1.366589912141640828e-01 +1.617789724392347250e-01 +1.617789724392347250e-01 +1.846942315551728386e-01 +1.846942315551727831e-01 +2.048079617451222800e-01 +2.048079617451223078e-01 +2.214932967206048875e-01 +2.214932967206049153e-01 +2.338420969515596082e-01 +2.338420969515596359e-01 +2.403111037175892362e-01 +2.403111037175892639e-01 +2.373965678969489346e-01 +2.373965678969490178e-01 +1.958673144850377190e-01 +1.958673144850377190e-01 +4.479763661313836132e-02 +4.479763661313831968e-02 +4.182663719869376441e-03 +4.182663719869376441e-03 +1.090647937398000296e-03 +1.090647937398000296e-03 +7.876783513847759747e-04 +7.876783513847759747e-04 +2.564153132969016445e-03 +2.564153132969016445e-03 +6.994555811957059270e-03 +6.994555811957060137e-03 +1.615679115076007824e-02 +1.615679115076007824e-02 +3.173693140264026435e-02 +3.173693140264026435e-02 +5.398270695003794944e-02 +5.398270695003793557e-02 +8.134758807191531937e-02 +8.134758807191529162e-02 +1.112135391285231784e-01 +1.112135391285231228e-01 +1.414554169813740103e-01 +1.414554169813740658e-01 +1.705636024392850336e-01 +1.705636024392849781e-01 +1.975134331916958141e-01 +1.975134331916957586e-01 +2.214929569170240586e-01 +2.214929569170240864e-01 +2.420869127739254112e-01 +2.420869127739255222e-01 +2.591321770916855116e-01 +2.591321770916856781e-01 +2.704307128313725994e-01 +2.704307128313726549e-01 +2.741933448333433443e-01 +2.741933448333433998e-01 +2.568499327578668279e-01 +2.568499327578667168e-01 +9.850584575914009267e-02 +9.850584575914000940e-02 +3.932394493596183538e-03 +3.932394493596181803e-03 +8.311962497912776977e-04 +8.311962497912776977e-04 +5.328597623791420232e-04 +5.328597623791420232e-04 +1.750134951744371226e-03 +1.750134951744371226e-03 +4.907263094111299732e-03 +4.907263094111299732e-03 +1.183416831723596054e-02 +1.183416831723595880e-02 +2.464032121681640833e-02 +2.464032121681640833e-02 +4.486408671390747921e-02 +4.486408671390748615e-02 +7.254900346100696951e-02 +7.254900346100696951e-02 +1.054478679766742422e-01 +1.054478679766742283e-01 +1.398467127207250749e-01 +1.398467127207250749e-01 +1.735807744393980945e-01 +1.735807744393980667e-01 +2.052172642235866040e-01 +2.052172642235866040e-01 +2.338415435484728078e-01 +2.338415435484728633e-01 +2.591320971332461420e-01 +2.591320971332461975e-01 +2.798645725826297559e-01 +2.798645725826298669e-01 +2.965617104930278192e-01 +2.965617104930279302e-01 +3.090226318968884445e-01 +3.090226318968883890e-01 +3.085424142381226198e-01 +3.085424142381225088e-01 +1.735338286015592202e-01 +1.735338286015590814e-01 +7.227259753144696548e-03 +7.227259753144685273e-03 +5.554783040650471040e-04 +5.554783040650471040e-04 +3.147970531959626240e-04 +3.147970531959626240e-04 +1.039506859183037370e-03 +1.039506859183037370e-03 +2.973763336067419914e-03 +2.973763336067420347e-03 +7.409491219088334100e-03 +7.409491219088334100e-03 +1.620931890789339275e-02 +1.620931890789339275e-02 +3.156489857502521590e-02 +3.156489857502522284e-02 +5.545028016263851833e-02 +5.545028016263851833e-02 +8.851215795612475679e-02 +8.851215795612478454e-02 +1.282668048574594766e-01 +1.282668048574594211e-01 +1.681131824024822696e-01 +1.681131824024822696e-01 +2.057026418858139416e-01 +2.057026418858139694e-01 +2.403104431523657436e-01 +2.403104431523657436e-01 +2.704305381265683095e-01 +2.704305381265683650e-01 +2.965617097303044902e-01 +2.965617097303046013e-01 +3.182896015308052662e-01 +3.182896015308052662e-01 +3.340831565257997848e-01 +3.340831565257997848e-01 +3.451934794710870835e-01 +3.451934794710869725e-01 +2.411783170948159638e-01 +2.411783170948158250e-01 +1.302185537050496063e-02 +1.302185537050494848e-02 +3.279045227196281174e-04 +3.279045227196281174e-04 +1.623368440866359403e-04 +1.623368440866359403e-04 +5.375841419514002005e-04 +5.375841419514002005e-04 +1.554104060753549923e-03 +1.554104060753549923e-03 +3.937992406960883797e-03 +3.937992406960883797e-03 +8.831953840527683938e-03 +8.831953840527683938e-03 +1.785300029290352516e-02 +1.785300029290352863e-02 +3.324891059345760208e-02 +3.324891059345760208e-02 +5.795645706712884010e-02 +5.795645706712885398e-02 +9.493348425490935605e-02 +9.493348425490932829e-02 +1.441848021593984674e-01 +1.441848021593984952e-01 +1.937187595136029006e-01 +1.937187595136028728e-01 +2.373962109950677746e-01 +2.373962109950677191e-01 +2.741932317658227758e-01 +2.741932317658227758e-01 +3.090226761265267852e-01 +3.090226761265268407e-01 +3.340832370287019626e-01 +3.340832370287019626e-01 +3.558459902664787844e-01 +3.558459902664787289e-01 +3.832034566260511532e-01 +3.832034566260510977e-01 +2.427982241635683602e-01 +2.427982241635683047e-01 +6.799226343597920309e-03 +6.799226343597916840e-03 +1.653210747998872432e-04 +1.653210747998872432e-04 +7.310859571484665210e-05 +7.310859571484665210e-05 +2.425170221551450781e-04 +2.425170221551450781e-04 +7.036989477970672811e-04 +7.036989477970672811e-04 +1.792356007629813605e-03 +1.792356007629813605e-03 +4.031983415764348565e-03 +4.031983415764348565e-03 +8.116159578801275531e-03 +8.116159578801275531e-03 +1.493512396753889429e-02 +1.493512396753889429e-02 +2.588662828108775557e-02 +2.588662828108775557e-02 +4.380896418098268602e-02 +4.380896418098268602e-02 +7.440553375438524464e-02 +7.440553375438524464e-02 +1.272659953386175435e-01 +1.272659953386176268e-01 +1.958664878124799502e-01 +1.958664878124798947e-01 +2.568497776471984628e-01 +2.568497776471984628e-01 +3.085424659321741592e-01 +3.085424659321741592e-01 +3.451935290658980082e-01 +3.451935290658980637e-01 +3.832027486913189285e-01 +3.832027486913188730e-01 +3.785081183925155246e-01 +3.785081183925154691e-01 +6.349274648644585850e-02 +6.349274648644585850e-02 +3.262993748132053760e-04 +3.262993748132053760e-04 +7.355575545791597403e-05 +7.355575545791597403e-05 +2.878417633875412273e-05 +2.878417633875412273e-05 +9.554248810953490172e-05 +9.554248810953490172e-05 +2.775169765631349409e-04 +2.775169765631349409e-04 +7.069969869371577498e-04 +7.069969869371577498e-04 +1.584610677719893845e-03 +1.584610677719893845e-03 +3.142121150415842432e-03 +3.142121150415842865e-03 +5.553101314729870161e-03 +5.553101314729870161e-03 +8.812175549613223205e-03 +8.812175549613223205e-03 +1.267269533704302661e-02 +1.267269533704302661e-02 +1.692810915879333980e-02 +1.692810915879333980e-02 +2.411867795861706576e-02 +2.411867795861706576e-02 +4.480371947703632590e-02 +4.480371947703630509e-02 +9.850675655352229831e-02 +9.850675655352229831e-02 +1.735338951668557450e-01 +1.735338951668556617e-01 +2.411789135273784968e-01 +2.411789135273784412e-01 +2.428169114374538751e-01 +2.428169114374538196e-01 +6.349055477697683469e-02 +6.349055477697680694e-02 +4.199002413641269422e-04 +4.199002413641268880e-04 +9.600310355639513039e-05 +9.600310355639513039e-05 +2.888823462211699076e-05 +2.888823462211699076e-05 +9.914690169635016724e-06 +9.914690169635016724e-06 +3.291269729790618955e-05 +3.291269729790618955e-05 +9.560149087867292317e-05 +9.560149087867292317e-05 +2.433022618710791655e-04 +2.433022618710791655e-04 +5.433631019703035608e-04 +5.433631019703035608e-04 +1.066864523189934568e-03 +1.066864523189934568e-03 +1.846075342423255739e-03 +1.846075342423255739e-03 +2.819088331979050129e-03 +2.819088331979050129e-03 +3.791975633204042965e-03 +3.791975633204042965e-03 +4.466137412513410815e-03 +4.466137412513409947e-03 +4.582548556093943987e-03 +4.582548556093943987e-03 +4.171467118629923476e-03 +4.171467118629923476e-03 +3.927088741292592868e-03 +3.927088741292592868e-03 +7.224825638779339143e-03 +7.224825638779332204e-03 +1.301542021939238770e-02 +1.301542021939237209e-02 +6.788302836763500901e-03 +6.788302836763498299e-03 +3.245016950693124249e-04 +3.245016950693124249e-04 +9.600205645207307125e-05 +9.600205645207307125e-05 +3.296714108898086905e-05 +3.296714108898086905e-05 +9.940075302262201066e-06 +9.940075302262201066e-06 +2.988664749935617133e-06 +2.988664749935617133e-06 +9.920924570081776942e-06 +9.920924570081776942e-06 +2.881317967875166792e-05 +2.881317967875166792e-05 +7.326904925447470470e-05 +7.326904925447470470e-05 +1.632250062373032053e-04 +1.632250062373032053e-04 +3.188906369253214084e-04 +3.188906369253214084e-04 +5.465710096525206807e-04 +5.465710096525206807e-04 +8.217080390470497602e-04 +8.217080390470497602e-04 +1.082455036809008711e-03 +1.082455036809008711e-03 +1.246557288649524772e-03 +1.246557288649524772e-03 +1.251481481860872201e-03 +1.251481481860872201e-03 +1.093564877026116872e-03 +1.093564877026116872e-03 +8.326500617539090527e-04 +8.326500617539090527e-04 +5.557655790220897719e-04 +5.557655790220897719e-04 +3.284118741657287600e-04 +3.284118741657287600e-04 +1.653473837208407882e-04 +1.653473837208407882e-04 +7.346977799051132585e-05 +7.346977799051132585e-05 +2.884874749537819899e-05 +2.884874749537819899e-05 +9.925858485812810467e-06 +9.925858485812810467e-06 +2.995136897139399850e-06 +2.995136897139399850e-06 +7.885053095920514014e-07 +7.885053095920514014e-07 +2.617389777183280353e-06 +2.617389777183280353e-06 +7.600937626803264477e-06 +7.600937626803264477e-06 +1.932222232485061964e-05 +1.932222232485061964e-05 +4.300068862743507863e-05 +4.300068862743507863e-05 +8.379561671722442803e-05 +8.379561671722442803e-05 +1.430549121404157324e-04 +1.430549121404157324e-04 +2.139464485629094059e-04 +2.139464485629094059e-04 +2.801103099579138004e-04 +2.801103099579138004e-04 +3.207747091802685621e-04 +3.207747091802685621e-04 +3.210589638091714037e-04 +3.210589638091714037e-04 +2.807612123041981246e-04 +2.807612123041981246e-04 +2.145820104749888502e-04 +2.145820104749888502e-04 +1.434580436888381453e-04 +1.434580436888381453e-04 +8.397143058324390489e-05 +8.397143058324390489e-05 +4.305480403940289269e-05 +4.305480403940289269e-05 +1.933576570876879773e-05 +1.933576570876879773e-05 +7.603459274283963977e-06 +7.603459274283963977e-06 +2.617733981131449207e-06 +2.617733981131449207e-06 +7.901141519359049630e-07 +7.901141519359049630e-07 +1.823784502805593729e-07 +1.823784502805593729e-07 +6.053837776196951026e-07 +6.053837776196951026e-07 +1.757960088509884190e-06 +1.757960088509884190e-06 +4.468283894874134908e-06 +4.468283894874134908e-06 +9.941094234632133530e-06 +9.941094234632133530e-06 +1.935925918764563875e-05 +1.935925918764563875e-05 +3.299912572322029374e-05 +3.299912572322029374e-05 +4.923826898070024017e-05 +4.923826898070024017e-05 +6.431094444417053335e-05 +6.431094444417053335e-05 +7.351330668568372426e-05 +7.351330668568372426e-05 +7.352753103316512964e-05 +7.352753103316512964e-05 +6.434357041419475033e-05 +6.434357041419475033e-05 +4.927003838266346287e-05 +4.927003838266346287e-05 +3.301895763488334651e-05 +3.301895763488334651e-05 +1.936816880633240542e-05 +1.936816880633240542e-05 +9.944162534031167691e-06 +9.944162534031167691e-06 +4.469028507916289829e-06 +4.469028507916289829e-06 +1.758094146113675732e-06 +1.758094146113675732e-06 +6.054025130129001712e-07 +6.054025130129001712e-07 +1.827442643655394671e-07 +1.827442643655394671e-07