diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.inl b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.inl index fd91291b..66e34b16 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.inl @@ -42,6 +42,7 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int center_nz_index = getSolverMatrixIndex(i_r, i_theta); const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta); @@ -63,7 +64,7 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ updateMatrixElement(solver_matrix, ptr, offset, row, col, val); offset = CenterStencil[StencilPosition::Left]; @@ -269,6 +270,7 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int center_nz_index = getSolverMatrixIndex(i_r, i_theta); const int left_nz_index = getSolverMatrixIndex(i_r, i_theta_AcrossOrigin); @@ -289,7 +291,7 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ updateMatrixElement(solver_matrix, ptr, offset, row, col, val); offset = CenterStencil[StencilPosition::Left]; @@ -424,6 +426,7 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -448,7 +451,7 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ updateMatrixElement(solver_matrix, ptr, offset, row, col, val); /* REMOVED: Moved to the right hand side to make the matrix symmetric */ @@ -604,6 +607,7 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -628,7 +632,7 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ updateMatrixElement(solver_matrix, ptr, offset, row, col, val); offset = CenterStencil[StencilPosition::Left]; diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.inl b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.inl index 9f952eb8..1e5bf9da 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.inl @@ -41,6 +41,7 @@ void DirectSolver_COO_MUMPS_Take::nodeBuildSolverMatrixTake( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int center_nz_index = getSolverMatrixIndex(i_r, i_theta); @@ -59,13 +60,12 @@ void DirectSolver_COO_MUMPS_Take::nodeBuildSolverMatrixTake( double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ - double center_value = - (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + double center_value = (coeff5 * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ @@ -168,6 +168,7 @@ void DirectSolver_COO_MUMPS_Take::nodeBuildSolverMatrixTake( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int center_nz_index = getSolverMatrixIndex(i_r, i_theta); @@ -184,13 +185,12 @@ void DirectSolver_COO_MUMPS_Take::nodeBuildSolverMatrixTake( double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ - double center_value = - (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + double center_value = (coeff5 * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ @@ -257,6 +257,7 @@ void DirectSolver_COO_MUMPS_Take::nodeBuildSolverMatrixTake( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int center_nz_index = getSolverMatrixIndex(i_r, i_theta); @@ -275,13 +276,12 @@ void DirectSolver_COO_MUMPS_Take::nodeBuildSolverMatrixTake( double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ - double center_value = - (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + double center_value = (coeff5 * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ @@ -364,6 +364,7 @@ void DirectSolver_COO_MUMPS_Take::nodeBuildSolverMatrixTake( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int center_nz_index = getSolverMatrixIndex(i_r, i_theta); @@ -382,13 +383,12 @@ void DirectSolver_COO_MUMPS_Take::nodeBuildSolverMatrixTake( double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ - double center_value = - (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + double center_value = (coeff5 * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ diff --git a/include/DirectSolver/DirectSolver-CSR-LU-Give/buildSolverMatrix.inl b/include/DirectSolver/DirectSolver-CSR-LU-Give/buildSolverMatrix.inl index 55cdb111..944c356b 100644 --- a/include/DirectSolver/DirectSolver-CSR-LU-Give/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolver-CSR-LU-Give/buildSolverMatrix.inl @@ -38,6 +38,7 @@ void DirectSolver_CSR_LU_Give::nodeBuildSolverMatrixGive(int i_r const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int center_index = grid.index(i_r, i_theta); const int left_index = grid.index(i_r - 1, i_theta); @@ -52,7 +53,7 @@ void DirectSolver_CSR_LU_Give::nodeBuildSolverMatrixGive(int i_r offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ updateMatrixElement(solver_matrix, offset, row, col, val); offset = CenterStencil[StencilPosition::Left]; @@ -250,14 +251,16 @@ void DirectSolver_CSR_LU_Give::nodeBuildSolverMatrixGive(int i_r assert(grid.ntheta() % 2 == 0); const int i_theta_AcrossOrigin = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); - double h1 = 2.0 * grid.radius(0); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta_M1); - double k2 = grid.angularSpacing(i_theta); + double h1 = 2.0 * grid.radius(0); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta_M1); + double k2 = grid.angularSpacing(i_theta); + double coeff1 = 0.5 * (k1 + k2) / h1; double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int center_index = grid.index(i_r, i_theta); const int left_index = grid.index(i_r, i_theta_AcrossOrigin); @@ -272,7 +275,7 @@ void DirectSolver_CSR_LU_Give::nodeBuildSolverMatrixGive(int i_r offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ updateMatrixElement(solver_matrix, offset, row, col, val); offset = CenterStencil[StencilPosition::Left]; @@ -404,6 +407,7 @@ void DirectSolver_CSR_LU_Give::nodeBuildSolverMatrixGive(int i_r const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -421,7 +425,7 @@ void DirectSolver_CSR_LU_Give::nodeBuildSolverMatrixGive(int i_r offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ updateMatrixElement(solver_matrix, offset, row, col, val); offset = CenterStencil[StencilPosition::Left]; @@ -564,6 +568,7 @@ void DirectSolver_CSR_LU_Give::nodeBuildSolverMatrixGive(int i_r const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -581,7 +586,7 @@ void DirectSolver_CSR_LU_Give::nodeBuildSolverMatrixGive(int i_r offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ updateMatrixElement(solver_matrix, offset, row, col, val); offset = CenterStencil[StencilPosition::Left]; diff --git a/include/DirectSolver/DirectSolver-CSR-LU-Take/buildSolverMatrix.inl b/include/DirectSolver/DirectSolver-CSR-LU-Take/buildSolverMatrix.inl index b9a31e36..db66a24d 100644 --- a/include/DirectSolver/DirectSolver-CSR-LU-Take/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolver-CSR-LU-Take/buildSolverMatrix.inl @@ -37,6 +37,7 @@ void DirectSolver_CSR_LU_Take::nodeBuildSolverMatrixTake( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int center_index = grid.index(i_r, i_theta); int left_index = grid.index(i_r - 1, i_theta); @@ -53,13 +54,12 @@ void DirectSolver_CSR_LU_Take::nodeBuildSolverMatrixTake( double bottom_value = -coeff3 * (att[center_index] + att[bottom_index]); /* Bottom */ double top_value = -coeff4 * (att[center_index] + att[top_index]); /* Top */ - double center_value = - (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * std::fabs(detDF[center_index]) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + double center_value = (coeff5 * coeff_beta[center_index] * std::fabs(detDF[center_index]) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); double bottom_left_value = -0.25 * (art[left_index] + art[bottom_index]); /* Bottom Left */ double bottom_right_value = +0.25 * (art[right_index] + art[bottom_index]); /* Bottom Right */ @@ -158,6 +158,7 @@ void DirectSolver_CSR_LU_Take::nodeBuildSolverMatrixTake( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int center_index = grid.index(i_r, i_theta); int left_index = grid.index(i_r, i_theta_AcrossOrigin); @@ -172,13 +173,12 @@ void DirectSolver_CSR_LU_Take::nodeBuildSolverMatrixTake( double bottom_value = -coeff3 * (att[center_index] + att[bottom_index]); /* Bottom */ double top_value = -coeff4 * (att[center_index] + att[top_index]); /* Top */ - double center_value = - (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF[center_index]) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + double center_value = (coeff5 * coeff_beta[center_index] * fabs(detDF[center_index]) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); double bottom_right_value = +0.25 * (art[right_index] + art[bottom_index]); /* Bottom Right */ double top_right_value = -0.25 * (art[right_index] + art[top_index]); /* Top Right */ diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.inl index f5310376..a0a1f146 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.inl @@ -70,6 +70,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -108,7 +109,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -207,7 +208,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -260,6 +261,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -302,7 +304,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -389,7 +391,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -521,7 +523,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ updateCOOCSRMatrixElement(inner_boundary_circle_matrix, ptr, offset, row, col, val); offset = CenterStencil[StencilPosition::Left]; @@ -620,6 +622,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -649,7 +652,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -712,7 +715,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -765,7 +768,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -838,6 +841,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -869,7 +873,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -922,7 +926,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -965,7 +969,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -1040,6 +1044,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -1070,7 +1075,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -1124,7 +1129,7 @@ void ExtrapolatedSmootherGive::nodeBuildAscGive( /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl index 373a5c74..2bfc0d70 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl @@ -76,6 +76,7 @@ void ExtrapolatedSmootherTake::nodeBuildInteriorBoundarySolverMa double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -109,10 +110,10 @@ void ExtrapolatedSmootherTake::nodeBuildInteriorBoundarySolverMa const int top = grid.index(i_r, i_theta_P1); const int right = grid.index(i_r + 1, i_theta); - const double center_value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) + + const double center_value = coeff5 * coeff_beta[center] * fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); - const double left_value = -coeff1 * (arr[center] + arr[left]); + const double left_value = -coeff1 * (arr[center] + arr[left]); /* Fill matrix row of (i,j) */ row = center_index; diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl index d487aa00..79cbabcd 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl @@ -52,6 +52,7 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -89,9 +90,9 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) + - coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + - coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta[center] * fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); /* Bottom */ @@ -128,9 +129,9 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) + - coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + - coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta[center] * fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); } else { @@ -190,6 +191,7 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -232,9 +234,9 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) + - coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + - coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta[center] * fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); /* Left */ @@ -276,9 +278,15 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; +<<<<<<< HEAD value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) + - coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + - coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); + coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + + coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); +======= + value = coeff5 * coeff_beta[center] * fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); +>>>>>>> 3d8ce52 (litz_coeff5 changes) updateMatrixElement(solver, batch, row, column, value); } else { @@ -306,6 +314,7 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -340,9 +349,15 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; +<<<<<<< HEAD value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) + - coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + - coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); + coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + + coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); +======= + value = coeff5 * coeff_beta[center] * fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); +>>>>>>> 3d8ce52 (litz_coeff5 changes) updateMatrixElement(solver, batch, row, column, value); /* Right */ @@ -363,9 +378,15 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; +<<<<<<< HEAD value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) + - coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + - coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); + coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + + coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); +======= + value = coeff5 * coeff_beta[center] * fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); +>>>>>>> 3d8ce52 (litz_coeff5 changes) updateMatrixElement(solver, batch, row, column, value); } else { @@ -399,6 +420,7 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -429,9 +451,15 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; +<<<<<<< HEAD value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) + - coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + - coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); + coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + + coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); +======= + value = coeff5 * coeff_beta[center] * fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); +>>>>>>> 3d8ce52 (litz_coeff5 changes) updateMatrixElement(solver, batch, row, column, value); /* Left */ @@ -457,9 +485,15 @@ void ExtrapolatedSmootherTake::nodeBuildTridiagonalSolverMatrice /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; +<<<<<<< HEAD value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) + - coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + - coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); + coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + + coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); +======= + value = coeff5 * coeff_beta[center] * fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); +>>>>>>> 3d8ce52 (litz_coeff5 changes) updateMatrixElement(solver, batch, row, column, value); } } diff --git a/include/Residual/ResidualGive/applyAGive.inl b/include/Residual/ResidualGive/applyAGive.inl index 2802e0ef..ae4ddafa 100644 --- a/include/Residual/ResidualGive/applyAGive.inl +++ b/include/Residual/ResidualGive/applyAGive.inl @@ -28,6 +28,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -38,7 +39,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet const int top = grid.index(i_r, i_theta_P1); /* Fill result(i,j) */ - result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] += (coeff5 * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ @@ -113,6 +114,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -124,7 +126,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet const int top = grid.index(i_r, i_theta_P1); /* Fill result(i,j) */ - result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] += (coeff5 * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ @@ -167,6 +169,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -177,7 +180,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet const int top = grid.index(i_r, i_theta_P1); /* Fill result(i,j) */ - result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] += (coeff5 * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ @@ -220,6 +223,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -230,7 +234,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet const int top = grid.index(i_r, i_theta_P1); /* Fill result(i,j) */ - result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] += (coeff5 * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index 38cbe669..8475791d 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -21,6 +21,7 @@ static inline void node_apply_a_take(int i_r, int i_theta, const PolarGrid& grid double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -35,18 +36,18 @@ static inline void node_apply_a_take(int i_r, int i_theta, const PolarGrid& grid const int right = grid.index(i_r + 1, i_theta); const int top_right = grid.index(i_r + 1, i_theta_P1); - result[center] = - 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) * x[center] /* beta_{i,j} */ + result[center] = coeff5 * coeff_beta[center] * std::fabs(detDF[center]) * x[center] /* beta_{i,j} */ - - coeff1 * (arr[center] + arr[left]) * (x[left] - x[center]) /* Left - Center: (Left) */ - - coeff2 * (arr[center] + arr[right]) * (x[right] - x[center]) /* Right - Center: (Right) */ - - coeff3 * (att[center] + att[bottom]) * (x[bottom] - x[center]) /* Bottom - Center: (Bottom) */ - - coeff4 * (att[center] + att[top]) * (x[top] - x[center]) /* Top - Center: (Top) */ + - coeff1 * (arr[center] + arr[left]) * (x[left] - x[center]) /* Left - Center: (Left) */ + - coeff2 * (arr[center] + arr[right]) * (x[right] - x[center]) /* Right - Center: (Right) */ + - + coeff3 * (att[center] + att[bottom]) * (x[bottom] - x[center]) /* Bottom - Center: (Bottom) */ + - coeff4 * (att[center] + att[top]) * (x[top] - x[center]) /* Top - Center: (Top) */ - - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */ - - 0.25 * (art[right] + art[top]) * x[top_right]; /* Top Right */ + - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ + + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */ + - 0.25 * (art[right] + art[top]) * x[top_right]; /* Top Right */ } /* -------------------------- */ /* Node on the inner boundary */ @@ -75,6 +76,7 @@ static inline void node_apply_a_take(int i_r, int i_theta, const PolarGrid& grid double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -89,8 +91,7 @@ static inline void node_apply_a_take(int i_r, int i_theta, const PolarGrid& grid const int top_right = grid.index(i_r + 1, i_theta_P1); result[center] = - 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) * - x[center] /* beta_{i,j} */ + coeff5 * coeff_beta[center] * std::fabs(detDF[center]) * x[center] /* beta_{i,j} */ - coeff1 * (arr[center] + arr[left]) * (x[left] - x[center]) /* Left - Center: (Left) */ - coeff2 * (arr[center] + arr[right]) * (x[right] - x[center]) /* Right - Center: (Right) */ diff --git a/include/Smoother/SmootherGive/buildMatrix.inl b/include/Smoother/SmootherGive/buildMatrix.inl index 7a77fc1c..3c7e4721 100644 --- a/include/Smoother/SmootherGive/buildMatrix.inl +++ b/include/Smoother/SmootherGive/buildMatrix.inl @@ -71,6 +71,7 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -105,7 +106,7 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -183,6 +184,7 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -210,7 +212,7 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -321,6 +323,7 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); /* left_matrix (across-the origin), center_matrix, right_matrix */ /* -| x | o | x | */ @@ -354,7 +357,7 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const offset = CenterStencil[StencilPosition::Center]; col = center_index; - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ + val = coeff5 * coeff_beta * fabs(detDF); /* beta_{i,j} */ updateCOOCSRMatrixElement(inner_boundary_circle_matrix, ptr, offset, row, col, val); offset = CenterStencil[StencilPosition::Left]; @@ -450,10 +453,11 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const double k1 = grid.angularSpacing(i_theta - 1); double k2 = grid.angularSpacing(i_theta); - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - double coeff4 = 0.5 * (h1 + h2) / k2; + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -481,7 +485,7 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; @@ -536,6 +540,7 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -564,7 +569,7 @@ void SmootherGive::nodeBuildAscGive(int i_r, int i_theta, const /* ---------------------------- */ /* Fill matrix row of (i,j) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ + value = coeff5 * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ updateMatrixElement(center_solver, center_batch, row, column, value); row = center_index; diff --git a/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl b/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl index 978acdc2..7bed23e6 100644 --- a/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl +++ b/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl @@ -76,6 +76,7 @@ void SmootherTake::nodeBuildInteriorBoundarySolverMatrix( const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -99,7 +100,7 @@ void SmootherTake::nodeBuildInteriorBoundarySolverMatrix( const double bottom_value = -coeff3 * (att[center] + att[bottom]); const double top_value = -coeff4 * (att[center] + att[top]); - const double center_value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) - + const double center_value = coeff5 * coeff_beta[center] * std::fabs(detDF[center]) - (left_value + right_value + bottom_value + top_value); /* Fill matrix row of (i,j) */ diff --git a/include/Smoother/SmootherTake/buildTridiagonalAsc.inl b/include/Smoother/SmootherTake/buildTridiagonalAsc.inl index 7134b2ae..e28229a4 100644 --- a/include/Smoother/SmootherTake/buildTridiagonalAsc.inl +++ b/include/Smoother/SmootherTake/buildTridiagonalAsc.inl @@ -49,6 +49,7 @@ void SmootherTake::nodeBuildTridiagonalSolverMatrices( const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -74,8 +75,8 @@ void SmootherTake::nodeBuildTridiagonalSolverMatrices( /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) - - (left_value + right_value + bottom_value + top_value); + value = coeff5 * coeff_beta[center] * std::fabs(detDF[center]) - + (left_value + right_value + bottom_value + top_value); updateMatrixElement(solver, batch, row, column, value); /* Bottom */ @@ -103,6 +104,7 @@ void SmootherTake::nodeBuildTridiagonalSolverMatrices( const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -128,8 +130,8 @@ void SmootherTake::nodeBuildTridiagonalSolverMatrices( /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) - - (left_value + right_value + bottom_value + top_value); + value = coeff5 * coeff_beta[center] * std::fabs(detDF[center]) - + (left_value + right_value + bottom_value + top_value); updateMatrixElement(solver, batch, row, column, value); /* Left */ @@ -190,6 +192,7 @@ void SmootherTake::nodeBuildTridiagonalSolverMatrices( const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -214,8 +217,8 @@ void SmootherTake::nodeBuildTridiagonalSolverMatrices( /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) - - (left_value + right_value + bottom_value + top_value); + value = coeff5 * coeff_beta[center] * std::fabs(detDF[center]) - + (left_value + right_value + bottom_value + top_value); updateMatrixElement(solver, batch, row, column, value); /* Right */ @@ -237,6 +240,7 @@ void SmootherTake::nodeBuildTridiagonalSolverMatrices( const double coeff2 = 0.5 * (k1 + k2) / h2; const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -262,8 +266,8 @@ void SmootherTake::nodeBuildTridiagonalSolverMatrices( /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) - - (left_value + right_value + bottom_value + top_value); + value = coeff5 * coeff_beta[center] * std::fabs(detDF[center]) - + (left_value + right_value + bottom_value + top_value); updateMatrixElement(solver, batch, row, column, value); /* Left */