diff --git a/tests_cpp/CMakeLists.txt b/tests_cpp/CMakeLists.txt index 50e43c21..6c3af1e4 100644 --- a/tests_cpp/CMakeLists.txt +++ b/tests_cpp/CMakeLists.txt @@ -27,6 +27,8 @@ add_subdirectory(eigen_rusanov_flux_jacobians_swe) add_subdirectory(create_vec_of_problems) +add_subdirectory(gradients) + # --------------------------------------------------------- # 1d problems # --------------------------------------------------------- diff --git a/tests_cpp/gradients/CMakeLists.txt b/tests_cpp/gradients/CMakeLists.txt new file mode 100644 index 00000000..5625b608 --- /dev/null +++ b/tests_cpp/gradients/CMakeLists.txt @@ -0,0 +1,14 @@ + +set(testname gradients) +set(exename ${testname}_exe) + +configure_file(info.dat info.dat COPYONLY) +configure_file(connectivity.dat connectivity.dat COPYONLY) +configure_file(coordinates.dat coordinates.dat COPYONLY) + +add_executable(${exename} ${CMAKE_CURRENT_SOURCE_DIR}/main.cc) +add_test(NAME ${testname} COMMAND ${exename}) +set_tests_properties(${testname} + PROPERTIES PASS_REGULAR_EXPRESSION "PASS" + FAIL_REGULAR_EXPRESSION "FAILED" + ) diff --git a/tests_cpp/gradients/connectivity.dat b/tests_cpp/gradients/connectivity.dat new file mode 100644 index 00000000..c73879ab --- /dev/null +++ b/tests_cpp/gradients/connectivity.dat @@ -0,0 +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 diff --git a/tests_cpp/gradients/coordinates.dat b/tests_cpp/gradients/coordinates.dat new file mode 100644 index 00000000..e6278725 --- /dev/null +++ b/tests_cpp/gradients/coordinates.dat @@ -0,0 +1,120 @@ + 0 0.05000000000000 0.04166666666667 + 1 0.15000000000000 0.04166666666667 + 2 0.25000000000000 0.04166666666667 + 3 0.35000000000000 0.04166666666667 + 4 0.45000000000000 0.04166666666667 + 5 0.55000000000000 0.04166666666667 + 6 0.65000000000000 0.04166666666667 + 7 0.75000000000000 0.04166666666667 + 8 0.85000000000000 0.04166666666667 + 9 0.95000000000000 0.04166666666667 + 10 0.05000000000000 0.12500000000000 + 11 0.15000000000000 0.12500000000000 + 12 0.25000000000000 0.12500000000000 + 13 0.35000000000000 0.12500000000000 + 14 0.45000000000000 0.12500000000000 + 15 0.55000000000000 0.12500000000000 + 16 0.65000000000000 0.12500000000000 + 17 0.75000000000000 0.12500000000000 + 18 0.85000000000000 0.12500000000000 + 19 0.95000000000000 0.12500000000000 + 20 0.05000000000000 0.20833333333333 + 21 0.15000000000000 0.20833333333333 + 22 0.25000000000000 0.20833333333333 + 23 0.35000000000000 0.20833333333333 + 24 0.45000000000000 0.20833333333333 + 25 0.55000000000000 0.20833333333333 + 26 0.65000000000000 0.20833333333333 + 27 0.75000000000000 0.20833333333333 + 28 0.85000000000000 0.20833333333333 + 29 0.95000000000000 0.20833333333333 + 30 0.05000000000000 0.29166666666667 + 31 0.15000000000000 0.29166666666667 + 32 0.25000000000000 0.29166666666667 + 33 0.35000000000000 0.29166666666667 + 34 0.45000000000000 0.29166666666667 + 35 0.55000000000000 0.29166666666667 + 36 0.65000000000000 0.29166666666667 + 37 0.75000000000000 0.29166666666667 + 38 0.85000000000000 0.29166666666667 + 39 0.95000000000000 0.29166666666667 + 40 0.05000000000000 0.37500000000000 + 41 0.15000000000000 0.37500000000000 + 42 0.25000000000000 0.37500000000000 + 43 0.35000000000000 0.37500000000000 + 44 0.45000000000000 0.37500000000000 + 45 0.55000000000000 0.37500000000000 + 46 0.65000000000000 0.37500000000000 + 47 0.75000000000000 0.37500000000000 + 48 0.85000000000000 0.37500000000000 + 49 0.95000000000000 0.37500000000000 + 50 0.05000000000000 0.45833333333333 + 51 0.15000000000000 0.45833333333333 + 52 0.25000000000000 0.45833333333333 + 53 0.35000000000000 0.45833333333333 + 54 0.45000000000000 0.45833333333333 + 55 0.55000000000000 0.45833333333333 + 56 0.65000000000000 0.45833333333333 + 57 0.75000000000000 0.45833333333333 + 58 0.85000000000000 0.45833333333333 + 59 0.95000000000000 0.45833333333333 + 60 0.05000000000000 0.54166666666667 + 61 0.15000000000000 0.54166666666667 + 62 0.25000000000000 0.54166666666667 + 63 0.35000000000000 0.54166666666667 + 64 0.45000000000000 0.54166666666667 + 65 0.55000000000000 0.54166666666667 + 66 0.65000000000000 0.54166666666667 + 67 0.75000000000000 0.54166666666667 + 68 0.85000000000000 0.54166666666667 + 69 0.95000000000000 0.54166666666667 + 70 0.05000000000000 0.62500000000000 + 71 0.15000000000000 0.62500000000000 + 72 0.25000000000000 0.62500000000000 + 73 0.35000000000000 0.62500000000000 + 74 0.45000000000000 0.62500000000000 + 75 0.55000000000000 0.62500000000000 + 76 0.65000000000000 0.62500000000000 + 77 0.75000000000000 0.62500000000000 + 78 0.85000000000000 0.62500000000000 + 79 0.95000000000000 0.62500000000000 + 80 0.05000000000000 0.70833333333333 + 81 0.15000000000000 0.70833333333333 + 82 0.25000000000000 0.70833333333333 + 83 0.35000000000000 0.70833333333333 + 84 0.45000000000000 0.70833333333333 + 85 0.55000000000000 0.70833333333333 + 86 0.65000000000000 0.70833333333333 + 87 0.75000000000000 0.70833333333333 + 88 0.85000000000000 0.70833333333333 + 89 0.95000000000000 0.70833333333333 + 90 0.05000000000000 0.79166666666667 + 91 0.15000000000000 0.79166666666667 + 92 0.25000000000000 0.79166666666667 + 93 0.35000000000000 0.79166666666667 + 94 0.45000000000000 0.79166666666667 + 95 0.55000000000000 0.79166666666667 + 96 0.65000000000000 0.79166666666667 + 97 0.75000000000000 0.79166666666667 + 98 0.85000000000000 0.79166666666667 + 99 0.95000000000000 0.79166666666667 + 100 0.05000000000000 0.87500000000000 + 101 0.15000000000000 0.87500000000000 + 102 0.25000000000000 0.87500000000000 + 103 0.35000000000000 0.87500000000000 + 104 0.45000000000000 0.87500000000000 + 105 0.55000000000000 0.87500000000000 + 106 0.65000000000000 0.87500000000000 + 107 0.75000000000000 0.87500000000000 + 108 0.85000000000000 0.87500000000000 + 109 0.95000000000000 0.87500000000000 + 110 0.05000000000000 0.95833333333333 + 111 0.15000000000000 0.95833333333333 + 112 0.25000000000000 0.95833333333333 + 113 0.35000000000000 0.95833333333333 + 114 0.45000000000000 0.95833333333333 + 115 0.55000000000000 0.95833333333333 + 116 0.65000000000000 0.95833333333333 + 117 0.75000000000000 0.95833333333333 + 118 0.85000000000000 0.95833333333333 + 119 0.95000000000000 0.95833333333333 diff --git a/tests_cpp/gradients/info.dat b/tests_cpp/gradients/info.dat new file mode 100644 index 00000000..cb765a86 --- /dev/null +++ b/tests_cpp/gradients/info.dat @@ -0,0 +1,12 @@ +dim 2 +xMin 0.00000000000000 +xMax 1.00000000000000 +yMin 0.00000000000000 +yMax 1.00000000000000 +dx 0.10000000000000 +dy 0.08333333333333 +sampleMeshSize 120 +stencilMeshSize 120 +stencilSize 5 +nx 10 +ny 12 diff --git a/tests_cpp/gradients/main.cc b/tests_cpp/gradients/main.cc new file mode 100644 index 00000000..4cb197d4 --- /dev/null +++ b/tests_cpp/gradients/main.cc @@ -0,0 +1,299 @@ +#include +#include "pressiodemoapps/mesh.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) +{ + for (int i=0; i +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; + +struct Face{ + int parentCellGID; + int parentCellGraphRow; + double x; + double y; + double gradValue; + int gradDirection; +}; + +template +auto one_sided_fd_normal_gradient_at_boundary_faces(const MeshType & mesh, + const Eigen::VectorXd & f, + int order) +{ + + const int ss = mesh.stencilSize(); + const auto & x = mesh.viewX(); + const auto & y = mesh.viewY(); + const auto dx = mesh.dx(); + const auto dy = mesh.dy(); + + 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); + 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 SkewRight = (order == 1) ? Skew::Right2 : Skew::Right3; + const auto SkewLeft = (order == 1) ? Skew::Left2 : Skew::Left3; + + Face face; + face.parentCellGID = cellGID; + face.parentCellGraphRow = rowInd; + + if (bdL){ + face.x = cellX-dx*0.5; + face.y = cellY; + face.gradValue = fd_gradient_x(SkewRight, ss, graphRow, dx, f); + 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 (bdR){ + face.x = cellX+dx*0.5; + face.y = cellY; + face.gradValue = fd_gradient_x(SkewLeft, ss, graphRow, dx, f); + 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); + } + } + + return result; +} + +void write_result(const std::vector & M, + const std::string & filename) +{ + std::ofstream file; file.open(filename); + for (auto & it : M){ + file << it.parentCellGID << " " + << it.parentCellGraphRow << " " + << it.x << " " << it.y << " " + << it.gradValue << " " + << it.gradDirection + << std::endl; + } +}; +//file << std::setprecision(14); + +int main() +{ + namespace pda = pressiodemoapps; + + auto mesh = pda::load_cellcentered_uniform_mesh_eigen("."); + + const auto & x = mesh.viewX(); + const auto & y = mesh.viewY(); + Eigen::VectorXd f(mesh.stencilMeshSize()); + sinef(x, y, f); + + 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"); + + // std::vector< std::array > fd_grads; + // const auto & G = mesh.graph(); + // for (int i=0; i