some clean-up in Hdiv-mass folder

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 \

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;
+        }
+      }
+    }
+  }
+}

examples/Hdiv-mass/conv_plot.py → 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()

examples/Hdiv-mass/conv_test.sh

@@ -24,27 +24,27 @@ 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
         run_flags[dm_plex_box_faces]=2,2,2
     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
 

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

examples/Hdiv-mass/include/post-processing.h

@@ -0,0 +1,12 @@
+#ifndef post_processing_h
+#define post_processing_h
+
+#include <ceed.h>
+#include <petsc.h>
+
+#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

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

examples/Hdiv-mass/include/setup-dm.h

@@ -6,11 +6,11 @@
 #include <petscdmplex.h>
 #include <petscsys.h>
 
-#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

examples/Hdiv-mass/include/setup-fe.h

@@ -0,0 +1,20 @@
+#ifndef setupfe_h
+#define setupfe_h
+
+#include <ceed.h>
+#include <petsc.h>
+#include <petscdmplex.h>
+#include <petscsys.h>
+
+#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

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

examples/Hdiv-mass/include/setup-matops.h

@@ -0,0 +1,13 @@
+#ifndef setup_matops_h
+#define setup_matops_h
+
+#include <ceed.h>
+#include <petsc.h>
+
+#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

examples/Hdiv-mass/include/setup-solvers.h

@@ -0,0 +1,16 @@
+#ifndef setup_solvers_h
+#define setup_solvers_h
+
+#include <ceed.h>
+#include <petsc.h>
+
+#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