diff --git a/examples/Hdiv-mixed/conv_plot.py b/examples/Hdiv-mixed/conv_plot.py index d6c7f6653b..d4bfcff84e 100644 --- a/examples/Hdiv-mixed/conv_plot.py +++ b/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) diff --git a/examples/Hdiv-mixed/conv_test.sh b/examples/Hdiv-mixed/conv_test.sh index 3acfecfc01..36c1d52723 100755 --- a/examples/Hdiv-mixed/conv_test.sh +++ b/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 diff --git a/examples/Hdiv-mixed/conv_test_result.csv b/examples/Hdiv-mixed/conv_test_result.csv index 3c3e88974f..50df6b824a 100644 --- a/examples/Hdiv-mixed/conv_test_result.csv +++ b/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 diff --git a/examples/Hdiv-mixed/convrate_mixed.png b/examples/Hdiv-mixed/convrate_mixed.png index 6c468b7216..3dab09a952 100644 Binary files a/examples/Hdiv-mixed/convrate_mixed.png and b/examples/Hdiv-mixed/convrate_mixed.png differ diff --git a/examples/Hdiv-mixed/include/setup-dm.h b/examples/Hdiv-mixed/include/setup-dm.h index 610d74efd0..40029bc7c4 100644 --- a/examples/Hdiv-mixed/include/setup-dm.h +++ b/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 diff --git a/examples/Hdiv-mixed/main.c b/examples/Hdiv-mixed/main.c index abd2739824..2a303d2ae0 100644 --- a/examples/Hdiv-mixed/main.c +++ b/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 // --------------------------------------------------------------------------- diff --git a/examples/Hdiv-mixed/problems/darcy2d.c b/examples/Hdiv-mixed/problems/darcy2d.c index 77dc8d2a74..16d61dbb62 100644 --- a/examples/Hdiv-mixed/problems/darcy2d.c +++ b/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; diff --git a/examples/Hdiv-mixed/src/setup-dm.c b/examples/Hdiv-mixed/src/setup-dm.c index c5304a17cb..eacb6c7789 100644 --- a/examples/Hdiv-mixed/src/setup-dm.c +++ b/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