From 5f63c1695cb98640c130f5275c9014bc088b4acf Mon Sep 17 00:00:00 2001 From: rezgarshakeri Date: Thu, 1 Dec 2022 11:41:59 -0700 Subject: [PATCH] some clean-up in Hdiv-mass folder --- examples/Hdiv-mass/Makefile | 2 +- examples/Hdiv-mass/basis/L2-P0.h | 58 ++++ .../Hdiv-mass/{conv_plot.py => conv_rate.py} | 28 +- examples/Hdiv-mass/conv_test.sh | 14 +- examples/Hdiv-mass/conv_test_result.csv | 10 - examples/Hdiv-mass/include/cl-options.h | 7 +- examples/Hdiv-mass/include/matops.h | 13 - examples/Hdiv-mass/include/petsc-macros.h | 17 -- examples/Hdiv-mass/include/post-processing.h | 12 + examples/Hdiv-mass/include/problems.h | 19 -- examples/Hdiv-mass/include/register-problem.h | 18 ++ examples/Hdiv-mass/include/setup-dm.h | 6 +- examples/Hdiv-mass/include/setup-fe.h | 20 ++ examples/Hdiv-mass/include/setup-libceed.h | 15 +- examples/Hdiv-mass/include/setup-matops.h | 13 + examples/Hdiv-mass/include/setup-solvers.h | 16 + examples/Hdiv-mass/include/structs.h | 78 ++--- examples/Hdiv-mass/main.c | 287 +++++------------- examples/Hdiv-mass/main.h | 7 +- examples/Hdiv-mass/problems/mass2d.c | 17 +- examples/Hdiv-mass/problems/mass3d.c | 19 +- .../Hdiv-mass/problems/register-problem.c | 33 ++ .../Hdiv-mass/qfunctions/poisson-error2d.h | 5 +- .../Hdiv-mass/qfunctions/poisson-error3d.h | 7 +- examples/Hdiv-mass/src/cl-options.c | 47 +-- examples/Hdiv-mass/src/matops.c | 106 ------- examples/Hdiv-mass/src/post-processing.c | 77 +++++ examples/Hdiv-mass/src/setup-dm.c | 89 ++---- examples/Hdiv-mass/src/setup-fe.c | 149 +++++++++ examples/Hdiv-mass/src/setup-libceed.c | 184 +++-------- examples/Hdiv-mass/src/setup-matops.c | 51 ++++ examples/Hdiv-mass/src/setup-solvers.c | 123 ++++++++ 32 files changed, 806 insertions(+), 741 deletions(-) create mode 100644 examples/Hdiv-mass/basis/L2-P0.h rename examples/Hdiv-mass/{conv_plot.py => conv_rate.py} (75%) delete mode 100644 examples/Hdiv-mass/conv_test_result.csv delete mode 100644 examples/Hdiv-mass/include/matops.h delete mode 100644 examples/Hdiv-mass/include/petsc-macros.h create mode 100644 examples/Hdiv-mass/include/post-processing.h delete mode 100644 examples/Hdiv-mass/include/problems.h create mode 100644 examples/Hdiv-mass/include/register-problem.h create mode 100644 examples/Hdiv-mass/include/setup-fe.h create mode 100644 examples/Hdiv-mass/include/setup-matops.h create mode 100644 examples/Hdiv-mass/include/setup-solvers.h create mode 100644 examples/Hdiv-mass/problems/register-problem.c delete mode 100644 examples/Hdiv-mass/src/matops.c create mode 100644 examples/Hdiv-mass/src/post-processing.c create mode 100644 examples/Hdiv-mass/src/setup-fe.c create mode 100644 examples/Hdiv-mass/src/setup-matops.c create mode 100644 examples/Hdiv-mass/src/setup-solvers.c diff --git a/examples/Hdiv-mass/Makefile b/examples/Hdiv-mass/Makefile index 6cfc950525..eecd42bcd2 100644 --- a/examples/Hdiv-mass/Makefile +++ b/examples/Hdiv-mass/Makefile @@ -69,7 +69,7 @@ print: $(PETSc.pc) $(ceed.pc) @true clean: - $(RM) -r $(OBJDIR) main *.vtu + $(RM) -r $(OBJDIR) main *.vtu *.csv $(PETSc.pc): $(if $(wildcard $@),,$(error \ diff --git a/examples/Hdiv-mass/basis/L2-P0.h b/examples/Hdiv-mass/basis/L2-P0.h new file mode 100644 index 0000000000..6149d4a34d --- /dev/null +++ b/examples/Hdiv-mass/basis/L2-P0.h @@ -0,0 +1,58 @@ +// Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. +// Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. +// All Rights reserved. See files LICENSE and NOTICE for details. +// +// This file is part of CEED, a collection of benchmarks, miniapps, software +// libraries and APIs for efficient high-order finite element and spectral +// element discretizations for exascale applications. For more information and +// source code availability see http://github.com/ceed. +// +// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, +// a collaborative effort of two U.S. Department of Energy organizations (Office +// of Science and the National Nuclear Security Administration) responsible for +// the planning and preparation of a capable exascale ecosystem, including +// software, applications, hardware, advanced system engineering and early +// testbed platforms, in support of the nation's exascale computing imperative. + +// Build L2 constant basis + +static void L2BasisP0(CeedInt dim, CeedInt Q, CeedScalar *q_ref, CeedScalar *q_weights, CeedScalar *interp, CeedQuadMode quad_mode) { + // Get 1D quadrature on [-1,1] + CeedScalar q_ref_1d[Q], q_weight_1d[Q]; + switch (quad_mode) { + case CEED_GAUSS: + CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d); + break; + case CEED_GAUSS_LOBATTO: + CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d); + break; + } + + // P0 L2 basis is just a constant + CeedScalar P0 = 1.0; + // Loop over quadrature points + if (dim == 2) { + for (CeedInt i = 0; i < Q; i++) { + for (CeedInt j = 0; j < Q; j++) { + CeedInt k1 = Q * i + j; + q_ref[k1] = q_ref_1d[j]; + q_ref[k1 + Q * Q] = q_ref_1d[i]; + q_weights[k1] = q_weight_1d[j] * q_weight_1d[i]; + interp[k1] = P0; + } + } + } else { + for (CeedInt k = 0; k < Q; k++) { + for (CeedInt i = 0; i < Q; i++) { + for (CeedInt j = 0; j < Q; j++) { + CeedInt k1 = Q * Q * k + Q * i + j; + q_ref[k1 + 0 * Q * Q] = q_ref_1d[j]; + q_ref[k1 + 1 * Q * Q] = q_ref_1d[i]; + q_ref[k1 + 2 * Q * Q] = q_ref_1d[k]; + q_weights[k1] = q_weight_1d[j] * q_weight_1d[i] * q_weight_1d[k]; + interp[k1] = P0; + } + } + } + } +} diff --git a/examples/Hdiv-mass/conv_plot.py b/examples/Hdiv-mass/conv_rate.py similarity index 75% rename from examples/Hdiv-mass/conv_plot.py rename to examples/Hdiv-mass/conv_rate.py index 6e9b9e623d..18aff73f3c 100644 --- a/examples/Hdiv-mass/conv_plot.py +++ b/examples/Hdiv-mass/conv_rate.py @@ -16,13 +16,16 @@ # software, applications, hardware, advanced system engineering and early # testbed platforms, in support of the nation's exascale computing imperative. +# After ./conv_test.sh you can get the table of convergence order by +# python conv_rate.py -f conv_test_result.csv + import pandas as pd import argparse from pylab import * from matplotlib import use -def plot(): +def convergence_rate(): # Define argparse for the input variables parser = argparse.ArgumentParser(description='Get input arguments') parser.add_argument('-f', @@ -40,23 +43,18 @@ def plot(): data = data.sort_values('run') E_u = data['error_u'] - #E_hdiv = data['error_hdiv'] h = 1/data['mesh_res'] - H2 = amin(E_u)* (h/amin(h))**2 # H = C h^2 + N = data['mesh_res'] + conv_u = [] + conv_u.append(0) - ax.loglog(h, E_u, 'o', color='black', label = 'Velocity') - #ax.loglog(h, E_hdiv, '*', color='red', label = 'Velocity in H(div)') - ax.loglog(h, H2, '--', color='black', label='O(h$^2$)') + for i in range(1,len(E_u)): + conv_u.append(log10(E_u[i]/E_u[i-1])/log10(h[i]/h[i-1])) - ax.legend(loc='upper left') - ax.set_xlabel('h') - ax.set_ylabel('L2 Error') - ax.set_title('Convergence by h Refinement') - #xlim(.06, .3) - fig.tight_layout() - plt.savefig('convrate_mass.png', bbox_inches='tight') - + result = {'Number of element/direction':N, 'convergence order of u':conv_u} + df = pd.DataFrame(result) + print(df) if __name__ == "__main__": - plot() + convergence_rate() diff --git a/examples/Hdiv-mass/conv_test.sh b/examples/Hdiv-mass/conv_test.sh index cdc2fddd74..3ba642a064 100755 --- a/examples/Hdiv-mass/conv_test.sh +++ b/examples/Hdiv-mass/conv_test.sh @@ -24,17 +24,17 @@ do d) dim=${OPTARG};; esac done -echo "Running convergence test in ${dim}D for mass problem"; +echo "Running convergence test in ${dim}D for Projection problem in H(div) space"; declare -A run_flags - run_flags[pc_type]=svd + #run_flags[pc_type]=svd if [[ $dim -eq 2 ]]; then run_flags[problem]=mass2d run_flags[dm_plex_dim]=$dim run_flags[dm_plex_box_faces]=2,2 run_flags[dm_plex_box_lower]=0,0 - run_flags[dm_plex_box_upper]=1,0.1 + run_flags[dm_plex_box_upper]=1,1 else run_flags[problem]=mass3d run_flags[dm_plex_dim]=$dim @@ -42,9 +42,9 @@ declare -A run_flags fi declare -A test_flags - test_flags[res_start]=2 - test_flags[res_stride]=1 - test_flags[res_end]=10 + test_flags[res_start]=4 + test_flags[res_stride]=2 + test_flags[res_end]=12 file_name=conv_test_result.csv @@ -64,7 +64,7 @@ for ((res=${test_flags[res_start]}; res<=${test_flags[res_end]}; res+=${test_fla args="$args -$arg ${run_flags[$arg]}" fi done - ./main $args | grep "L2 Error" | awk -v i="$i" -v res="$res" '{ printf "%d,%d,%.5f\n", i, res, $4}' >> $file_name + ./main $args | grep "L2 error of u" | awk -v i="$i" -v res="$res" '{ printf "%d,%d,%e\n", i, res, $6}' >> $file_name i=$((i+1)) done diff --git a/examples/Hdiv-mass/conv_test_result.csv b/examples/Hdiv-mass/conv_test_result.csv deleted file mode 100644 index 296e9812fd..0000000000 --- a/examples/Hdiv-mass/conv_test_result.csv +++ /dev/null @@ -1,10 +0,0 @@ -run,mesh_res,error_u -0,2,0.06228 -1,3,0.02836 -2,4,0.01616 -3,5,0.01043 -4,6,0.00728 -5,7,0.00537 -6,8,0.00413 -7,9,0.00327 -8,10,0.00265 diff --git a/examples/Hdiv-mass/include/cl-options.h b/examples/Hdiv-mass/include/cl-options.h index a48941521e..79f4fb51bc 100644 --- a/examples/Hdiv-mass/include/cl-options.h +++ b/examples/Hdiv-mass/include/cl-options.h @@ -1,12 +1,9 @@ #ifndef cloptions_h #define cloptions_h -#include "../include/structs.h" - -// Register problems to be available on the command line -PetscErrorCode RegisterProblems_Hdiv(AppCtx app_ctx); +#include "structs.h" // Process general command line options -PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx); +PetscErrorCode ProcessCommandLineOptions(AppCtx app_ctx); #endif // cloptions_h diff --git a/examples/Hdiv-mass/include/matops.h b/examples/Hdiv-mass/include/matops.h deleted file mode 100644 index f3f662215b..0000000000 --- a/examples/Hdiv-mass/include/matops.h +++ /dev/null @@ -1,13 +0,0 @@ -#ifndef matops_h -#define matops_h - -#include -#include - -#include "structs.h" - -PetscErrorCode ApplyLocal_Ceed(User user, Vec X, Vec Y); -PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y); -PetscErrorCode ComputeError(User user, Vec X, CeedVector target, CeedScalar *l2_error); - -#endif // matops_h diff --git a/examples/Hdiv-mass/include/petsc-macros.h b/examples/Hdiv-mass/include/petsc-macros.h deleted file mode 100644 index f8c63ebdc4..0000000000 --- a/examples/Hdiv-mass/include/petsc-macros.h +++ /dev/null @@ -1,17 +0,0 @@ -#ifndef petsc_macros -#define petsc_macros - -#if PETSC_VERSION_LT(3, 14, 0) -#define DMPlexGetClosureIndices(a, b, c, d, e, f, g, h, i) DMPlexGetClosureIndices(a, b, c, d, f, g, i) -#define DMPlexRestoreClosureIndices(a, b, c, d, e, f, g, h, i) DMPlexRestoreClosureIndices(a, b, c, d, f, g, i) -#endif - -#if PETSC_VERSION_LT(3, 14, 0) -#define DMAddBoundary(a, b, c, d, e, f, g, h, i, j, k, l, m, n) DMAddBoundary(a, b, c, e, h, i, j, k, f, g, m) -#elif PETSC_VERSION_LT(3, 16, 0) -#define DMAddBoundary(a, b, c, d, e, f, g, h, i, j, k, l, m, n) DMAddBoundary(a, b, c, e, h, i, j, k, l, f, g, m) -#else -#define DMAddBoundary(a, b, c, d, e, f, g, h, i, j, k, l, m, n) DMAddBoundary(a, b, c, d, f, g, h, i, j, k, l, m, n) -#endif - -#endif diff --git a/examples/Hdiv-mass/include/post-processing.h b/examples/Hdiv-mass/include/post-processing.h new file mode 100644 index 0000000000..6338a2f2c4 --- /dev/null +++ b/examples/Hdiv-mass/include/post-processing.h @@ -0,0 +1,12 @@ +#ifndef post_processing_h +#define post_processing_h + +#include +#include + +#include "setup-fe.h" +#include "structs.h" + +PetscErrorCode PrintOutput(DM dm, Ceed ceed, AppCtx app_ctx, KSP ksp, Vec X, CeedScalar l2_error_u); + +#endif // post_processing_h \ No newline at end of file diff --git a/examples/Hdiv-mass/include/problems.h b/examples/Hdiv-mass/include/problems.h deleted file mode 100644 index 19e0499bb6..0000000000 --- a/examples/Hdiv-mass/include/problems.h +++ /dev/null @@ -1,19 +0,0 @@ -#ifndef problems_h -#define problems_h - -#include "../include/structs.h" - -// ----------------------------------------------------------------------------- -// Set up problems function prototype -// ----------------------------------------------------------------------------- -// 1) poisson-quad2d -PetscErrorCode Hdiv_POISSON_MASS2D(ProblemData *problem_data, void *ctx); - -// 2) poisson-hex3d -PetscErrorCode Hdiv_POISSON_MASS3D(ProblemData *problem_data, void *ctx); - -// 3) poisson-prism3d - -// 4) richard - -#endif // problems_h diff --git a/examples/Hdiv-mass/include/register-problem.h b/examples/Hdiv-mass/include/register-problem.h new file mode 100644 index 0000000000..ab15dfbe50 --- /dev/null +++ b/examples/Hdiv-mass/include/register-problem.h @@ -0,0 +1,18 @@ +#ifndef register_problem_h +#define register_problem_h + +#include "structs.h" + +// Register problems to be available on the command line +PetscErrorCode RegisterProblems_Hdiv(AppCtx app_ctx); + +// ----------------------------------------------------------------------------- +// Set up problems function prototype +// ----------------------------------------------------------------------------- +// 1) poisson-quad2d +PetscErrorCode Hdiv_POISSON_MASS2D(ProblemData problem_data, void *ctx); + +// 2) poisson-hex3d +PetscErrorCode Hdiv_POISSON_MASS3D(ProblemData problem_data, void *ctx); + +#endif // register_problem_h diff --git a/examples/Hdiv-mass/include/setup-dm.h b/examples/Hdiv-mass/include/setup-dm.h index 53ded5c4ce..192db3d16d 100644 --- a/examples/Hdiv-mass/include/setup-dm.h +++ b/examples/Hdiv-mass/include/setup-dm.h @@ -6,11 +6,11 @@ #include #include -#include "../include/structs.h" +#include "structs.h" // --------------------------------------------------------------------------- -// Set-up DM +// Create DM // --------------------------------------------------------------------------- -PetscErrorCode CreateDistributedDM(MPI_Comm comm, ProblemData *problem_data, DM *dm); +PetscErrorCode CreateDM(MPI_Comm comm, Ceed ceed, DM *dm); #endif // setupdm_h diff --git a/examples/Hdiv-mass/include/setup-fe.h b/examples/Hdiv-mass/include/setup-fe.h new file mode 100644 index 0000000000..f30914931b --- /dev/null +++ b/examples/Hdiv-mass/include/setup-fe.h @@ -0,0 +1,20 @@ +#ifndef setupfe_h +#define setupfe_h + +#include +#include +#include +#include + +#include "structs.h" + +// --------------------------------------------------------------------------- +// Setup H(div) FE space +// --------------------------------------------------------------------------- +CeedMemType MemTypeP2C(PetscMemType mtype); +PetscErrorCode SetupFEHdiv(AppCtx app_ctx, ProblemData problem_data, DM dm); +CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type); +PetscInt Involute(PetscInt i); +PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr); +PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, CeedInt P, CeedElemRestriction *elem_restr_u, CeedElemRestriction *elem_restr_p); +#endif // setupfe_h diff --git a/examples/Hdiv-mass/include/setup-libceed.h b/examples/Hdiv-mass/include/setup-libceed.h index 7be8d75464..c28c543108 100644 --- a/examples/Hdiv-mass/include/setup-libceed.h +++ b/examples/Hdiv-mass/include/setup-libceed.h @@ -1,19 +1,10 @@ #ifndef setuplibceed_h #define setuplibceed_h -#include "../include/structs.h" +#include "setup-fe.h" +#include "structs.h" -// Convert PETSc MemType to libCEED MemType -CeedMemType MemTypeP2C(PetscMemType mtype); // Destroy libCEED objects PetscErrorCode CeedDataDestroy(CeedData ceed_data); -// Utility function - essential BC dofs are encoded in closure indices as -(i+1) -PetscInt Involute(PetscInt i); -// Utility function to create local CEED restriction from DMPlex -PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr); -// Utility function to create local CEED Oriented restriction from DMPlex -PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, CeedInt P, CeedElemRestriction *elem_restr_oriented); -// Set up libCEED for a given degree -PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData *problem_data, PetscInt U_g_size, PetscInt U_loc_size, CeedData ceed_data, - CeedVector rhs_ceed, CeedVector *target); +PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData problem_data, CeedData ceed_data, CeedVector rhs_ceed); #endif // setuplibceed_h diff --git a/examples/Hdiv-mass/include/setup-matops.h b/examples/Hdiv-mass/include/setup-matops.h new file mode 100644 index 0000000000..4240f28ef5 --- /dev/null +++ b/examples/Hdiv-mass/include/setup-matops.h @@ -0,0 +1,13 @@ +#ifndef setup_matops_h +#define setup_matops_h + +#include +#include + +#include "setup-fe.h" +#include "structs.h" + +PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, OperatorApplyContext op_apply_ctx); +PetscErrorCode ApplyAddLocalCeedOp(Vec X, Vec Y, OperatorApplyContext op_apply_ctx); + +#endif // setup_matops_h \ No newline at end of file diff --git a/examples/Hdiv-mass/include/setup-solvers.h b/examples/Hdiv-mass/include/setup-solvers.h new file mode 100644 index 0000000000..71110cb53b --- /dev/null +++ b/examples/Hdiv-mass/include/setup-solvers.h @@ -0,0 +1,16 @@ +#ifndef setup_solvers_h +#define setup_solvers_h + +#include +#include + +#include "petscvec.h" +#include "structs.h" + +PetscErrorCode SetupResidualOperatorCtx(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_residual); +PetscErrorCode SetupErrorOperatorCtx(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_error_u); +PetscErrorCode ApplyMatOp(Mat A, Vec X, Vec Y); +PetscErrorCode PDESolver(CeedData ceed_data, AppCtx app_ctx, KSP ksp, Vec rhs, Vec *X); +PetscErrorCode ComputeL2Error(Vec X, PetscScalar *l2_error, OperatorApplyContext op_error_ctx); +PetscErrorCode CtxVecDestroy(ProblemData problem_data, AppCtx app_ctx); +#endif // setup_solvers_h \ No newline at end of file diff --git a/examples/Hdiv-mass/include/structs.h b/examples/Hdiv-mass/include/structs.h index 69f65d916b..a39fb829b0 100644 --- a/examples/Hdiv-mass/include/structs.h +++ b/examples/Hdiv-mass/include/structs.h @@ -4,78 +4,48 @@ #include #include +// PETSc operator contexts +typedef struct OperatorApplyContext_ *OperatorApplyContext; +struct OperatorApplyContext_ { + MPI_Comm comm; + DM dm; + Vec X_loc, Y_loc; + CeedVector x_ceed, y_ceed; + CeedOperator op_apply; + Ceed ceed; +}; + // Application context from user command line options typedef struct AppCtx_ *AppCtx; struct AppCtx_ { + char ceed_resource[PETSC_MAX_PATH_LEN]; // libCEED backend + MPI_Comm comm; // libCEED arguments PetscInt degree; PetscInt q_extra; // Problem type arguments - PetscFunctionList problems; - char problem_name[PETSC_MAX_PATH_LEN]; + PetscFunctionList problems; + char problem_name[PETSC_MAX_PATH_LEN]; + OperatorApplyContext ctx_residual, ctx_error_u; }; // libCEED data struct typedef struct CeedData_ *CeedData; struct CeedData_ { - CeedBasis basis_x, basis_u; - CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i; - CeedQFunction qf_residual, qf_error; - CeedOperator op_residual, op_error; - CeedVector x_ceed, y_ceed; - CeedQFunctionContext pq2d_context; -}; - -// 1) poisson-quad2d -#ifndef PHYSICS_POISSONQUAD2D_STRUCT -#define PHYSICS_POISSONQUAD2D_STRUCT -typedef struct PQ2DContext_ *PQ2DContext; -struct PQ2DContext_ { - CeedScalar kappa; -}; -#endif - -// 2) poisson-hex3d -#ifndef PHYSICS_POISSONHEX3D_STRUCT -#define PHYSICS_POISSONHEX3D_STRUCT -typedef struct PH3DContext_ *PH3DContext; -struct PH3DContext_ { - CeedScalar kappa; -}; -#endif - -// 3) poisson-prism3d - -// 4) richard - -// Struct that contains all enums and structs used for the physics of all problems -typedef struct Physics_ *Physics; -struct Physics_ { - PQ2DContext pq2d_ctx; - PH3DContext ph3d_ctx; -}; - -// PETSc user data -typedef struct User_ *User; -struct User_ { - MPI_Comm comm; - Vec X_loc, Y_loc; - CeedVector x_ceed, y_ceed; - CeedOperator op_apply, op_error; - DM dm; - Ceed ceed; - AppCtx app_ctx; - Physics phys; + CeedBasis basis_x, basis_u, basis_p; + CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_p; + CeedQFunction qf_residual, qf_error; + CeedOperator op_residual, op_error; + CeedVector x_ceed, y_ceed; }; // Problem specific data -typedef struct { +typedef struct ProblemData_ *ProblemData; +struct ProblemData_ { CeedQFunctionUser setup_rhs, residual, setup_error; const char *setup_rhs_loc, *residual_loc, *setup_error_loc; CeedQuadMode quadrature_mode; CeedInt elem_node, dim; - PetscErrorCode (*setup_ctx)(Ceed, CeedData, Physics); - -} ProblemData; +}; #endif // structs_h \ No newline at end of file diff --git a/examples/Hdiv-mass/main.c b/examples/Hdiv-mass/main.c index 338e745cfb..abc929eef8 100644 --- a/examples/Hdiv-mass/main.c +++ b/examples/Hdiv-mass/main.c @@ -14,10 +14,7 @@ // software, applications, hardware, advanced system engineering and early // testbed platforms, in support of the nation's exascale computing imperative. -// libCEED + PETSc Example: Mixed-Poisson in H(div) space -// -// This example demonstrates a simple usage of libCEED with PETSc to solve -// elasticity problems. +// libCEED + PETSc Example: Projection problem in H(div) space // // The code uses higher level communication protocols in DMPlex. // @@ -27,7 +24,8 @@ // ./main -pc_type svd -problem mass2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 // ./main -pc_type svd -problem mass3d -dm_plex_dim 3 -dm_plex_box_faces 4,4,4 -const char help[] = "Solve H(div)-mixed problem using PETSc and libCEED\n"; +#include +const char help[] = "Solve Projection problem in H(div) space using PETSc and libCEED\n"; #include "main.h" @@ -35,59 +33,34 @@ int main(int argc, char **argv) { // --------------------------------------------------------------------------- // Initialize PETSc // --------------------------------------------------------------------------- - PetscInt ierr; - ierr = PetscInitialize(&argc, &argv, NULL, help); - if (ierr) return ierr; + PetscCall(PetscInitialize(&argc, &argv, NULL, help)); + MPI_Comm comm = PETSC_COMM_WORLD; // --------------------------------------------------------------------------- // Create structs // --------------------------------------------------------------------------- AppCtx app_ctx; - ierr = PetscCalloc1(1, &app_ctx); - CHKERRQ(ierr); + PetscCall(PetscCalloc1(1, &app_ctx)); - ProblemData *problem_data = NULL; - ierr = PetscCalloc1(1, &problem_data); - CHKERRQ(ierr); - - User user; - ierr = PetscCalloc1(1, &user); - CHKERRQ(ierr); + ProblemData problem_data = NULL; + PetscCall(PetscCalloc1(1, &problem_data)); CeedData ceed_data; - ierr = PetscCalloc1(1, &ceed_data); - CHKERRQ(ierr); - - Physics phys_ctx; - ierr = PetscCalloc1(1, &phys_ctx); - CHKERRQ(ierr); + PetscCall(PetscCalloc1(1, &ceed_data)); - user->app_ctx = app_ctx; - user->phys = phys_ctx; + OperatorApplyContext ctx_residual, ctx_error_u; + PetscCall(PetscCalloc1(1, &ctx_residual)); + PetscCall(PetscCalloc1(1, &ctx_error_u)); + // Context for residual + app_ctx->ctx_residual = ctx_residual; + // Context for computing error + app_ctx->ctx_error_u = ctx_error_u; + app_ctx->comm = comm; // --------------------------------------------------------------------------- // Process command line options // --------------------------------------------------------------------------- - // -- Register problems to be available on the command line - ierr = RegisterProblems_Hdiv(app_ctx); - CHKERRQ(ierr); - - // -- Process general command line options - MPI_Comm comm = PETSC_COMM_WORLD; - ierr = ProcessCommandLineOptions(comm, app_ctx); - CHKERRQ(ierr); - - // --------------------------------------------------------------------------- - // Choose the problem from the list of registered problems - // --------------------------------------------------------------------------- - { - PetscErrorCode (*p)(ProblemData *, void *); - ierr = PetscFunctionListFind(app_ctx->problems, app_ctx->problem_name, &p); - CHKERRQ(ierr); - if (!p) SETERRQ(PETSC_COMM_SELF, 1, "Problem '%s' not found", app_ctx->problem_name); - ierr = (*p)(problem_data, &user); - CHKERRQ(ierr); - } + PetscCall(ProcessCommandLineOptions(app_ctx)); // --------------------------------------------------------------------------- // Initialize libCEED @@ -95,202 +68,104 @@ int main(int argc, char **argv) { // -- Initialize backend Ceed ceed; CeedInit("/cpu/self/ref/serial", &ceed); - CeedMemType mem_type_backend; - CeedGetPreferredMemType(ceed, &mem_type_backend); + // CeedInit(app_ctx->ceed_resource, &ceed); + // --------------------------------------------------------------------------- - // Set-up DM + // Choose the problem from the list of registered problems // --------------------------------------------------------------------------- - // PETSc objects - DM dm; - VecType vec_type; - ierr = CreateDistributedDM(comm, problem_data, &dm); - CHKERRQ(ierr); - ierr = DMGetVecType(dm, &vec_type); - CHKERRQ(ierr); - if (!vec_type) { // Not yet set by user -dm_vec_type - switch (mem_type_backend) { - case CEED_MEM_HOST: - vec_type = VECSTANDARD; - break; - case CEED_MEM_DEVICE: { - const char *resolved; - CeedGetResource(ceed, &resolved); - if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA; - else if (strstr(resolved, "/gpu/hip/occa")) vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678 - else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP; - else vec_type = VECSTANDARD; - } - } - ierr = DMSetVecType(dm, vec_type); - CHKERRQ(ierr); + PetscCall(RegisterProblems_Hdiv(app_ctx)); + { + PetscErrorCode (*p)(ProblemData, void *); + PetscCall(PetscFunctionListFind(app_ctx->problems, app_ctx->problem_name, &p)); + if (!p) SETERRQ(PETSC_COMM_SELF, 1, "Problem '%s' not found", app_ctx->problem_name); + PetscCall((*p)(problem_data, &app_ctx)); } + // --------------------------------------------------------------------------- - // Create global, local solution, local rhs vector + // Create DM and Setup FE space // --------------------------------------------------------------------------- - Vec U_g, U_loc; - PetscInt U_l_size, U_g_size, U_loc_size; - // Create global and local solution vectors - ierr = DMCreateGlobalVector(dm, &U_g); - CHKERRQ(ierr); - ierr = VecGetSize(U_g, &U_g_size); - CHKERRQ(ierr); - // Local size for matShell - ierr = VecGetLocalSize(U_g, &U_l_size); - CHKERRQ(ierr); - // Create local unknown vector U_loc - ierr = DMCreateLocalVector(dm, &U_loc); - CHKERRQ(ierr); - // Local size for libCEED - ierr = VecGetSize(U_loc, &U_loc_size); - CHKERRQ(ierr); + DM dm; + PetscCall(CreateDM(comm, ceed, &dm)); + PetscCall(SetupFEHdiv(app_ctx, problem_data, dm)); - // Get RHS vector + // --------------------------------------------------------------------------- + // Create zero rhs local vector + // --------------------------------------------------------------------------- + CeedVector rhs_ceed; Vec rhs_loc; PetscScalar *r; - CeedVector rhs_ceed, target; PetscMemType mem_type; - ierr = VecDuplicate(U_loc, &rhs_loc); - CHKERRQ(ierr); - ierr = VecZeroEntries(rhs_loc); - CHKERRQ(ierr); - ierr = VecGetArrayAndMemType(rhs_loc, &r, &mem_type); - CHKERRQ(ierr); - CeedVectorCreate(ceed, U_l_size, &rhs_ceed); + PetscInt rhs_l_size; + // Create global and local solution vectors + PetscCall(DMCreateLocalVector(dm, &rhs_loc)); + PetscCall(VecGetSize(rhs_loc, &rhs_l_size)); + PetscCall(VecZeroEntries(rhs_loc)); + PetscCall(VecGetArrayAndMemType(rhs_loc, &r, &mem_type)); + CeedVectorCreate(ceed, rhs_l_size, &rhs_ceed); CeedVectorSetArray(rhs_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, r); // --------------------------------------------------------------------------- - // Setup libCEED + // Setup libCEED qfunctions and operators // --------------------------------------------------------------------------- - // -- Set up libCEED objects - ierr = SetupLibceed(dm, ceed, app_ctx, problem_data, U_g_size, U_loc_size, ceed_data, rhs_ceed, &target); - CHKERRQ(ierr); + PetscCall(SetupLibceed(dm, ceed, app_ctx, problem_data, ceed_data, rhs_ceed)); // --------------------------------------------------------------------------- - // Gather RHS + // Setup rhs global vector entries with the computed rhs_ceed // --------------------------------------------------------------------------- Vec rhs; + PetscCall(DMCreateGlobalVector(dm, &rhs)); CeedVectorTakeArray(rhs_ceed, MemTypeP2C(mem_type), NULL); - ierr = VecRestoreArrayAndMemType(rhs_loc, &r); - CHKERRQ(ierr); - ierr = VecDuplicate(U_g, &rhs); - CHKERRQ(ierr); - ierr = VecZeroEntries(rhs); - CHKERRQ(ierr); - ierr = DMLocalToGlobal(dm, rhs_loc, ADD_VALUES, rhs); - CHKERRQ(ierr); - // VecView(rhs, PETSC_VIEWER_STDOUT_WORLD); + PetscCall(VecRestoreArrayAndMemType(rhs_loc, &r)); + PetscCall(VecZeroEntries(rhs)); + PetscCall(DMLocalToGlobal(dm, rhs_loc, ADD_VALUES, rhs)); + CeedVectorDestroy(&rhs_ceed); // --------------------------------------------------------------------------- - // Setup Mat, KSP + // Solve A*X=rhs; setup-solver.c // --------------------------------------------------------------------------- - user->comm = comm; - user->dm = dm; - user->X_loc = U_loc; - ierr = VecDuplicate(U_loc, &user->Y_loc); - CHKERRQ(ierr); - user->x_ceed = ceed_data->x_ceed; - user->y_ceed = ceed_data->y_ceed; - user->op_apply = ceed_data->op_residual; - user->op_error = ceed_data->op_error; - user->ceed = ceed; - // Operator - Mat mat; - ierr = MatCreateShell(comm, U_l_size, U_l_size, U_g_size, U_g_size, user, &mat); - CHKERRQ(ierr); - ierr = MatShellSetOperation(mat, MATOP_MULT, (void (*)(void))MatMult_Ceed); - CHKERRQ(ierr); - ierr = MatShellSetVecType(mat, vec_type); - CHKERRQ(ierr); - + Vec X; KSP ksp; - ierr = KSPCreate(comm, &ksp); - CHKERRQ(ierr); - ierr = KSPSetOperators(ksp, mat, mat); - CHKERRQ(ierr); - ierr = KSPSetFromOptions(ksp); - CHKERRQ(ierr); - ierr = KSPSetUp(ksp); - CHKERRQ(ierr); - ierr = KSPSolve(ksp, rhs, U_g); - CHKERRQ(ierr); - // printf("U_g\n"); - // VecView(U_g, PETSC_VIEWER_STDOUT_WORLD); - // --------------------------------------------------------------------------- - // Compute pointwise L2 maximum error - // --------------------------------------------------------------------------- - CeedScalar l2_error; - ierr = ComputeError(user, U_g, target, &l2_error); - CHKERRQ(ierr); + PetscCall(SetupResidualOperatorCtx(app_ctx->comm, dm, ceed, ceed_data, app_ctx->ctx_residual)); + // Create global and local solution vectors + PetscCall(DMCreateGlobalVector(dm, &X)); + PetscCall(KSPCreate(app_ctx->comm, &ksp)); + PetscCall(KSPSetDM(ksp, dm)); + PetscCall(KSPSetDMActive(ksp, PETSC_FALSE)); + PetscCall(PDESolver(ceed_data, app_ctx, ksp, rhs, &X)); + + // --------------------------------------------------------------------------- + // Compute L2 error of mms problem; setup-solver.c + // --------------------------------------------------------------------------- + CeedScalar l2_error_u = 0.0; + PetscCall(SetupErrorOperatorCtx(app_ctx->comm, dm, ceed, ceed_data, app_ctx->ctx_error_u)); + PetscCall(ComputeL2Error(X, &l2_error_u, app_ctx->ctx_error_u)); // --------------------------------------------------------------------------- - // Output results + // Print solver iterations and final norms; post-processing // --------------------------------------------------------------------------- - KSPType ksp_type; - KSPConvergedReason reason; - PetscReal rnorm; - PetscInt its; - ierr = KSPGetType(ksp, &ksp_type); - CHKERRQ(ierr); - ierr = KSPGetConvergedReason(ksp, &reason); - CHKERRQ(ierr); - ierr = KSPGetIterationNumber(ksp, &its); - CHKERRQ(ierr); - ierr = KSPGetResidualNorm(ksp, &rnorm); - CHKERRQ(ierr); - ierr = PetscPrintf(comm, - " KSP:\n" - " KSP Type : %s\n" - " KSP Convergence : %s\n" - " Total KSP Iterations : %" PetscInt_FMT - "\n" - " Final rnorm : %e\n" - " L2 Error : %e\n", - ksp_type, KSPConvergedReasons[reason], its, (double)rnorm, (double)l2_error); - CHKERRQ(ierr); + PetscCall(PrintOutput(dm, ceed, app_ctx, ksp, X, l2_error_u)); // --------------------------------------------------------------------------- // Free objects // --------------------------------------------------------------------------- // Free PETSc objects - ierr = DMDestroy(&dm); - CHKERRQ(ierr); - ierr = VecDestroy(&U_g); - CHKERRQ(ierr); - ierr = VecDestroy(&U_loc); - CHKERRQ(ierr); - ierr = VecDestroy(&rhs); - CHKERRQ(ierr); - ierr = VecDestroy(&rhs_loc); - CHKERRQ(ierr); - ierr = VecDestroy(&user->Y_loc); - CHKERRQ(ierr); - ierr = MatDestroy(&mat); - CHKERRQ(ierr); - ierr = KSPDestroy(&ksp); - CHKERRQ(ierr); - + PetscCall(DMDestroy(&dm)); + PetscCall(VecDestroy(&X)); + PetscCall(VecDestroy(&rhs)); + PetscCall(VecDestroy(&rhs_loc)); + PetscCall(KSPDestroy(&ksp)); + PetscCall(CtxVecDestroy(problem_data, app_ctx)); // -- Function list - ierr = PetscFunctionListDestroy(&app_ctx->problems); - CHKERRQ(ierr); - + PetscCall(PetscFunctionListDestroy(&app_ctx->problems)); // -- Structs - ierr = PetscFree(app_ctx); - CHKERRQ(ierr); - ierr = PetscFree(problem_data); - CHKERRQ(ierr); - ierr = PetscFree(user); - CHKERRQ(ierr); - ierr = PetscFree(phys_ctx->pq2d_ctx); - CHKERRQ(ierr); - ierr = PetscFree(phys_ctx); - CHKERRQ(ierr); - + PetscCall(CeedDataDestroy(ceed_data)); + PetscCall(PetscFree(app_ctx)); + PetscCall(PetscFree(ctx_residual)); + PetscCall(PetscFree(ctx_error_u)); + PetscCall(PetscFree(problem_data)); // Free libCEED objects CeedVectorDestroy(&rhs_ceed); - CeedVectorDestroy(&target); - ierr = CeedDataDestroy(ceed_data); - CHKERRQ(ierr); CeedDestroy(&ceed); return PetscFinalize(); diff --git a/examples/Hdiv-mass/main.h b/examples/Hdiv-mass/main.h index c5ecde00d9..f433f47a27 100644 --- a/examples/Hdiv-mass/main.h +++ b/examples/Hdiv-mass/main.h @@ -3,9 +3,12 @@ #define MAIN_H #include "include/cl-options.h" -#include "include/matops.h" -#include "include/problems.h" +#include "include/post-processing.h" +#include "include/register-problem.h" #include "include/setup-dm.h" +#include "include/setup-fe.h" #include "include/setup-libceed.h" +#include "include/setup-matops.h" +#include "include/setup-solvers.h" #endif // MAIN_H \ No newline at end of file diff --git a/examples/Hdiv-mass/problems/mass2d.c b/examples/Hdiv-mass/problems/mass2d.c index 2e58f87045..0cdef98f85 100644 --- a/examples/Hdiv-mass/problems/mass2d.c +++ b/examples/Hdiv-mass/problems/mass2d.c @@ -17,22 +17,16 @@ /// @file /// Utility functions for setting up POISSON_QUAD2D -#include "../include/problems.h" +#include "../include/register-problem.h" #include "../include/setup-libceed.h" #include "../qfunctions/poisson-error2d.h" #include "../qfunctions/poisson-mass2d.h" #include "../qfunctions/poisson-rhs2d.h" // Hdiv_POISSON_MASS2D is registered in cl-option.c -PetscErrorCode Hdiv_POISSON_MASS2D(ProblemData *problem_data, void *ctx) { - User user = *(User *)ctx; - MPI_Comm comm = PETSC_COMM_WORLD; - PetscInt ierr; +PetscErrorCode Hdiv_POISSON_MASS2D(ProblemData problem_data, void *ctx) { PetscFunctionBeginUser; - ierr = PetscCalloc1(1, &user->phys->pq2d_ctx); - CHKERRQ(ierr); - // ------------------------------------------------------ // SET UP POISSON_QUAD2D // ------------------------------------------------------ @@ -46,12 +40,5 @@ PetscErrorCode Hdiv_POISSON_MASS2D(ProblemData *problem_data, void *ctx) { problem_data->setup_error = SetupError2D; problem_data->setup_error_loc = SetupError2D_loc; - // ------------------------------------------------------ - // Command line Options - // ------------------------------------------------------ - PetscOptionsBegin(comm, NULL, "Options for Hdiv-mass problem", NULL); - - PetscOptionsEnd(); - PetscFunctionReturn(0); } diff --git a/examples/Hdiv-mass/problems/mass3d.c b/examples/Hdiv-mass/problems/mass3d.c index 2997ebb572..b2b847d19f 100644 --- a/examples/Hdiv-mass/problems/mass3d.c +++ b/examples/Hdiv-mass/problems/mass3d.c @@ -15,24 +15,16 @@ // testbed platforms, in support of the nation's exascale computing imperative. /// @file -/// Utility functions for setting up POISSON_QUAD2D -#include "../include/problems.h" +#include "../include/register-problem.h" #include "../include/setup-libceed.h" #include "../qfunctions/poisson-error3d.h" #include "../qfunctions/poisson-mass3d.h" #include "../qfunctions/poisson-rhs3d.h" -// Hdiv_POISSON_MASS2D is registered in cl-option.c -PetscErrorCode Hdiv_POISSON_MASS3D(ProblemData *problem_data, void *ctx) { - User user = *(User *)ctx; - MPI_Comm comm = PETSC_COMM_WORLD; - PetscInt ierr; +PetscErrorCode Hdiv_POISSON_MASS3D(ProblemData problem_data, void *ctx) { PetscFunctionBeginUser; - ierr = PetscCalloc1(1, &user->phys->ph3d_ctx); - CHKERRQ(ierr); - // ------------------------------------------------------ // SET UP POISSON_QUAD3D // ------------------------------------------------------ @@ -46,12 +38,5 @@ PetscErrorCode Hdiv_POISSON_MASS3D(ProblemData *problem_data, void *ctx) { problem_data->setup_error = SetupError3D; problem_data->setup_error_loc = SetupError3D_loc; - // ------------------------------------------------------ - // Command line Options - // ------------------------------------------------------ - PetscOptionsBegin(comm, NULL, "Options for Hdiv-mass problem", NULL); - - PetscOptionsEnd(); - PetscFunctionReturn(0); } diff --git a/examples/Hdiv-mass/problems/register-problem.c b/examples/Hdiv-mass/problems/register-problem.c new file mode 100644 index 0000000000..0e66ad3a15 --- /dev/null +++ b/examples/Hdiv-mass/problems/register-problem.c @@ -0,0 +1,33 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at +// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights +// reserved. See files LICENSE and NOTICE for details. +// +// This file is part of CEED, a collection of benchmarks, miniapps, software +// libraries and APIs for efficient high-order finite element and spectral +// element discretizations for exascale applications. For more information and +// source code availability see http://github.com/ceed. +// +// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, +// a collaborative effort of two U.S. Department of Energy organizations (Office +// of Science and the National Nuclear Security Administration) responsible for +// the planning and preparation of a capable exascale ecosystem, including +// software, applications, hardware, advanced system engineering and early +// testbed platforms, in support of the nation's exascale computing imperative. + +/// @file +/// Command line option processing for H(div) example using PETSc + +#include "../include/register-problem.h" + +// Register problems to be available on the command line +PetscErrorCode RegisterProblems_Hdiv(AppCtx app_ctx) { + app_ctx->problems = NULL; + + PetscFunctionBeginUser; + // 1) poisson-quad2d (Hdiv_POISSON_MASS2D is created in poisson-mass2d.c) + PetscCall(PetscFunctionListAdd(&app_ctx->problems, "mass2d", Hdiv_POISSON_MASS2D)); + // 2) poisson-hex3d + PetscCall(PetscFunctionListAdd(&app_ctx->problems, "mass3d", Hdiv_POISSON_MASS3D)); + + PetscFunctionReturn(0); +} diff --git a/examples/Hdiv-mass/qfunctions/poisson-error2d.h b/examples/Hdiv-mass/qfunctions/poisson-error2d.h index f48574f93c..e7f19719dc 100644 --- a/examples/Hdiv-mass/qfunctions/poisson-error2d.h +++ b/examples/Hdiv-mass/qfunctions/poisson-error2d.h @@ -45,9 +45,10 @@ CEED_QFUNCTION(SetupError2D)(void *ctx, const CeedInt Q, const CeedScalar *const CeedScalar u1[2] = {u[0][i], u[1][i]}, uh[2]; AlphaMatVecMult2x2(1 / det_J, J, u1, uh); // Error - error[i + 0 * Q] = (uh[0] - target[i + 0 * Q]) * (uh[0] - target[i + 0 * Q]) * w[i] * det_J; - error[i + 1 * Q] = (uh[1] - target[i + 1 * Q]) * (uh[1] - target[i + 1 * Q]) * w[i] * det_J; + CeedScalar err2_ux = (uh[0] - target[i + 0 * Q]) * (uh[0] - target[i + 0 * Q]); + CeedScalar err2_uy = (uh[1] - target[i + 1 * Q]) * (uh[1] - target[i + 1 * Q]); + error[i + 0 * Q] = (err2_ux + err2_uy) * w[i] * det_J; } // End of Quadrature Point Loop return 0; diff --git a/examples/Hdiv-mass/qfunctions/poisson-error3d.h b/examples/Hdiv-mass/qfunctions/poisson-error3d.h index 2d1c821165..b8b09d5025 100644 --- a/examples/Hdiv-mass/qfunctions/poisson-error3d.h +++ b/examples/Hdiv-mass/qfunctions/poisson-error3d.h @@ -46,10 +46,11 @@ CEED_QFUNCTION(SetupError3D)(void *ctx, const CeedInt Q, const CeedScalar *const CeedScalar u1[3] = {u[0][i], u[1][i], u[2][i]}, uh[3]; AlphaMatVecMult3x3(1 / det_J, J, u1, uh); // Error - error[i + 0 * Q] = (uh[0] - target[i + 0 * Q]) * (uh[0] - target[i + 0 * Q]) * w[i] * det_J; - error[i + 1 * Q] = (uh[1] - target[i + 1 * Q]) * (uh[1] - target[i + 1 * Q]) * w[i] * det_J; - error[i + 2 * Q] = (uh[2] - target[i + 2 * Q]) * (uh[2] - target[i + 2 * Q]) * w[i] * det_J; + CeedScalar err2_ux = (uh[0] - target[i + 0 * Q]) * (uh[0] - target[i + 0 * Q]); + CeedScalar err2_uy = (uh[1] - target[i + 1 * Q]) * (uh[1] - target[i + 1 * Q]); + CeedScalar err2_uz = (uh[2] - target[i + 2 * Q]) * (uh[2] - target[i + 2 * Q]); + error[i + 0 * Q] = (err2_ux + err2_uy + err2_uz) * w[i] * det_J; } // End of Quadrature Point Loop return 0; diff --git a/examples/Hdiv-mass/src/cl-options.c b/examples/Hdiv-mass/src/cl-options.c index 105de000df..4b785b430a 100644 --- a/examples/Hdiv-mass/src/cl-options.c +++ b/examples/Hdiv-mass/src/cl-options.c @@ -19,46 +19,31 @@ #include "../include/cl-options.h" -#include "../include/problems.h" - -// Register problems to be available on the command line -PetscErrorCode RegisterProblems_Hdiv(AppCtx app_ctx) { - app_ctx->problems = NULL; - PetscErrorCode ierr; +// Process general command line options +PetscErrorCode ProcessCommandLineOptions(AppCtx app_ctx) { + PetscBool problem_flag = PETSC_FALSE; + PetscBool ceed_flag = PETSC_FALSE; PetscFunctionBeginUser; - // 1) poisson-quad2d (Hdiv_POISSON_MASS2D is created in poisson-mass2d.c) - ierr = PetscFunctionListAdd(&app_ctx->problems, "mass2d", Hdiv_POISSON_MASS2D); - CHKERRQ(ierr); - // 2) poisson-hex3d - ierr = PetscFunctionListAdd(&app_ctx->problems, "mass3d", Hdiv_POISSON_MASS3D); - CHKERRQ(ierr); - - // 3) poisson-prism3d - // 4) richard + PetscOptionsBegin(app_ctx->comm, NULL, "H(div) examples in PETSc with libCEED", NULL); - PetscFunctionReturn(0); -} + PetscCall(PetscOptionsString("-ceed", "CEED resource specifier", NULL, app_ctx->ceed_resource, app_ctx->ceed_resource, + sizeof(app_ctx->ceed_resource), &ceed_flag)); -// Process general command line options -PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx) { - PetscBool problem_flag = PETSC_FALSE; - PetscErrorCode ierr; - PetscFunctionBeginUser; - - PetscOptionsBegin(comm, NULL, "H(div) examples in PETSc with libCEED", NULL); + // Provide default ceed resource if not specified + if (!ceed_flag) { + const char *ceed_resource = "/cpu/self"; + strncpy(app_ctx->ceed_resource, ceed_resource, 10); + } - ierr = PetscOptionsFList("-problem", "Problem to solve", NULL, app_ctx->problems, app_ctx->problem_name, app_ctx->problem_name, - sizeof(app_ctx->problem_name), &problem_flag); - CHKERRQ(ierr); + PetscCall(PetscOptionsFList("-problem", "Problem to solve", NULL, app_ctx->problems, app_ctx->problem_name, app_ctx->problem_name, + sizeof(app_ctx->problem_name), &problem_flag)); app_ctx->degree = 1; - ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", NULL, app_ctx->degree, &app_ctx->degree, NULL); - CHKERRQ(ierr); + PetscCall(PetscOptionsInt("-degree", "Polynomial degree of finite elements", NULL, app_ctx->degree, &app_ctx->degree, NULL)); app_ctx->q_extra = 0; - ierr = PetscOptionsInt("-q_extra", "Number of extra quadrature points", NULL, app_ctx->q_extra, &app_ctx->q_extra, NULL); - CHKERRQ(ierr); + PetscCall(PetscOptionsInt("-q_extra", "Number of extra quadrature points", NULL, app_ctx->q_extra, &app_ctx->q_extra, NULL)); // Provide default problem if not specified if (!problem_flag) { diff --git a/examples/Hdiv-mass/src/matops.c b/examples/Hdiv-mass/src/matops.c deleted file mode 100644 index 45c39f0160..0000000000 --- a/examples/Hdiv-mass/src/matops.c +++ /dev/null @@ -1,106 +0,0 @@ -#include "../include/matops.h" - -#include "../include/setup-libceed.h" - -// ----------------------------------------------------------------------------- -// This function uses libCEED to compute the action of the Laplacian with -// Dirichlet boundary conditions -// ----------------------------------------------------------------------------- -PetscErrorCode ApplyLocal_Ceed(User user, Vec X, Vec Y) { - PetscErrorCode ierr; - PetscScalar *x, *y; - PetscMemType x_mem_type, y_mem_type; - - PetscFunctionBeginUser; - - // Global-to-local - ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc); - CHKERRQ(ierr); - - // Setup libCEED vectors - ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, &x_mem_type); - CHKERRQ(ierr); - ierr = VecGetArrayAndMemType(user->Y_loc, &y, &y_mem_type); - CHKERRQ(ierr); - CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); - CeedVectorSetArray(user->y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); - - // Apply libCEED operator - CeedOperatorApply(user->op_apply, user->x_ceed, user->y_ceed, CEED_REQUEST_IMMEDIATE); - - // Restore PETSc vectors - CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); - CeedVectorTakeArray(user->y_ceed, MemTypeP2C(y_mem_type), NULL); - ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x); - CHKERRQ(ierr); - ierr = VecRestoreArrayAndMemType(user->Y_loc, &y); - CHKERRQ(ierr); - - // Local-to-global - ierr = VecZeroEntries(Y); - CHKERRQ(ierr); - ierr = DMLocalToGlobal(user->dm, user->Y_loc, ADD_VALUES, Y); - CHKERRQ(ierr); - - PetscFunctionReturn(0); -}; - -// ----------------------------------------------------------------------------- -// This function wraps the libCEED operator for a MatShell -// ----------------------------------------------------------------------------- -PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y) { - PetscErrorCode ierr; - User user; - - PetscFunctionBeginUser; - - ierr = MatShellGetContext(A, &user); - CHKERRQ(ierr); - - // libCEED for local action of residual evaluator - ierr = ApplyLocal_Ceed(user, X, Y); - CHKERRQ(ierr); - - PetscFunctionReturn(0); -}; - -// ----------------------------------------------------------------------------- -// This function calculates the error in the final solution -// ----------------------------------------------------------------------------- -PetscErrorCode ComputeError(User user, Vec X, CeedVector target, CeedScalar *l2_error) { - PetscErrorCode ierr; - PetscScalar *x; - PetscMemType mem_type; - CeedVector collocated_error; - CeedSize length; - - PetscFunctionBeginUser; - CeedVectorGetLength(target, &length); - CeedVectorCreate(user->ceed, length, &collocated_error); - - // Global-to-local - ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc); - CHKERRQ(ierr); - - // Setup CEED vector - ierr = VecGetArrayAndMemType(user->X_loc, &x, &mem_type); - CHKERRQ(ierr); - CeedVectorSetArray(user->x_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, x); - - // Apply CEED operator - CeedOperatorApply(user->op_error, user->x_ceed, collocated_error, CEED_REQUEST_IMMEDIATE); - // Restore PETSc vector - CeedVectorTakeArray(user->x_ceed, MemTypeP2C(mem_type), NULL); - ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x); - CHKERRQ(ierr); - - CeedScalar error; - CeedVectorNorm(collocated_error, CEED_NORM_1, &error); - *l2_error = sqrt(error); - // Cleanup - CeedVectorDestroy(&collocated_error); - - PetscFunctionReturn(0); -}; - -// ----------------------------------------------------------------------------- diff --git a/examples/Hdiv-mass/src/post-processing.c b/examples/Hdiv-mass/src/post-processing.c new file mode 100644 index 0000000000..c1281a1146 --- /dev/null +++ b/examples/Hdiv-mass/src/post-processing.c @@ -0,0 +1,77 @@ +#include "../include/post-processing.h" + +// ----------------------------------------------------------------------------- +// This function print the output +// ----------------------------------------------------------------------------- +PetscErrorCode PrintOutput(DM dm, Ceed ceed, AppCtx app_ctx, KSP ksp, Vec X, CeedScalar l2_error_u) { + PetscFunctionBeginUser; + + const char *used_resource; + CeedMemType mem_type_backend; + CeedGetResource(ceed, &used_resource); + CeedGetPreferredMemType(ceed, &mem_type_backend); + char hostname[PETSC_MAX_PATH_LEN]; + PetscCall(PetscGetHostName(hostname, sizeof hostname)); + PetscMPIInt comm_size; + PetscCall(MPI_Comm_size(app_ctx->comm, &comm_size)); + PetscCall(PetscPrintf(app_ctx->comm, + "\n-- Mixed-Elasticity Example - libCEED + PETSc --\n" + " MPI:\n" + " Hostname : %s\n" + " Total ranks : %d\n" + " libCEED:\n" + " libCEED Backend : %s\n" + " libCEED Backend MemType : %s\n", + hostname, comm_size, used_resource, CeedMemTypes[mem_type_backend])); + + MatType mat_type; + VecType vec_type; + PetscCall(DMGetMatType(dm, &mat_type)); + PetscCall(DMGetVecType(dm, &vec_type)); + PetscCall(PetscPrintf(app_ctx->comm, + " PETSc:\n" + " DM MatType : %s\n" + " DM VecType : %s\n", + mat_type, vec_type)); + PetscInt X_l_size, X_g_size; + PetscCall(VecGetSize(X, &X_g_size)); + PetscCall(VecGetLocalSize(X, &X_l_size)); + PetscInt c_start, c_end; + PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); + DMPolytopeType cell_type; + PetscCall(DMPlexGetCellType(dm, c_start, &cell_type)); + CeedElemTopology elem_topo = ElemTopologyP2C(cell_type); + PetscCall(PetscPrintf(app_ctx->comm, + " Problem:\n" + " Problem Name : %s\n" + " Mesh:\n" + " Local Elements : %" PetscInt_FMT "\n" + " Element topology : %s\n", + app_ctx->problem_name, c_end - c_start, CeedElemTopologies[elem_topo])); + PetscInt ksp_its; + KSPType ksp_type; + KSPConvergedReason ksp_reason; + PetscReal ksp_rnorm; + PC pc; + PCType pc_type; + PetscCall(KSPGetPC(ksp, &pc)); + PetscCall(PCGetType(pc, &pc_type)); + PetscCall(KSPGetType(ksp, &ksp_type)); + PetscCall(KSPGetConvergedReason(ksp, &ksp_reason)); + PetscCall(KSPGetIterationNumber(ksp, &ksp_its)); + PetscCall(KSPGetResidualNorm(ksp, &ksp_rnorm)); + PetscCall(PetscPrintf(app_ctx->comm, + " KSP:\n" + " KSP Type : %s\n" + " PC Type : %s\n" + " KSP Convergence : %s\n" + " Total KSP Iterations : %" PetscInt_FMT "\n" + " Final rnorm : %e\n", + ksp_type, pc_type, KSPConvergedReasons[ksp_reason], ksp_its, (double)ksp_rnorm)); + + PetscCall(PetscPrintf(app_ctx->comm, + " L2 Error (MMS):\n" + " L2 error of u : %e\n", + (double)l2_error_u)); + PetscFunctionReturn(0); +}; \ No newline at end of file diff --git a/examples/Hdiv-mass/src/setup-dm.c b/examples/Hdiv-mass/src/setup-dm.c index a48ee93606..ce925c71a8 100644 --- a/examples/Hdiv-mass/src/setup-dm.c +++ b/examples/Hdiv-mass/src/setup-dm.c @@ -1,68 +1,45 @@ #include "../include/setup-dm.h" +#include "petscerror.h" + // --------------------------------------------------------------------------- -// Set-up DM +// Create DM // --------------------------------------------------------------------------- -PetscErrorCode CreateDistributedDM(MPI_Comm comm, ProblemData *problem_data, DM *dm) { - PetscErrorCode ierr; - PetscSection sec; - PetscInt dofs_per_face; - PetscInt p_start, p_end; - PetscInt c_start, c_end; // cells - PetscInt f_start, f_end; // faces - PetscInt v_start, v_end; // vertices - +PetscErrorCode CreateDM(MPI_Comm comm, Ceed ceed, DM *dm) { PetscFunctionBeginUser; + CeedMemType mem_type_backend; + CeedGetPreferredMemType(ceed, &mem_type_backend); + + VecType vec_type = NULL; + MatType mat_type = NULL; + switch (mem_type_backend) { + case CEED_MEM_HOST: + vec_type = VECSTANDARD; + break; + case CEED_MEM_DEVICE: { + const char *resolved; + CeedGetResource(ceed, &resolved); + if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA; + else if (strstr(resolved, "/gpu/hip/occa")) vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678 + else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP; + else vec_type = VECSTANDARD; + } + } + if (strstr(vec_type, VECCUDA)) mat_type = MATAIJCUSPARSE; + else if (strstr(vec_type, VECKOKKOS)) mat_type = MATAIJKOKKOS; + else mat_type = MATAIJ; + // Create DMPLEX - ierr = DMCreate(comm, dm); - CHKERRQ(ierr); - ierr = DMSetType(*dm, DMPLEX); - CHKERRQ(ierr); + PetscCall(DMCreate(comm, dm)); + PetscCall(DMSetType(*dm, DMPLEX)); + PetscCall(DMSetMatType(*dm, mat_type)); + PetscCall(DMSetVecType(*dm, vec_type)); // Set Tensor elements - ierr = PetscOptionsSetValue(NULL, "-dm_plex_simplex", "0"); - CHKERRQ(ierr); + PetscCall(PetscOptionsSetValue(NULL, "-dm_plex_simplex", "0")); // Set CL options - ierr = DMSetFromOptions(*dm); - CHKERRQ(ierr); - ierr = DMViewFromOptions(*dm, NULL, "-dm_view"); - CHKERRQ(ierr); - - // Get plex limits - ierr = DMPlexGetChart(*dm, &p_start, &p_end); - CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(*dm, 0, &c_start, &c_end); - CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(*dm, 1, &f_start, &f_end); - CHKERRQ(ierr); - ierr = DMPlexGetDepthStratum(*dm, 0, &v_start, &v_end); - CHKERRQ(ierr); - // Create section - ierr = PetscSectionCreate(comm, &sec); - CHKERRQ(ierr); - ierr = PetscSectionSetNumFields(sec, 1); - CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 0, "Velocity"); - CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 0, 1); - CHKERRQ(ierr); - ierr = PetscSectionSetChart(sec, p_start, p_end); - CHKERRQ(ierr); - // Setup dofs per face - for (PetscInt f = f_start; f < f_end; f++) { - ierr = DMPlexGetConeSize(*dm, f, &dofs_per_face); - CHKERRQ(ierr); - ierr = PetscSectionSetFieldDof(sec, f, 0, dofs_per_face); - CHKERRQ(ierr); - ierr = PetscSectionSetDof(sec, f, dofs_per_face); - CHKERRQ(ierr); - } - ierr = PetscSectionSetUp(sec); - CHKERRQ(ierr); - ierr = DMSetSection(*dm, sec); - CHKERRQ(ierr); - ierr = PetscSectionDestroy(&sec); - CHKERRQ(ierr); + PetscCall(DMSetFromOptions(*dm)); + PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); PetscFunctionReturn(0); }; \ No newline at end of file diff --git a/examples/Hdiv-mass/src/setup-fe.c b/examples/Hdiv-mass/src/setup-fe.c new file mode 100644 index 0000000000..09bf4d78cd --- /dev/null +++ b/examples/Hdiv-mass/src/setup-fe.c @@ -0,0 +1,149 @@ +#include "../include/setup-fe.h" + +#include "petscerror.h" +// ----------------------------------------------------------------------------- +// Convert PETSc MemType to libCEED MemType +// ----------------------------------------------------------------------------- +CeedMemType MemTypeP2C(PetscMemType mem_type) { return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; } + +// --------------------------------------------------------------------------- +// Setup FE for H(div) space +// --------------------------------------------------------------------------- +PetscErrorCode SetupFEHdiv(AppCtx app_ctx, ProblemData problem_data, DM dm) { + PetscSection sec; + PetscInt dofs_per_face; + PetscInt p_start, p_end; + PetscInt c_start, c_end; // cells + PetscInt f_start, f_end; // faces + PetscInt v_start, v_end; // vertices + + PetscFunctionBeginUser; + + // Get plex limits + PetscCall(DMPlexGetChart(dm, &p_start, &p_end)); + PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); + PetscCall(DMPlexGetHeightStratum(dm, 1, &f_start, &f_end)); + PetscCall(DMPlexGetDepthStratum(dm, 0, &v_start, &v_end)); + // Create section + PetscCall(PetscSectionCreate(app_ctx->comm, &sec)); + PetscCall(PetscSectionSetNumFields(sec, 1)); + PetscCall(PetscSectionSetFieldName(sec, 0, "Velocity")); + PetscCall(PetscSectionSetFieldComponents(sec, 0, 1)); + PetscCall(PetscSectionSetChart(sec, p_start, p_end)); + // Setup dofs per face + for (PetscInt f = f_start; f < f_end; f++) { + PetscCall(DMPlexGetConeSize(dm, f, &dofs_per_face)); + PetscCall(PetscSectionSetFieldDof(sec, f, 0, dofs_per_face)); + PetscCall(PetscSectionSetDof(sec, f, dofs_per_face)); + } + PetscCall(PetscSectionSetUp(sec)); + PetscCall(DMSetSection(dm, sec)); + PetscCall(PetscSectionDestroy(&sec)); + + PetscFunctionReturn(0); +}; + +// ----------------------------------------------------------------------------- +// Utility function - convert from DMPolytopeType to CeedElemTopology +// ----------------------------------------------------------------------------- +CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type) { + switch (cell_type) { + case DM_POLYTOPE_TRIANGLE: + return CEED_TOPOLOGY_TRIANGLE; + case DM_POLYTOPE_QUADRILATERAL: + return CEED_TOPOLOGY_QUAD; + case DM_POLYTOPE_TETRAHEDRON: + return CEED_TOPOLOGY_TET; + case DM_POLYTOPE_HEXAHEDRON: + return CEED_TOPOLOGY_HEX; + default: + return 0; + } +}; + +// ----------------------------------------------------------------------------- +// Utility function - essential BC dofs are encoded in closure indices as -(i+1) +// ----------------------------------------------------------------------------- +PetscInt Involute(PetscInt i) { return i >= 0 ? i : -(i + 1); }; + +// ----------------------------------------------------------------------------- +// Get CEED restriction data from DMPlex +// ----------------------------------------------------------------------------- +PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) { + PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets; + + PetscFunctionBeginUser; + + PetscCall(DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets)); + CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, elem_restr_offsets, elem_restr); + PetscCall(PetscFree(elem_restr_offsets)); + + PetscFunctionReturn(0); +}; + +// ----------------------------------------------------------------------------- +// Get Oriented CEED restriction data from DMPlex +// ----------------------------------------------------------------------------- +PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, CeedInt P, CeedElemRestriction *elem_restr_u, CeedElemRestriction *elem_restr_p) { + PetscSection section; + PetscInt p, num_elem, num_dof, *restr_indices_u, *restr_indices_p, elem_offset, num_fields, dim, c_start, c_end; + Vec U_loc; + const PetscInt *ornt; // this is for orientation of dof + + PetscFunctionBeginUser; + + PetscCall(DMGetDimension(dm, &dim)); + PetscCall(DMGetLocalSection(dm, §ion)); + PetscCall(PetscSectionGetNumFields(section, &num_fields)); + PetscInt num_comp[num_fields], field_offsets[num_fields + 1]; + field_offsets[0] = 0; + for (PetscInt f = 0; f < num_fields; f++) { + PetscCall(PetscSectionGetFieldComponents(section, f, &num_comp[f])); + field_offsets[f + 1] = field_offsets[f] + num_comp[f]; + } + PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); + num_elem = c_end - c_start; + PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &restr_indices_u)); + PetscCall(PetscMalloc1(num_elem, &restr_indices_p)); + bool *orient_indices; // to flip the dof + PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &orient_indices)); + + for (p = 0, elem_offset = 0; p < num_elem; p++) { + restr_indices_p[p] = p; // each cell has on P0 dof + PetscInt num_indices, *indices, faces_per_elem, dofs_per_face; + PetscCall(DMPlexGetClosureIndices(dm, section, section, p, PETSC_TRUE, &num_indices, &indices, NULL, NULL)); + PetscCall(DMPlexGetConeOrientation(dm, p, &ornt)); + // Get number of faces per element + PetscCall(DMPlexGetConeSize(dm, p, &faces_per_elem)); + dofs_per_face = faces_per_elem - 2; + for (PetscInt f = 0; f < faces_per_elem; f++) { + for (PetscInt i = 0; i < dofs_per_face; i++) { + PetscInt ii = dofs_per_face * f + i; + // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. + PetscInt loc = Involute(indices[ii * num_comp[0]]); + restr_indices_u[elem_offset] = loc; + // Set orientation + orient_indices[elem_offset] = ornt[f] < 0; + elem_offset++; + } + } + PetscCall(DMPlexRestoreClosureIndices(dm, section, section, p, PETSC_TRUE, &num_indices, &indices, NULL, NULL)); + } + // if (elem_offset != num_elem*dim*PetscPowInt(P, dim)) + // SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, + // "ElemRestriction of size (%" PetscInt_FMT ",%" PetscInt_FMT ") + // initialized %" PetscInt_FMT "nodes", num_elem, + // dim*PetscPowInt(P, dim),elem_offset); + PetscCall(DMGetLocalVector(dm, &U_loc)); + PetscCall(VecGetLocalSize(U_loc, &num_dof)); + PetscCall(DMRestoreLocalVector(dm, &U_loc)); + // dof per element in Hdiv is dim*P^dim, for linear element P=2 + CeedElemRestrictionCreateOriented(ceed, num_elem, dim * PetscPowInt(P, dim), field_offsets[num_fields], 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, + restr_indices_u, orient_indices, elem_restr_u); + CeedElemRestrictionCreate(ceed, num_elem, 1, 1, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices_p, elem_restr_p); + PetscCall(PetscFree(restr_indices_u)); + PetscCall(PetscFree(orient_indices)); + PetscCall(PetscFree(restr_indices_p)); + + PetscFunctionReturn(0); +}; \ No newline at end of file diff --git a/examples/Hdiv-mass/src/setup-libceed.c b/examples/Hdiv-mass/src/setup-libceed.c index fd5c29b7b6..82ac374751 100644 --- a/examples/Hdiv-mass/src/setup-libceed.c +++ b/examples/Hdiv-mass/src/setup-libceed.c @@ -1,19 +1,15 @@ #include "../include/setup-libceed.h" +#include + #include "../basis/Hdiv-hex.h" #include "../basis/Hdiv-quad.h" -#include "../include/petsc-macros.h" +#include "../basis/L2-P0.h" -// ----------------------------------------------------------------------------- -// Convert PETSc MemType to libCEED MemType -// ----------------------------------------------------------------------------- -CeedMemType MemTypeP2C(PetscMemType mem_type) { return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; } // ----------------------------------------------------------------------------- // Destroy libCEED objects // ----------------------------------------------------------------------------- PetscErrorCode CeedDataDestroy(CeedData ceed_data) { - PetscErrorCode ierr; - PetscFunctionBegin; // Vectors @@ -23,133 +19,30 @@ PetscErrorCode CeedDataDestroy(CeedData ceed_data) { CeedElemRestrictionDestroy(&ceed_data->elem_restr_x); CeedElemRestrictionDestroy(&ceed_data->elem_restr_u); CeedElemRestrictionDestroy(&ceed_data->elem_restr_u_i); + CeedElemRestrictionDestroy(&ceed_data->elem_restr_p); // Bases CeedBasisDestroy(&ceed_data->basis_x); CeedBasisDestroy(&ceed_data->basis_u); + CeedBasisDestroy(&ceed_data->basis_p); // QFunctions CeedQFunctionDestroy(&ceed_data->qf_residual); CeedQFunctionDestroy(&ceed_data->qf_error); // Operators CeedOperatorDestroy(&ceed_data->op_residual); CeedOperatorDestroy(&ceed_data->op_error); - ierr = PetscFree(ceed_data); - CHKERRQ(ierr); - - PetscFunctionReturn(0); -}; - -// ----------------------------------------------------------------------------- -// Utility function - essential BC dofs are encoded in closure indices as -(i+1) -// ----------------------------------------------------------------------------- -PetscInt Involute(PetscInt i) { return i >= 0 ? i : -(i + 1); }; - -// ----------------------------------------------------------------------------- -// Get CEED restriction data from DMPlex -// ----------------------------------------------------------------------------- -PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) { - PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets; - PetscErrorCode ierr; - - PetscFunctionBeginUser; - - ierr = DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets); - CHKERRQ(ierr); + PetscCall(PetscFree(ceed_data)); - CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, elem_restr_offsets, elem_restr); - ierr = PetscFree(elem_restr_offsets); - CHKERRQ(ierr); - - PetscFunctionReturn(0); -}; - -// ----------------------------------------------------------------------------- -// Get Oriented CEED restriction data from DMPlex -// ----------------------------------------------------------------------------- -PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, CeedInt P, CeedElemRestriction *elem_restr_oriented) { - PetscSection section; - PetscInt p, num_elem, num_dof, *restr_indices, elem_offset, num_fields, dim, c_start, c_end; - Vec U_loc; - PetscErrorCode ierr; - const PetscInt *ornt; // this is for orientation of dof - PetscFunctionBeginUser; - ierr = DMGetDimension(dm, &dim); - CHKERRQ(ierr); - ierr = DMGetLocalSection(dm, §ion); - CHKERRQ(ierr); - ierr = PetscSectionGetNumFields(section, &num_fields); - CHKERRQ(ierr); - PetscInt num_comp[num_fields], field_offsets[num_fields + 1]; - field_offsets[0] = 0; - for (PetscInt f = 0; f < num_fields; f++) { - ierr = PetscSectionGetFieldComponents(section, f, &num_comp[f]); - CHKERRQ(ierr); - field_offsets[f + 1] = field_offsets[f] + num_comp[f]; - } - ierr = DMPlexGetHeightStratum(dm, 0, &c_start, &c_end); - CHKERRQ(ierr); - num_elem = c_end - c_start; - ierr = PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &restr_indices); - CHKERRQ(ierr); - bool *orient_indices; // to flip the dof - ierr = PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &orient_indices); - CHKERRQ(ierr); - for (p = 0, elem_offset = 0; p < num_elem; p++) { - PetscInt num_indices, *indices, faces_per_elem, dofs_per_face; - ierr = DMPlexGetClosureIndices(dm, section, section, p, PETSC_TRUE, &num_indices, &indices, NULL, NULL); - CHKERRQ(ierr); - - ierr = DMPlexGetConeOrientation(dm, p, &ornt); - CHKERRQ(ierr); - // Get number of faces per element - ierr = DMPlexGetConeSize(dm, p, &faces_per_elem); - CHKERRQ(ierr); - dofs_per_face = faces_per_elem - 2; - for (PetscInt f = 0; f < faces_per_elem; f++) { - for (PetscInt i = 0; i < dofs_per_face; i++) { - PetscInt ii = dofs_per_face * f + i; - // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. - PetscInt loc = Involute(indices[ii * num_comp[0]]); - restr_indices[elem_offset] = loc; - // Set orientation - orient_indices[elem_offset] = ornt[f] < 0; - elem_offset++; - } - } - ierr = DMPlexRestoreClosureIndices(dm, section, section, p, PETSC_TRUE, &num_indices, &indices, NULL, NULL); - CHKERRQ(ierr); - } - // if (elem_offset != num_elem*dim*PetscPowInt(P, dim)) - // SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, - // "ElemRestriction of size (%" PetscInt_FMT ",%" PetscInt_FMT ") - // initialized %" PetscInt_FMT "nodes", num_elem, - // dim*PetscPowInt(P, dim),elem_offset); - ierr = DMGetLocalVector(dm, &U_loc); - CHKERRQ(ierr); - ierr = VecGetLocalSize(U_loc, &num_dof); - CHKERRQ(ierr); - ierr = DMRestoreLocalVector(dm, &U_loc); - CHKERRQ(ierr); - // dof per element in Hdiv is dim*P^dim, for linear element P=2 - CeedElemRestrictionCreateOriented(ceed, num_elem, dim * PetscPowInt(P, dim), field_offsets[num_fields], 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, - restr_indices, orient_indices, elem_restr_oriented); - ierr = PetscFree(restr_indices); - CHKERRQ(ierr); - ierr = PetscFree(orient_indices); - CHKERRQ(ierr); PetscFunctionReturn(0); }; // ----------------------------------------------------------------------------- // Set up libCEED on the fine grid for a given degree // ----------------------------------------------------------------------------- -PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData *problem_data, PetscInt U_g_size, PetscInt U_loc_size, CeedData ceed_data, - CeedVector rhs_ceed, CeedVector *target) { - int ierr; +PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData problem_data, CeedData ceed_data, CeedVector rhs_ceed) { CeedInt P = app_ctx->degree + 1; // Number of quadratures in 1D, q_extra is set in cl-options.c - CeedInt Q = P + 1 + app_ctx->q_extra; - CeedInt dim, num_comp_x, num_comp_u; - // CeedInt elem_node = problem_data->elem_node; + CeedInt Q = P + 1 + app_ctx->q_extra; + CeedInt dim, num_comp_x, num_comp_u; DM dm_coord; Vec coords; PetscInt c_start, c_end, num_elem; @@ -169,58 +62,52 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData *probl CeedInt num_qpts = PetscPowInt(Q, dim); CeedInt P_u = dim * PetscPowInt(P, dim); // dof per element CeedScalar q_ref[dim * num_qpts], q_weights[num_qpts]; - CeedScalar div[P_u * num_qpts], interp[dim * P_u * num_qpts]; + CeedScalar div[P_u * num_qpts], interp[dim * P_u * num_qpts], interp_p[num_qpts], *grad = NULL; if (dim == 2) { HdivBasisQuad(Q, q_ref, q_weights, interp, div, problem_data->quadrature_mode); CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, num_comp_u, P_u, num_qpts, interp, div, q_ref, q_weights, &ceed_data->basis_u); + L2BasisP0(dim, Q, q_ref, q_weights, interp_p, problem_data->quadrature_mode); + CeedBasisCreateH1(ceed, CEED_TOPOLOGY_QUAD, 1, 1, num_qpts, interp_p, grad, q_ref, q_weights, &ceed_data->basis_p); } else { HdivBasisHex(Q, q_ref, q_weights, interp, div, problem_data->quadrature_mode); CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_HEX, num_comp_u, P_u, num_qpts, interp, div, q_ref, q_weights, &ceed_data->basis_u); + L2BasisP0(dim, Q, q_ref, q_weights, interp_p, problem_data->quadrature_mode); + CeedBasisCreateH1(ceed, CEED_TOPOLOGY_HEX, 1, 1, num_qpts, interp_p, grad, q_ref, q_weights, &ceed_data->basis_p); } - CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, problem_data->quadrature_mode, &ceed_data->basis_x); // --------------------------------------------------------------------------- // libCEED restrictions // --------------------------------------------------------------------------- - ierr = DMGetCoordinateDM(dm, &dm_coord); - CHKERRQ(ierr); - ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); - CHKERRQ(ierr); + PetscCall(DMGetCoordinateDM(dm, &dm_coord)); + PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL)); + CeedInt height = 0; // 0 means no boundary conditions DMLabel domain_label = 0; PetscInt value = 0; // -- Coordinate restriction - ierr = CreateRestrictionFromPlex(ceed, dm_coord, height, domain_label, value, &ceed_data->elem_restr_x); - CHKERRQ(ierr); - // -- Solution and projected true solution restriction - ierr = CreateRestrictionFromPlexOriented(ceed, dm, P, &ceed_data->elem_restr_u); - CHKERRQ(ierr); - // -- Geometric ceed_data restriction - ierr = DMPlexGetHeightStratum(dm, 0, &c_start, &c_end); - CHKERRQ(ierr); + PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, height, domain_label, value, &ceed_data->elem_restr_x)); + // -- Solution restriction, Error restriction + PetscCall(CreateRestrictionFromPlexOriented(ceed, dm, P, &ceed_data->elem_restr_u, &ceed_data->elem_restr_p)); + PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); num_elem = c_end - c_start; - + // -- Target restriction for MMS CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, dim, dim * num_elem * num_qpts, CEED_STRIDES_BACKEND, &ceed_data->elem_restr_u_i); - // --------------------------------------------------------------------------- // Element coordinates // --------------------------------------------------------------------------- - ierr = DMGetCoordinatesLocal(dm, &coords); - CHKERRQ(ierr); - ierr = VecGetArrayRead(coords, &coordArray); - CHKERRQ(ierr); - + PetscCall(DMGetCoordinatesLocal(dm, &coords)); + PetscCall(VecGetArrayRead(coords, &coordArray)); CeedElemRestrictionCreateVector(ceed_data->elem_restr_x, &x_coord, NULL); CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *)coordArray); - ierr = VecRestoreArrayRead(coords, &coordArray); - CHKERRQ(ierr); + PetscCall(VecRestoreArrayRead(coords, &coordArray)); // --------------------------------------------------------------------------- // Setup RHS and true solution // --------------------------------------------------------------------------- - CeedVectorCreate(ceed, num_elem * num_qpts * dim, target); + CeedVector target; + CeedVectorCreate(ceed, num_elem * num_qpts * dim, &target); // Create the q-function that sets up the RHS and true solution CeedQFunctionCreateInterior(ceed, 1, problem_data->setup_rhs, problem_data->setup_rhs_loc, &qf_setup_rhs); CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); @@ -233,7 +120,7 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData *probl CeedOperatorSetField(op_setup_rhs, "x", ceed_data->elem_restr_x, ceed_data->basis_x, CEED_VECTOR_ACTIVE); CeedOperatorSetField(op_setup_rhs, "weight", CEED_ELEMRESTRICTION_NONE, ceed_data->basis_x, CEED_VECTOR_NONE); CeedOperatorSetField(op_setup_rhs, "dx", ceed_data->elem_restr_x, ceed_data->basis_x, CEED_VECTOR_ACTIVE); - CeedOperatorSetField(op_setup_rhs, "true_soln", ceed_data->elem_restr_u_i, CEED_BASIS_COLLOCATED, *target); + CeedOperatorSetField(op_setup_rhs, "true_soln", ceed_data->elem_restr_u_i, CEED_BASIS_COLLOCATED, target); CeedOperatorSetField(op_setup_rhs, "rhs", ceed_data->elem_restr_u, ceed_data->basis_u, CEED_VECTOR_ACTIVE); // Setup RHS and true solution @@ -242,9 +129,9 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData *probl // --------------------------------------------------------------------------- // Persistent libCEED vectors // --------------------------------------------------------------------------- - // -- Operator action variables - CeedVectorCreate(ceed, U_loc_size, &ceed_data->x_ceed); - CeedVectorCreate(ceed, U_loc_size, &ceed_data->y_ceed); + // -- Operator action variables: we use them in setup-solvers.c + CeedElemRestrictionCreateVector(ceed_data->elem_restr_u, &ceed_data->x_ceed, NULL); + CeedElemRestrictionCreateVector(ceed_data->elem_restr_u, &ceed_data->y_ceed, NULL); // Local residual evaluator // --------------------------------------------------------------------------- @@ -277,14 +164,16 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData *probl CeedQFunctionAddInput(qf_error, "u", dim, CEED_EVAL_INTERP); CeedQFunctionAddInput(qf_error, "true_soln", dim, CEED_EVAL_NONE); CeedQFunctionAddInput(qf_error, "weight", 1, CEED_EVAL_WEIGHT); - CeedQFunctionAddOutput(qf_error, "error", dim, CEED_EVAL_NONE); + // CeedQFunctionAddOutput(qf_error, "error", 1, CEED_EVAL_NONE); + CeedQFunctionAddOutput(qf_error, "error", 1, CEED_EVAL_INTERP); // Create the operator that builds the error CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_error); CeedOperatorSetField(op_error, "dx", ceed_data->elem_restr_x, ceed_data->basis_x, x_coord); CeedOperatorSetField(op_error, "u", ceed_data->elem_restr_u, ceed_data->basis_u, CEED_VECTOR_ACTIVE); - CeedOperatorSetField(op_error, "true_soln", ceed_data->elem_restr_u_i, CEED_BASIS_COLLOCATED, *target); + CeedOperatorSetField(op_error, "true_soln", ceed_data->elem_restr_u_i, CEED_BASIS_COLLOCATED, target); CeedOperatorSetField(op_error, "weight", CEED_ELEMRESTRICTION_NONE, ceed_data->basis_x, CEED_VECTOR_NONE); - CeedOperatorSetField(op_error, "error", ceed_data->elem_restr_u_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); + // CeedOperatorSetField(op_error, "error", ceed_data->elem_restr_e_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_error, "error", ceed_data->elem_restr_p, ceed_data->basis_p, CEED_VECTOR_ACTIVE); // -- Save libCEED data to apply operator in matops.c ceed_data->qf_error = qf_error; ceed_data->op_error = op_error; @@ -292,6 +181,7 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData *probl CeedQFunctionDestroy(&qf_setup_rhs); CeedOperatorDestroy(&op_setup_rhs); CeedVectorDestroy(&x_coord); + CeedVectorDestroy(&target); PetscFunctionReturn(0); }; diff --git a/examples/Hdiv-mass/src/setup-matops.c b/examples/Hdiv-mass/src/setup-matops.c new file mode 100644 index 0000000000..27fb84f3ba --- /dev/null +++ b/examples/Hdiv-mass/src/setup-matops.c @@ -0,0 +1,51 @@ +#include "../include/setup-matops.h" + +#include + +// ----------------------------------------------------------------------------- +// Apply the local action of a libCEED operator and store result in PETSc vector +// i.e. compute A X = Y +// ----------------------------------------------------------------------------- +PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, OperatorApplyContext op_apply_ctx) { + PetscFunctionBeginUser; + + // Zero target vector + PetscCall(VecZeroEntries(Y)); + + // Sum into target vector + PetscCall(ApplyAddLocalCeedOp(X, Y, op_apply_ctx)); + + PetscFunctionReturn(0); +} + +PetscErrorCode ApplyAddLocalCeedOp(Vec X, Vec Y, OperatorApplyContext op_apply_ctx) { + PetscScalar *x, *y; + PetscMemType x_mem_type, y_mem_type; + + PetscFunctionBeginUser; + + // Global-to-local + PetscCall(DMGlobalToLocal(op_apply_ctx->dm, X, INSERT_VALUES, op_apply_ctx->X_loc)); + + // Setup libCEED vectors + PetscCall(VecGetArrayReadAndMemType(op_apply_ctx->X_loc, (const PetscScalar **)&x, &x_mem_type)); + PetscCall(VecGetArrayAndMemType(op_apply_ctx->Y_loc, &y, &y_mem_type)); + CeedVectorSetArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); + CeedVectorSetArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); + + // Apply libCEED operator + CeedOperatorApply(op_apply_ctx->op_apply, op_apply_ctx->x_ceed, op_apply_ctx->y_ceed, CEED_REQUEST_IMMEDIATE); + + // Restore PETSc vectors + CeedVectorTakeArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), NULL); + CeedVectorTakeArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), NULL); + PetscCall(VecRestoreArrayReadAndMemType(op_apply_ctx->X_loc, (const PetscScalar **)&x)); + PetscCall(VecRestoreArrayAndMemType(op_apply_ctx->Y_loc, &y)); + + // Local-to-global + PetscCall(DMLocalToGlobal(op_apply_ctx->dm, op_apply_ctx->Y_loc, ADD_VALUES, Y)); + + PetscFunctionReturn(0); +}; + +// ----------------------------------------------------------------------------- \ No newline at end of file diff --git a/examples/Hdiv-mass/src/setup-solvers.c b/examples/Hdiv-mass/src/setup-solvers.c new file mode 100644 index 0000000000..fac3e91970 --- /dev/null +++ b/examples/Hdiv-mass/src/setup-solvers.c @@ -0,0 +1,123 @@ +#include "../include/setup-solvers.h" + +#include "../include/setup-libceed.h" +#include "../include/setup-matops.h" +#include "petscvec.h" + +// ----------------------------------------------------------------------------- +// Setup operator context data +// ----------------------------------------------------------------------------- +PetscErrorCode SetupResidualOperatorCtx(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_residual) { + PetscFunctionBeginUser; + + ctx_residual->comm = comm; + ctx_residual->dm = dm; + PetscCall(DMCreateLocalVector(dm, &ctx_residual->X_loc)); + PetscCall(VecDuplicate(ctx_residual->X_loc, &ctx_residual->Y_loc)); + ctx_residual->x_ceed = ceed_data->x_ceed; + ctx_residual->y_ceed = ceed_data->y_ceed; + ctx_residual->ceed = ceed; + ctx_residual->op_apply = ceed_data->op_residual; + + PetscFunctionReturn(0); +} + +PetscErrorCode SetupErrorOperatorCtx(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_error_u) { + PetscFunctionBeginUser; + + ctx_error_u->comm = comm; + ctx_error_u->dm = dm; + PetscCall(DMCreateLocalVector(dm, &ctx_error_u->X_loc)); + PetscCall(VecDuplicate(ctx_error_u->X_loc, &ctx_error_u->Y_loc)); + ctx_error_u->x_ceed = ceed_data->x_ceed; + ctx_error_u->y_ceed = ceed_data->y_ceed; + ctx_error_u->ceed = ceed; + ctx_error_u->op_apply = ceed_data->op_error; + + PetscFunctionReturn(0); +} + +// ----------------------------------------------------------------------------- +// This function wraps the libCEED operator for a MatShell +// ----------------------------------------------------------------------------- +PetscErrorCode ApplyMatOp(Mat A, Vec X, Vec Y) { + OperatorApplyContext op_apply_ctx; + + PetscFunctionBeginUser; + + PetscCall(MatShellGetContext(A, &op_apply_ctx)); + + // libCEED for local action of residual evaluator + PetscCall(ApplyLocalCeedOp(X, Y, op_apply_ctx)); + + PetscFunctionReturn(0); +}; + +// --------------------------------------------------------------------------- +// Setup Solver +// --------------------------------------------------------------------------- +PetscErrorCode PDESolver(CeedData ceed_data, AppCtx app_ctx, KSP ksp, Vec rhs, Vec *X) { + PetscInt X_l_size, X_g_size; + + PetscFunctionBeginUser; + + // Create global unknown solution U + PetscCall(VecGetSize(*X, &X_g_size)); + // Local size for matShell + PetscCall(VecGetLocalSize(*X, &X_l_size)); + + // --------------------------------------------------------------------------- + // Setup SNES + // --------------------------------------------------------------------------- + // Operator + Mat mat_op; + VecType vec_type; + PetscCall(DMGetVecType(app_ctx->ctx_residual->dm, &vec_type)); + // -- Form Action of Jacobian on delta_u + PetscCall(MatCreateShell(app_ctx->comm, X_l_size, X_l_size, X_g_size, X_g_size, app_ctx->ctx_residual, &mat_op)); + PetscCall(MatShellSetOperation(mat_op, MATOP_MULT, (void (*)(void))ApplyMatOp)); + PetscCall(MatShellSetVecType(mat_op, vec_type)); + + PC pc; + PetscCall(KSPGetPC(ksp, &pc)); + PetscCall(PCSetType(pc, PCNONE)); + PetscCall(KSPSetType(ksp, KSPCG)); + // PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); + PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); + + PetscCall(KSPSetOperators(ksp, mat_op, mat_op)); + PetscCall(KSPSetFromOptions(ksp)); + PetscCall(VecZeroEntries(*X)); + PetscCall(KSPSolve(ksp, rhs, *X)); + + PetscCall(MatDestroy(&mat_op)); + + PetscFunctionReturn(0); +}; + +// ----------------------------------------------------------------------------- +// This function calculates the error in the final solution +// ----------------------------------------------------------------------------- +PetscErrorCode ComputeL2Error(Vec X, PetscScalar *l2_error, OperatorApplyContext op_error_ctx) { + Vec E; + PetscFunctionBeginUser; + + PetscCall(VecDuplicate(X, &E)); + PetscCall(ApplyLocalCeedOp(X, E, op_error_ctx)); + PetscScalar error_sq = 1.0; + PetscCall(VecSum(E, &error_sq)); + *l2_error = sqrt(error_sq); + PetscCall(VecDestroy(&E)); + PetscFunctionReturn(0); +}; + +PetscErrorCode CtxVecDestroy(ProblemData problem_data, AppCtx app_ctx) { + PetscFunctionBegin; + + PetscCall(VecDestroy(&app_ctx->ctx_residual->Y_loc)); + PetscCall(VecDestroy(&app_ctx->ctx_residual->X_loc)); + PetscCall(VecDestroy(&app_ctx->ctx_error_u->Y_loc)); + PetscCall(VecDestroy(&app_ctx->ctx_error_u->X_loc)); + + PetscFunctionReturn(0); +} \ No newline at end of file