diff --git a/include/pressiodemoapps/gradient.hpp b/include/pressiodemoapps/gradient.hpp new file mode 100644 index 00000000..4ce6e134 --- /dev/null +++ b/include/pressiodemoapps/gradient.hpp @@ -0,0 +1,84 @@ +/* +//@HEADER +// ************************************************************************ +// +// mesh.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_GRADIENT_PUBLIC_HPP_ +#define PRESSIODEMOAPPS_GRADIENT_PUBLIC_HPP_ + +#include "mesh.hpp" +#include "schemes_info.hpp" +#include "impl/gradient_2d.hpp" + +namespace pressiodemoapps{ + +constexpr BoundaryFacesGradientScheme defaultGradSchemeForBoundaryFaces = + BoundaryFacesGradientScheme::OneSidedFdAutoStencil; + +template +class GradientEvaluator +{ +public: + explicit GradientEvaluator(const MeshType & mesh) + : m_impl(mesh, defaultGradSchemeForBoundaryFaces) + {} + + template + void compute(const StateType & state){ + m_impl.compute(state); + } + + const auto & queryFace(int cellGID, FacePosition fp) const{ + return m_impl.queryFace(cellGID, fp); + } + +private: + using face_t = impl::Face; + impl::GradientInternal m_impl; +}; + +}//end namespace pressiodemoapps +#endif diff --git a/include/pressiodemoapps/impl/gradient_2d.hpp b/include/pressiodemoapps/impl/gradient_2d.hpp new file mode 100644 index 00000000..1de4f5fa --- /dev/null +++ b/include/pressiodemoapps/impl/gradient_2d.hpp @@ -0,0 +1,276 @@ +/* +//@HEADER +// ************************************************************************ +// +// functor_fill_stencil.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_IMPL_GRADIENTS_2D_HPP_ +#define PRESSIODEMOAPPS_IMPL_GRADIENTS_2D_HPP_ + +namespace pressiodemoapps{ namespace impl{ + +template +ScalarType fd_face_normal_gradient_for_cell_centered_state_2d(const StateType & f, + const CellConnectivityVectorLike & cc, + char axis, + GradFdMode mode, + const ScalarType & h, + int nDofPerCell = 1, + int dofShift = 0) +{ + + assert(axis == 'x' || axis == 'y'); + const int graphNumCols = cc.size(); + assert(graphNumCols ==5 || graphNumCols==9); + + const int i_05 = cc[0]; // this always exists + const int i_p15 = (axis=='x') ? cc[3] : (axis=='y') ? cc[2] : -1; + const int i_m15 = (axis=='x') ? cc[1] : (axis=='y') ? cc[4] : -1; + const int i_p30 = (graphNumCols==9) ? ((axis=='x') ? cc[7] : ((axis=='y') ? cc[6] : -1)) : -1; + const int i_m30 = (graphNumCols==9) ? ((axis=='x') ? cc[5] : ((axis=='y') ? cc[8] : -1)) : -1; + + const auto f_05 = f(i_05*nDofPerCell + dofShift); + const auto f_p15 = (i_p15 != -1) ? f(i_p15*nDofPerCell + dofShift) : 0; + const auto f_m15 = (i_m15 != -1) ? f(i_m15*nDofPerCell + dofShift) : 0; + const auto f_p30 = (i_p30 != -1) ? f(i_p30*nDofPerCell + dofShift) : 0; + const auto f_m30 = (i_m30 != -1) ? f(i_m30*nDofPerCell + dofShift) : 0; + + // https://web.media.mit.edu/~crtaylor/calculator.html + switch(mode){ + case GradFdMode::ForwardTwoPt: return (-f_05 + f_p15)/h; + case GradFdMode::BackwardTwoPt: return ( f_05 - f_m15)/h; + case GradFdMode::CenterThreePt: return (-f_m15 + f_p15)/(2*h); + + case GradFdMode::ForwardThreePt: return (-2.*f_05 + 3.*f_p15 - 1.*f_p30)/h; + case GradFdMode::BackwardThreePt: return ( 2.*f_05 - 3.*f_m15 + 1.*f_m30)/h; + case GradFdMode::CenterFivePt: return (-f_m30 - 8.*f_m15 + 8.*f_p15 - f_p30)/(12.*h); + + default: + return 0; + } +} + +template +struct Face{ + static constexpr std::size_t N_ = N; + std::array centerCoords = {}; + std::array normalGradValue = {}; +}; + +struct MyKey{ + int parentCellGID; + FacePosition pos; + + bool operator==(const MyKey & p) const { + return parentCellGID == p.parentCellGID && pos == p.pos; + } +}; + +struct hash_fn{ + std::size_t operator() (const MyKey & k) const{ + const std::size_t h1 = std::hash()(k.parentCellGID); + const std::size_t h2 = std::hash()(static_cast(k.pos)); + return h1 ^ h2; + } +}; + +template< + class MeshType, + class FaceType = Face > +class GradientInternal{ + using sc_t = typename MeshType::scalar_type; + +public: + GradientInternal(const MeshType & mesh, + BoundaryFacesGradientScheme schemeAtBoundaryFaces) + : m_meshObj(mesh), + m_schemeAtBoundaryFaces(schemeAtBoundaryFaces) + { + initialize(mesh); + } + + template + void compute(const StateType & state){ + if (m_schemeAtBoundaryFaces == BoundaryFacesGradientScheme::OneSidedFdAutoStencil){ + this->computeOneSidedFdAutoStencil(state); + } + else{ + throw std::runtime_error("invalid choice"); + } + } + + const auto & queryFace(int cellGID, FacePosition fp) const{ + const MyKey key{cellGID, fp}; + std::cout << m_data.count(key) << '\n'; + assert(m_data.count(key) == 1); + auto it = m_data.find(key); + return it->second; + } + +private: + template + void computeOneSidedFdAutoStencil(const StateType & state) + { + constexpr std::size_t nDofPerCell = FaceType::N_; + const int ss = m_meshObj.get().stencilSize(); + const auto SkewRight = (ss == 3) ? GradFdMode::ForwardTwoPt + : GradFdMode::ForwardThreePt; + const auto SkewLeft = (ss == 3) ? GradFdMode::BackwardTwoPt + : GradFdMode::BackwardThreePt; + + const auto dx = m_meshObj.get().dx(); + const auto dy = m_meshObj.get().dy(); + + const auto & rowsForLoop = m_meshObj.get().graphRowsOfCellsStrictlyOnBd(); + const auto & G = m_meshObj.get().graph(); + for (auto rowInd : rowsForLoop){ + const bool bdL = m_meshObj.get().cellHasLeftFaceOnBoundary2d(rowInd); + const bool bdF = m_meshObj.get().cellHasFrontFaceOnBoundary2d(rowInd); + const bool bdR = m_meshObj.get().cellHasRightFaceOnBoundary2d(rowInd); + const bool bdB = m_meshObj.get().cellHasBackFaceOnBoundary2d(rowInd); + //std::cout << "row = " << rowInd << " " << G(rowInd,0) << std::endl; + + const auto currCellConnec = G.row(rowInd); + const int cellGID = currCellConnec(0); + + if (bdL){ + MyKey key{cellGID, FacePosition::Left}; + assert(m_data.count(key) == 1); + auto & faceObj= m_data[key]; + for (int j=0; j(1)/static_cast(2); + const auto dxHalf = mesh.dx()*half; + const auto dyHalf = mesh.dy()*half; + + for (auto rowInd : rowsForLoop){ + const bool bdL = mesh.cellHasLeftFaceOnBoundary2d(rowInd); + const bool bdF = mesh.cellHasFrontFaceOnBoundary2d(rowInd); + const bool bdR = mesh.cellHasRightFaceOnBoundary2d(rowInd); + const bool bdB = mesh.cellHasBackFaceOnBoundary2d(rowInd); + //std::cout << "row = " << rowInd << " " << G(rowInd,0) << std::endl; + + const auto graphRow = G.row(rowInd); + const int cellGID = graphRow(0); + const auto cellX = x(cellGID); + const auto cellY = y(cellGID); + const auto cellZ = z(cellGID); + + if (bdL) { + m_data[MyKey{cellGID, FacePosition::Left}] = FaceType{{cellX-dxHalf, cellY, cellZ}}; + } + + // if (bdF) { + // m_data[MyKey{cellGID, FacePosition::Front}] = FaceType{{cellX, cellY+dyHalf, cellZ}}; + // } + + if (bdR) { + m_data[MyKey{cellGID, FacePosition::Right}] = FaceType{{cellX+dxHalf, cellY, cellZ}}; + } + // if (bdB) { + // m_data[MyKey{cellGID, FacePosition::Back}] = FaceType{{cellX, cellY-dyHalf, cellZ}}; + // } + } + } + +private: + BoundaryFacesGradientScheme m_schemeAtBoundaryFaces; + std::reference_wrapper m_meshObj; + std::unordered_map m_data = {}; +}; + +}} +#endif diff --git a/include/pressiodemoapps/schemes_info.hpp b/include/pressiodemoapps/schemes_info.hpp index 44cdd7d1..3363ead2 100644 --- a/include/pressiodemoapps/schemes_info.hpp +++ b/include/pressiodemoapps/schemes_info.hpp @@ -109,6 +109,20 @@ bool stencilSizeCompatibleWithViscousFluxReconstruction(ViscousFluxReconstructio return false; } +enum class FacePosition{ + Left, Front, Right, Back, Bottom, Top +}; + +enum class GradFdMode{ + ForwardTwoPt, BackwardTwoPt, CenterThreePt, + ForwardThreePt, BackwardThreePt, CenterFivePt +}; + +enum class BoundaryFacesGradientScheme{ + OneSidedFdAutoStencil, // this selects the stencil cells automatically using the connectivity + LSQAutoStencil +}; + }//end namespace pressiodemoapps #include "./weno.hpp" diff --git a/tests_cpp/gradients/connectivity.dat b/tests_cpp/gradients/connectivity.dat index c73879ab..a091cc28 100644 --- a/tests_cpp/gradients/connectivity.dat +++ b/tests_cpp/gradients/connectivity.dat @@ -1,120 +1,120 @@ - 0 -1 10 1 -1 -1 20 2 -1 - 1 0 11 2 -1 -1 21 3 -1 - 2 1 12 3 -1 0 22 4 -1 - 3 2 13 4 -1 1 23 5 -1 - 4 3 14 5 -1 2 24 6 -1 - 5 4 15 6 -1 3 25 7 -1 - 6 5 16 7 -1 4 26 8 -1 - 7 6 17 8 -1 5 27 9 -1 - 8 7 18 9 -1 6 28 -1 -1 - 9 8 19 -1 -1 7 29 -1 -1 - 10 -1 20 11 0 -1 30 12 -1 - 11 10 21 12 1 -1 31 13 -1 - 12 11 22 13 2 10 32 14 -1 - 13 12 23 14 3 11 33 15 -1 - 14 13 24 15 4 12 34 16 -1 - 15 14 25 16 5 13 35 17 -1 - 16 15 26 17 6 14 36 18 -1 - 17 16 27 18 7 15 37 19 -1 - 18 17 28 19 8 16 38 -1 -1 - 19 18 29 -1 9 17 39 -1 -1 - 20 -1 30 21 10 -1 40 22 0 - 21 20 31 22 11 -1 41 23 1 - 22 21 32 23 12 20 42 24 2 - 23 22 33 24 13 21 43 25 3 - 24 23 34 25 14 22 44 26 4 - 25 24 35 26 15 23 45 27 5 - 26 25 36 27 16 24 46 28 6 - 27 26 37 28 17 25 47 29 7 - 28 27 38 29 18 26 48 -1 8 - 29 28 39 -1 19 27 49 -1 9 - 30 -1 40 31 20 -1 50 32 10 - 31 30 41 32 21 -1 51 33 11 - 32 31 42 33 22 30 52 34 12 - 33 32 43 34 23 31 53 35 13 - 34 33 44 35 24 32 54 36 14 - 35 34 45 36 25 33 55 37 15 - 36 35 46 37 26 34 56 38 16 - 37 36 47 38 27 35 57 39 17 - 38 37 48 39 28 36 58 -1 18 - 39 38 49 -1 29 37 59 -1 19 - 40 -1 50 41 30 -1 60 42 20 - 41 40 51 42 31 -1 61 43 21 - 42 41 52 43 32 40 62 44 22 - 43 42 53 44 33 41 63 45 23 - 44 43 54 45 34 42 64 46 24 - 45 44 55 46 35 43 65 47 25 - 46 45 56 47 36 44 66 48 26 - 47 46 57 48 37 45 67 49 27 - 48 47 58 49 38 46 68 -1 28 - 49 48 59 -1 39 47 69 -1 29 - 50 -1 60 51 40 -1 70 52 30 - 51 50 61 52 41 -1 71 53 31 - 52 51 62 53 42 50 72 54 32 - 53 52 63 54 43 51 73 55 33 - 54 53 64 55 44 52 74 56 34 - 55 54 65 56 45 53 75 57 35 - 56 55 66 57 46 54 76 58 36 - 57 56 67 58 47 55 77 59 37 - 58 57 68 59 48 56 78 -1 38 - 59 58 69 -1 49 57 79 -1 39 - 60 -1 70 61 50 -1 80 62 40 - 61 60 71 62 51 -1 81 63 41 - 62 61 72 63 52 60 82 64 42 - 63 62 73 64 53 61 83 65 43 - 64 63 74 65 54 62 84 66 44 - 65 64 75 66 55 63 85 67 45 - 66 65 76 67 56 64 86 68 46 - 67 66 77 68 57 65 87 69 47 - 68 67 78 69 58 66 88 -1 48 - 69 68 79 -1 59 67 89 -1 49 - 70 -1 80 71 60 -1 90 72 50 - 71 70 81 72 61 -1 91 73 51 - 72 71 82 73 62 70 92 74 52 - 73 72 83 74 63 71 93 75 53 - 74 73 84 75 64 72 94 76 54 - 75 74 85 76 65 73 95 77 55 - 76 75 86 77 66 74 96 78 56 - 77 76 87 78 67 75 97 79 57 - 78 77 88 79 68 76 98 -1 58 - 79 78 89 -1 69 77 99 -1 59 - 80 -1 90 81 70 -1 100 82 60 - 81 80 91 82 71 -1 101 83 61 - 82 81 92 83 72 80 102 84 62 - 83 82 93 84 73 81 103 85 63 - 84 83 94 85 74 82 104 86 64 - 85 84 95 86 75 83 105 87 65 - 86 85 96 87 76 84 106 88 66 - 87 86 97 88 77 85 107 89 67 - 88 87 98 89 78 86 108 -1 68 - 89 88 99 -1 79 87 109 -1 69 - 90 -1 100 91 80 -1 110 92 70 - 91 90 101 92 81 -1 111 93 71 - 92 91 102 93 82 90 112 94 72 - 93 92 103 94 83 91 113 95 73 - 94 93 104 95 84 92 114 96 74 - 95 94 105 96 85 93 115 97 75 - 96 95 106 97 86 94 116 98 76 - 97 96 107 98 87 95 117 99 77 - 98 97 108 99 88 96 118 -1 78 - 99 98 109 -1 89 97 119 -1 79 - 100 -1 110 101 90 -1 -1 102 80 - 101 100 111 102 91 -1 -1 103 81 - 102 101 112 103 92 100 -1 104 82 - 103 102 113 104 93 101 -1 105 83 - 104 103 114 105 94 102 -1 106 84 - 105 104 115 106 95 103 -1 107 85 - 106 105 116 107 96 104 -1 108 86 - 107 106 117 108 97 105 -1 109 87 - 108 107 118 109 98 106 -1 -1 88 - 109 108 119 -1 99 107 -1 -1 89 - 110 -1 -1 111 100 -1 -1 112 90 - 111 110 -1 112 101 -1 -1 113 91 - 112 111 -1 113 102 110 -1 114 92 - 113 112 -1 114 103 111 -1 115 93 - 114 113 -1 115 104 112 -1 116 94 - 115 114 -1 116 105 113 -1 117 95 - 116 115 -1 117 106 114 -1 118 96 - 117 116 -1 118 107 115 -1 119 97 - 118 117 -1 119 108 116 -1 -1 98 - 119 118 -1 -1 109 117 -1 -1 99 + 0 -1 10 1 -1 + 1 0 11 2 -1 + 2 1 12 3 -1 + 3 2 13 4 -1 + 4 3 14 5 -1 + 5 4 15 6 -1 + 6 5 16 7 -1 + 7 6 17 8 -1 + 8 7 18 9 -1 + 9 8 19 -1 -1 + 10 -1 20 11 0 + 11 10 21 12 1 + 12 11 22 13 2 + 13 12 23 14 3 + 14 13 24 15 4 + 15 14 25 16 5 + 16 15 26 17 6 + 17 16 27 18 7 + 18 17 28 19 8 + 19 18 29 -1 9 + 20 -1 30 21 10 + 21 20 31 22 11 + 22 21 32 23 12 + 23 22 33 24 13 + 24 23 34 25 14 + 25 24 35 26 15 + 26 25 36 27 16 + 27 26 37 28 17 + 28 27 38 29 18 + 29 28 39 -1 19 + 30 -1 40 31 20 + 31 30 41 32 21 + 32 31 42 33 22 + 33 32 43 34 23 + 34 33 44 35 24 + 35 34 45 36 25 + 36 35 46 37 26 + 37 36 47 38 27 + 38 37 48 39 28 + 39 38 49 -1 29 + 40 -1 50 41 30 + 41 40 51 42 31 + 42 41 52 43 32 + 43 42 53 44 33 + 44 43 54 45 34 + 45 44 55 46 35 + 46 45 56 47 36 + 47 46 57 48 37 + 48 47 58 49 38 + 49 48 59 -1 39 + 50 -1 60 51 40 + 51 50 61 52 41 + 52 51 62 53 42 + 53 52 63 54 43 + 54 53 64 55 44 + 55 54 65 56 45 + 56 55 66 57 46 + 57 56 67 58 47 + 58 57 68 59 48 + 59 58 69 -1 49 + 60 -1 70 61 50 + 61 60 71 62 51 + 62 61 72 63 52 + 63 62 73 64 53 + 64 63 74 65 54 + 65 64 75 66 55 + 66 65 76 67 56 + 67 66 77 68 57 + 68 67 78 69 58 + 69 68 79 -1 59 + 70 -1 80 71 60 + 71 70 81 72 61 + 72 71 82 73 62 + 73 72 83 74 63 + 74 73 84 75 64 + 75 74 85 76 65 + 76 75 86 77 66 + 77 76 87 78 67 + 78 77 88 79 68 + 79 78 89 -1 69 + 80 -1 90 81 70 + 81 80 91 82 71 + 82 81 92 83 72 + 83 82 93 84 73 + 84 83 94 85 74 + 85 84 95 86 75 + 86 85 96 87 76 + 87 86 97 88 77 + 88 87 98 89 78 + 89 88 99 -1 79 + 90 -1 100 91 80 + 91 90 101 92 81 + 92 91 102 93 82 + 93 92 103 94 83 + 94 93 104 95 84 + 95 94 105 96 85 + 96 95 106 97 86 + 97 96 107 98 87 + 98 97 108 99 88 + 99 98 109 -1 89 + 100 -1 110 101 90 + 101 100 111 102 91 + 102 101 112 103 92 + 103 102 113 104 93 + 104 103 114 105 94 + 105 104 115 106 95 + 106 105 116 107 96 + 107 106 117 108 97 + 108 107 118 109 98 + 109 108 119 -1 99 + 110 -1 -1 111 100 + 111 110 -1 112 101 + 112 111 -1 113 102 + 113 112 -1 114 103 + 114 113 -1 115 104 + 115 114 -1 116 105 + 116 115 -1 117 106 + 117 116 -1 118 107 + 118 117 -1 119 108 + 119 118 -1 -1 109 diff --git a/tests_cpp/gradients/info.dat b/tests_cpp/gradients/info.dat index cb765a86..0fde6d91 100644 --- a/tests_cpp/gradients/info.dat +++ b/tests_cpp/gradients/info.dat @@ -7,6 +7,6 @@ dx 0.10000000000000 dy 0.08333333333333 sampleMeshSize 120 stencilMeshSize 120 -stencilSize 5 +stencilSize 3 nx 10 ny 12 diff --git a/tests_cpp/gradients/main.cc b/tests_cpp/gradients/main.cc index 4cb197d4..08113329 100644 --- a/tests_cpp/gradients/main.cc +++ b/tests_cpp/gradients/main.cc @@ -1,57 +1,8 @@ #include #include "pressiodemoapps/mesh.hpp" +#include "pressiodemoapps/gradient.hpp" #include - // double x0 = cellX; - // double x1 = cellX + dx; - // double y1 = cellY; - // double f0 = f(row[0]); - // double f1 = f(row[3]); - // double alfa = -(f0-f1)/dx; - // double beta = f0 - alfa*x0; - // double fleft = alfa * bdFaceX + beta; - // gradVals.push_back( (-fleft + f1)/(2.*dx) ); - -template -bool cell_has_left_face_on_boundary(const MeshType & mesh, - int graphRow){ - // true if the first order neighboring cell index is -1 - const auto & G = mesh.graph(); - return (G(graphRow,1)==-1); -} - -template -bool cell_has_front_face_on_boundary(const MeshType & mesh, - int graphRow){ - const auto & G = mesh.graph(); - return (G(graphRow,2)==-1); -} - -template -bool cell_has_right_face_on_boundary(const MeshType & mesh, - int graphRow){ - const auto & G = mesh.graph(); - return (G(graphRow,3)==-1); -} - -template -bool cell_has_back_face_on_boundary(const MeshType & mesh, - int graphRow){ - const auto & G = mesh.graph(); - return (G(graphRow,4)==-1); -} - -template -bool cell_is_strictly_next_to_boundary(const MeshType & mesh, - int graphRow) -{ - const auto bL = cell_has_left_face_on_boundary(mesh, graphRow); - const auto bF = cell_has_front_face_on_boundary(mesh, graphRow); - const auto bR = cell_has_right_face_on_boundary(mesh, graphRow); - const auto bB = cell_has_back_face_on_boundary(mesh, graphRow); - return bL || bF || bR || bB; -} - template void sinef(const T & x, const T & y, StateType & f) { @@ -60,72 +11,6 @@ void sinef(const T & x, const T & y, StateType & f) } } -enum class Skew{Right2, Left2, Center3, Right3, Left3, Center5}; - -template -double fd_gradient_x(Skew direction, - int stencilSize, - const Connect & conn, - const double h, - const Eigen::VectorXd & f) -{ - if (stencilSize != 3 && stencilSize!=5){ - throw std::runtime_error("invalid stencil size"); - } - - // https://web.media.mit.edu/~crtaylor/calculator.html - switch(direction){ - case Skew::Right2: - return (-f(conn[0]) + f(conn[3]))/h; - case Skew::Left2: - return ( f(conn[0]) - f(conn[1]))/h; - case Skew::Center3: - return (-f(conn[1]) + f(conn[3]))/(2*h); - - case Skew::Right3: - return (-2.*f(conn[0]) + 3.*f(conn[3]) - 1.*f(conn[7]))/h; - case Skew::Left3: - return ( 2.*f(conn[0]) - 3.*f(conn[1]) + 1.*f(conn[5]))/h; - case Skew::Center5: - return (-f(conn[5]) - 8.*f(conn[1]) + 8.*f(conn[3]) - f(conn[7]))/(12.*h); - - default: - return 0; - } -} - -template -double fd_gradient_y(Skew direction, - int stencilSize, - const Connect & conn, - const double h, - const Eigen::VectorXd & f) -{ - if (stencilSize != 3 && stencilSize!=5){ - throw std::runtime_error("invalid stencil size"); - } - - // https://web.media.mit.edu/~crtaylor/calculator.html - switch(direction){ - case Skew::Right2: - return (-f(conn[0]) + f(conn[2]))/h; - case Skew::Left2: - return ( f(conn[0]) - f(conn[4]))/h; - case Skew::Center3: - return (-f(conn[4]) + f(conn[2]))/(2*h); - - case Skew::Right3: - return (-2.*f(conn[0]) + 3.*f(conn[2]) - 1.*f(conn[6]))/h; - case Skew::Left3: - return ( 2.*f(conn[0]) - 3.*f(conn[4]) + 1.*f(conn[8]))/h; - case Skew::Center5: - return (-f(conn[8]) - 8.*f(conn[4]) + 8.*f(conn[2]) - f(conn[6]))/(12.*h); - - default: - return 0; - } -} - constexpr int grad_x = 0; constexpr int grad_y = 1; @@ -143,6 +28,7 @@ auto one_sided_fd_normal_gradient_at_boundary_faces(const MeshType & mesh, const Eigen::VectorXd & f, int order) { + namespace pda = pressiodemoapps; const int ss = mesh.stencilSize(); const auto & x = mesh.viewX(); @@ -150,20 +36,15 @@ auto one_sided_fd_normal_gradient_at_boundary_faces(const MeshType & mesh, const auto dx = mesh.dx(); const auto dy = mesh.dy(); + const auto & rowsForLoop = mesh.graphRowsOfCellsStrictlyOnBd(); const auto & G = mesh.graph(); - std::vector rowsForLoop; - for (auto it : mesh.graphRowsOfCellsNearBd()){ - const bool b = cell_is_strictly_next_to_boundary(mesh, it); - if (b) rowsForLoop.push_back(it); - } - std::vector result; for (auto rowInd : rowsForLoop) { - const bool bdL = cell_has_left_face_on_boundary(mesh, rowInd); - const bool bdF = cell_has_front_face_on_boundary(mesh, rowInd); - const bool bdR = cell_has_right_face_on_boundary(mesh, rowInd); - const bool bdB = cell_has_back_face_on_boundary(mesh, rowInd); + const bool bdL = mesh.cellHasLeftFaceOnBoundary2d(rowInd); + const bool bdF = mesh.cellHasFrontFaceOnBoundary2d(rowInd); + const bool bdR = mesh.cellHasRightFaceOnBoundary2d(rowInd); + const bool bdB = mesh.cellHasBackFaceOnBoundary2d(rowInd); std::cout << "row = " << rowInd << " " << G(rowInd,0) << std::endl; const auto graphRow = G.row(rowInd); @@ -171,9 +52,8 @@ auto one_sided_fd_normal_gradient_at_boundary_faces(const MeshType & mesh, const auto cellX = x(cellGID); const auto cellY = y(cellGID); - const auto SkewRight = (order == 1) ? Skew::Right2 : Skew::Right3; - const auto SkewLeft = (order == 1) ? Skew::Left2 : Skew::Left3; - + const auto SkewRight = (order == 1) ? pda::GradFdMode::ForwardTwoPt : pda::GradFdMode::ForwardThreePt; + const auto SkewLeft = (order == 1) ? pda::GradFdMode::BackwardTwoPt : pda::GradFdMode::BackwardThreePt; Face face; face.parentCellGID = cellGID; face.parentCellGraphRow = rowInd; @@ -181,73 +61,45 @@ auto one_sided_fd_normal_gradient_at_boundary_faces(const MeshType & mesh, if (bdL){ face.x = cellX-dx*0.5; face.y = cellY; - face.gradValue = fd_gradient_x(SkewRight, ss, graphRow, dx, f); + face.gradValue = pda::impl::fd_face_normal_gradient_for_cell_centered_state_2d(f, graphRow, 'x', + SkewRight, dx); face.gradDirection = grad_x; result.push_back(face); - - if (bdB){ - face.x = cellX; - face.y = cellY-dy*0.5; - face.gradValue = fd_gradient_y(SkewRight, ss, graphRow, dy, f); - face.gradDirection = grad_y; - result.push_back(face); - } - - if (bdF){ - face.x = cellX; - face.y = cellY+dy*0.5; - face.gradValue = fd_gradient_y(SkewLeft, ss, graphRow, dy, f); - face.gradDirection = grad_y; - result.push_back(face); - } } - if (bdF && !bdL && !bdR){ - face.x = cellX; - face.y = cellY+dy*0.5; - face.gradValue = fd_gradient_y(SkewLeft, ss, graphRow, dy, f); - face.gradDirection = grad_y; - result.push_back(face); - } + // if (bdF){ + // face.x = cellX; + // face.y = cellY+dy*0.5; + // face.gradValue = pda::impl::fd_face_normal_gradient_for_cell_centered_state_2d(f, graphRow, 'y', + // SkewLeft, dy); + // face.gradDirection = grad_y; + // result.push_back(face); + // } if (bdR){ face.x = cellX+dx*0.5; face.y = cellY; - face.gradValue = fd_gradient_x(SkewLeft, ss, graphRow, dx, f); + face.gradValue = pda::impl::fd_face_normal_gradient_for_cell_centered_state_2d(f, graphRow, 'x', + SkewLeft, dx); face.gradDirection = grad_x; result.push_back(face); - - if (bdB){ - face.x = cellX; - face.y = cellY-dy*0.5; - face.gradValue = fd_gradient_y(SkewRight, ss, graphRow, dy, f); - face.gradDirection = grad_y; - result.push_back(face); - } - - if (bdF){ - face.x = cellX; - face.y = cellY+dy*0.5; - face.gradValue = fd_gradient_y(SkewLeft, ss, graphRow, dy, f); - face.gradDirection = grad_y; - result.push_back(face); - } } - if (bdB && !bdL && !bdR){ - face.x = cellX; - face.y = cellY-dy*0.5; - face.gradValue = fd_gradient_y(SkewRight, ss, graphRow, dy, f); - face.gradDirection = grad_y; - result.push_back(face); - } + // if (bdB){ + // face.x = cellX; + // face.y = cellY-dy*0.5; + // face.gradValue = pda::impl::fd_face_normal_gradient_for_cell_centered_state_2d(f, graphRow, 'y', + // SkewRight, dy); + // face.gradDirection = grad_y; + // result.push_back(face); + // } } return result; } void write_result(const std::vector & M, - const std::string & filename) + const std::string & filename) { std::ofstream file; file.open(filename); for (auto & it : M){ @@ -258,8 +110,8 @@ void write_result(const std::vector & M, << it.gradDirection << std::endl; } + file.close(); }; -//file << std::setprecision(14); int main() { @@ -275,22 +127,45 @@ int main() const auto fd_result1 = one_sided_fd_normal_gradient_at_boundary_faces(mesh, f, 1); write_result(fd_result1, "fd_results_1.txt"); - const auto fd_result2 = one_sided_fd_normal_gradient_at_boundary_faces(mesh, f, 2); - write_result(fd_result2, "fd_results_2.txt"); + // const auto fd_result2 = one_sided_fd_normal_gradient_at_boundary_faces(mesh, f, 2); + // write_result(fd_result2, "fd_results_2.txt"); - // std::vector< std::array > fd_grads; - // const auto & G = mesh.graph(); - // for (int i=0; i