Added random vertices distribution

examples/Hdiv-mixed/conv_plot.py

@@ -46,6 +46,7 @@ def plot():
     E_p = data['error_p']
     #E_hdiv = data['error_hdiv']
     h = 1/data['mesh_res']
+    N = data['mesh_res']
     H1 =  amin(E_p)* (h/amin(h)) # H = C h^1
     H2 =  amin(E_u)* (h/amin(h))**2  # H = C h^2
 
@@ -62,6 +63,18 @@ def plot():
     #xlim(.06, .3)
     fig.tight_layout()
     plt.savefig('convrate_mixed.png', bbox_inches='tight')
+
+    conv_u = []
+    conv_p = []
+    conv_u.append(0)
+    conv_p.append(0)
+    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]))
+        conv_p.append(log10(E_p[i]/E_p[i-1])/log10(h[i]/h[i-1]))
+
+    result = {'Number of element':N, 'order of u':conv_u, 'order of p':conv_p}
+    df = pd.DataFrame(result)
+    print(df)
 
 
 

examples/Hdiv-mixed/conv_test.sh

@@ -42,13 +42,13 @@ declare -A run_flags
         run_flags[dm_plex_dim]=$dim
         run_flags[dm_plex_box_faces]=2,2,2
         run_flags[dm_plex_box_lower]=0,0,0
-        run_flags[dm_plex_box_upper]=1,0.1,1
+        run_flags[dm_plex_box_upper]=1,0.2,1
     fi
 
 declare -A test_flags
-    test_flags[res_start]=2
+    test_flags[res_start]=4
     test_flags[res_stride]=1
-    test_flags[res_end]=10
+    test_flags[res_end]=8
 
 file_name=conv_test_result.csv
 

examples/Hdiv-mixed/conv_test_result.csv

@@ -1,10 +1,6 @@
 run,mesh_res,error_u,error_p
-0,2,23.34221,0.01392
-1,3,14.88356,0.00793
-2,4,11.04429,0.00485
-3,5,8.90114,0.00337
-4,6,6.49004,0.00260
-5,7,4.87041,0.00207
-6,8,3.99474,0.00165
-7,9,3.24332,0.00136
-8,10,2.65430,0.00114
+0,4,9.10868,0.02425
+1,5,7.66382,0.02087
+2,6,5.80589,0.01738
+3,7,5.13953,0.01512
+4,8,4.74095,0.01357

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

@@ -12,5 +12,6 @@
 // ---------------------------------------------------------------------------
 PetscErrorCode CreateDM(MPI_Comm comm, VecType vec_type, DM *dm);
 PetscErrorCode PerturbVerticesSmooth(DM dm);
+PetscErrorCode PerturbVerticesRandom(DM dm);
 
 #endif // setupdm_h

examples/Hdiv-mixed/main.c

@@ -30,6 +30,7 @@
 // (boundary is not working)
 //   ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1
 //   ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1,2,3,4
+//   ./main -pc_type svd -problem darcy3d -dm_plex_dim 3 -dm_plex_box_faces 4,4,4 -view_solution -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,0.1,1
 const char help[] = "Solve H(div)-mixed problem using PETSc and libCEED\n";
 
 #include "main.h"
@@ -83,8 +84,8 @@ int main(int argc, char **argv) {
   // ---------------------------------------------------------------------------
   // -- Initialize backend
   Ceed ceed;
-  //CeedInit("/cpu/self/ref/serial", &ceed);
-  CeedInit(app_ctx->ceed_resource, &ceed);
+  CeedInit("/cpu/self/ref/serial", &ceed);
+  //CeedInit(app_ctx->ceed_resource, &ceed);
   CeedMemType mem_type_backend;
   CeedGetPreferredMemType(ceed, &mem_type_backend);
 
@@ -118,9 +119,11 @@ int main(int argc, char **argv) {
   PetscCall( CreateDM(comm, vec_type, &dm_H1) );
   // TODO: add mesh option
   // perturb to have smooth random mesh
-  PetscCall( PerturbVerticesSmooth(dm) );
-  PetscCall( PerturbVerticesSmooth(dm_H1) );
+  //PetscCall( PerturbVerticesSmooth(dm) );
+  //PetscCall( PerturbVerticesSmooth(dm_H1) );
 
+  PetscCall(PerturbVerticesRandom(dm) );
+  PetscCall(PerturbVerticesRandom(dm_H1) );
   // ---------------------------------------------------------------------------
   // Process command line options
   // ---------------------------------------------------------------------------

examples/Hdiv-mixed/problems/darcy2d.c

@@ -91,7 +91,6 @@ PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, DM dm,
   PetscCall( DMGetBoundingBox(dm, domain_min, domain_max) );
   for (PetscInt i=0; i<2; i++) domain_size[i] = domain_max[i] - domain_min[i];
 
-  printf("lx %f, ly %f \n",domain_size[0], domain_size[1]);
   darcy_ctx->kappa = kappa;
   darcy_ctx->rho_a0 = rho_a0;
   darcy_ctx->g = g;

examples/Hdiv-mixed/src/setup-dm.c

@@ -68,3 +68,59 @@ PetscErrorCode PerturbVerticesSmooth(DM dm) {
   PetscCall( VecRestoreArray(coordinates,&coords) );
   PetscFunctionReturn(0);
 }
+
+
+PetscErrorCode PerturbVerticesRandom(DM dm) {
+
+  PetscFunctionBegin;
+  Vec          coordinates;
+  PetscSection coordSection;
+  PetscScalar *coords;
+  PetscInt     v,vStart,vEnd,offset, dim;
+  PetscReal    x, y, z;
+
+  PetscCall( DMGetDimension(dm,&dim) );
+  PetscCall( DMGetCoordinateSection(dm, &coordSection) );
+  PetscCall( DMGetCoordinatesLocal(dm, &coordinates) );
+  PetscCall( DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd) );
+  PetscCall( VecGetArray(coordinates,&coords) );
+  PetscInt c_end, c_start, num_elem;
+  PetscCall( DMPlexGetHeightStratum(dm, 0, &c_start, &c_end) );
+  num_elem = c_end - c_start;
+
+  for(v=vStart; v<vEnd; v++) {
+    PetscCall( PetscSectionGetOffset(coordSection,v,&offset) );
+    if(dim==2) {
+      PetscScalar nx = sqrt(num_elem);
+      PetscReal domain_min[2], domain_max[2], domain_size[2];
+      PetscCall( DMGetBoundingBox(dm, domain_min, domain_max) );
+      for (PetscInt i=0; i<2; i++) domain_size[i] = domain_max[i] - domain_min[i];
+      PetscReal hx = domain_size[0]/nx, hy = domain_size[1]/nx;
+      x = coords[offset]; y = coords[offset+1];
+      // perturb randomly O(h*sqrt(2)/3)
+      PetscReal rx = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hx*0.471404);
+      PetscReal ry = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hy*0.471404);
+      PetscReal t = ((PetscReal)rand())/((PetscReal)RAND_MAX)*PETSC_PI;
+      coords[offset  ] = x + rx*PetscCosReal(t);
+      coords[offset+1] = y + ry*PetscSinReal(t);
+    } else {
+      PetscScalar nx = cbrt(num_elem);
+      PetscReal domain_min[3], domain_max[3], domain_size[3];
+      PetscCall( DMGetBoundingBox(dm, domain_min, domain_max) );
+      for (PetscInt i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i];
+      PetscReal hx = domain_size[0]/nx, hy = domain_size[1]/nx,
+                hz = domain_size[2]/nx;
+      x = coords[offset]; y = coords[offset+1], z = coords[offset+2];
+      // This is because 'boundary' is broken in 3D
+      PetscReal rx = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hx*0.471404);
+      PetscReal ry = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hy*0.471404);
+      PetscReal rz = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hz*0.471404);
+      PetscReal t = ((PetscReal)rand())/((PetscReal)RAND_MAX)*PETSC_PI;
+      coords[offset  ] = x + rx*PetscCosReal(t);
+      coords[offset+1] = y + ry*PetscCosReal(t);
+      coords[offset+2] = z + rz*PetscSinReal(t);
+    }
+  }
+  PetscCall( VecRestoreArray(coordinates,&coords) );
+  PetscFunctionReturn(0);
+}