Skip to content

Commit a5937fa

Browse files
rezgarshakerirezgarshakeri
committed
Added random vertices distribution
1 parent fe521b6 commit a5937fa

File tree

10 files changed

+108
-53
lines changed

10 files changed

+108
-53
lines changed

‎examples/Hdiv-mixed/conv_plot.py‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,7 @@ def plot():
4646
E_p = data['error_p']
4747
#E_hdiv = data['error_hdiv']
4848
h = 1/data['mesh_res']
49+
N = data['mesh_res']
4950
H1 = amin(E_p)* (h/amin(h)) # H = C h^1
5051
H2 = amin(E_u)* (h/amin(h))**2 # H = C h^2
5152

@@ -62,6 +63,18 @@ def plot():
6263
#xlim(.06, .3)
6364
fig.tight_layout()
6465
plt.savefig('convrate_mixed.png', bbox_inches='tight')
66+
67+
conv_u = []
68+
conv_p = []
69+
conv_u.append(0)
70+
conv_p.append(0)
71+
for i in range(1,len(E_u)):
72+
conv_u.append(log10(E_u[i]/E_u[i-1])/log10(h[i]/h[i-1]))
73+
conv_p.append(log10(E_p[i]/E_p[i-1])/log10(h[i]/h[i-1]))
74+
75+
result = {'Number of element':N, 'order of u':conv_u, 'order of p':conv_p}
76+
df = pd.DataFrame(result)
77+
print(df)
6578

6679

6780

‎examples/Hdiv-mixed/conv_test.sh‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -36,19 +36,19 @@ declare -A run_flags
3636
run_flags[dm_plex_dim]=$dim
3737
run_flags[dm_plex_box_faces]=2,2
3838
run_flags[dm_plex_box_lower]=0,0
39-
run_flags[dm_plex_box_upper]=1,0.1
39+
run_flags[dm_plex_box_upper]=1,1
4040
else
4141
run_flags[problem]=darcy3d
4242
run_flags[dm_plex_dim]=$dim
4343
run_flags[dm_plex_box_faces]=2,2,2
4444
run_flags[dm_plex_box_lower]=0,0,0
45-
run_flags[dm_plex_box_upper]=1,0.1,1
45+
run_flags[dm_plex_box_upper]=1,1,1
4646
fi
4747

4848
declare -A test_flags
49-
test_flags[res_start]=2
50-
test_flags[res_stride]=1
51-
test_flags[res_end]=10
49+
test_flags[res_start]=4
50+
test_flags[res_stride]=2
51+
test_flags[res_end]=12
5252

5353
file_name=conv_test_result.csv
5454

‎examples/Hdiv-mixed/conv_test_result.csv‎

Lines changed: 5 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,6 @@
11
run,mesh_res,error_u,error_p
2-
0,2,23.34221,0.01392
3-
1,3,14.88356,0.00793
4-
2,4,11.04429,0.00485
5-
3,5,8.90114,0.00337
6-
4,6,6.49004,0.00260
7-
5,7,4.87041,0.00207
8-
6,8,3.99474,0.00165
9-
7,9,3.24332,0.00136
10-
8,10,2.65430,0.00114
2+
0,4,26.30005,0.03133
3+
1,6,11.97420,0.01464
4+
2,8,6.79226,0.00838
5+
3,10,4.36393,0.00540
6+
4,12,3.03689,0.00377

‎examples/Hdiv-mixed/convrate_mixed.png‎

-2.74 KB
Loading

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,5 +12,6 @@
1212
// ---------------------------------------------------------------------------
1313
PetscErrorCode CreateDM(MPI_Comm comm, VecType vec_type, DM *dm);
1414
PetscErrorCode PerturbVerticesSmooth(DM dm);
15+
PetscErrorCode PerturbVerticesRandom(DM dm);
1516

1617
#endif // setupdm_h

‎examples/Hdiv-mixed/main.c‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -24,12 +24,14 @@
2424
// Build with: make
2525
// Run with:
2626
// ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4
27+
// ./main -pc_type none -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -ksp_type minres
2728
// ./main -pc_type svd -problem darcy3d -dm_plex_dim 3 -dm_plex_box_faces 4,4,4
2829
// ./main -pc_type svd -problem darcy3d -dm_plex_filename /path to the mesh file
2930
// ./main -pc_type svd -problem richard2d -dm_plex_dim 2 -dm_plex_box_faces 4,4
3031
// (boundary is not working)
3132
// ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1
3233
// ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1,2,3,4
34+
// ./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
3335
const char help[] = "Solve H(div)-mixed problem using PETSc and libCEED\n";
3436

3537
#include "main.h"
@@ -83,8 +85,8 @@ int main(int argc, char **argv) {
8385
// ---------------------------------------------------------------------------
8486
// -- Initialize backend
8587
Ceed ceed;
86-
//CeedInit("/cpu/self/ref/serial", &ceed);
87-
CeedInit(app_ctx->ceed_resource, &ceed);
88+
CeedInit("/cpu/self/ref/serial", &ceed);
89+
//CeedInit(app_ctx->ceed_resource, &ceed);
8890
CeedMemType mem_type_backend;
8991
CeedGetPreferredMemType(ceed, &mem_type_backend);
9092

@@ -118,9 +120,11 @@ int main(int argc, char **argv) {
118120
PetscCall( CreateDM(comm, vec_type, &dm_H1) );
119121
// TODO: add mesh option
120122
// perturb to have smooth random mesh
121-
PetscCall( PerturbVerticesSmooth(dm) );
122-
PetscCall( PerturbVerticesSmooth(dm_H1) );
123+
//PetscCall( PerturbVerticesSmooth(dm) );
124+
//PetscCall( PerturbVerticesSmooth(dm_H1) );
123125

126+
//PetscCall(PerturbVerticesRandom(dm) );
127+
//PetscCall(PerturbVerticesRandom(dm_H1) );
124128
// ---------------------------------------------------------------------------
125129
// Process command line options
126130
// ---------------------------------------------------------------------------

‎examples/Hdiv-mixed/problems/darcy2d.c‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,6 @@ PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, DM dm,
9191
PetscCall( DMGetBoundingBox(dm, domain_min, domain_max) );
9292
for (PetscInt i=0; i<2; i++) domain_size[i] = domain_max[i] - domain_min[i];
9393

94-
printf("lx %f, ly %f \n",domain_size[0], domain_size[1]);
9594
darcy_ctx->kappa = kappa;
9695
darcy_ctx->rho_a0 = rho_a0;
9796
darcy_ctx->g = g;

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

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -68,3 +68,59 @@ PetscErrorCode PerturbVerticesSmooth(DM dm) {
6868
PetscCall( VecRestoreArray(coordinates,&coords) );
6969
PetscFunctionReturn(0);
7070
}
71+
72+
73+
PetscErrorCode PerturbVerticesRandom(DM dm) {
74+
75+
PetscFunctionBegin;
76+
Vec coordinates;
77+
PetscSection coordSection;
78+
PetscScalar *coords;
79+
PetscInt v,vStart,vEnd,offset, dim;
80+
PetscReal x, y, z;
81+
82+
PetscCall( DMGetDimension(dm,&dim) );
83+
PetscCall( DMGetCoordinateSection(dm, &coordSection) );
84+
PetscCall( DMGetCoordinatesLocal(dm, &coordinates) );
85+
PetscCall( DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd) );
86+
PetscCall( VecGetArray(coordinates,&coords) );
87+
PetscInt c_end, c_start, num_elem;
88+
PetscCall( DMPlexGetHeightStratum(dm, 0, &c_start, &c_end) );
89+
num_elem = c_end - c_start;
90+
91+
for(v=vStart; v<vEnd; v++) {
92+
PetscCall( PetscSectionGetOffset(coordSection,v,&offset) );
93+
if(dim==2) {
94+
PetscScalar nx = sqrt(num_elem);
95+
PetscReal domain_min[2], domain_max[2], domain_size[2];
96+
PetscCall( DMGetBoundingBox(dm, domain_min, domain_max) );
97+
for (PetscInt i=0; i<2; i++) domain_size[i] = domain_max[i] - domain_min[i];
98+
PetscReal hx = domain_size[0]/nx, hy = domain_size[1]/nx;
99+
x = coords[offset]; y = coords[offset+1];
100+
// perturb randomly O(h*sqrt(2)/3)
101+
PetscReal rx = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hx*0.471404);
102+
PetscReal ry = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hy*0.471404);
103+
PetscReal t = ((PetscReal)rand())/((PetscReal)RAND_MAX)*PETSC_PI;
104+
coords[offset ] = x + rx*PetscCosReal(t);
105+
coords[offset+1] = y + ry*PetscSinReal(t);
106+
} else {
107+
PetscScalar nx = cbrt(num_elem);
108+
PetscReal domain_min[3], domain_max[3], domain_size[3];
109+
PetscCall( DMGetBoundingBox(dm, domain_min, domain_max) );
110+
for (PetscInt i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i];
111+
PetscReal hx = domain_size[0]/nx, hy = domain_size[1]/nx,
112+
hz = domain_size[2]/nx;
113+
x = coords[offset]; y = coords[offset+1], z = coords[offset+2];
114+
// This is because 'boundary' is broken in 3D
115+
PetscReal rx = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hx*0.471404);
116+
PetscReal ry = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hy*0.471404);
117+
PetscReal rz = ((PetscReal)rand())/((PetscReal)RAND_MAX)*(hz*0.471404);
118+
PetscReal t = ((PetscReal)rand())/((PetscReal)RAND_MAX)*PETSC_PI;
119+
coords[offset ] = x + rx*PetscCosReal(t);
120+
coords[offset+1] = y + ry*PetscCosReal(t);
121+
coords[offset+2] = z + rz*PetscSinReal(t);
122+
}
123+
}
124+
PetscCall( VecRestoreArray(coordinates,&coords) );
125+
PetscFunctionReturn(0);
126+
}

‎interface/ceed-basis.c‎

Lines changed: 18 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -1389,35 +1389,25 @@ int CeedBasisGetQWeights(CeedBasis basis, const CeedScalar **q_weight) {
13891389
@ref Advanced
13901390
**/
13911391
int CeedBasisGetInterp(CeedBasis basis, const CeedScalar **interp) {
1392-
switch (basis->basis_space) {
1393-
case 1: // H^1 discretization
1394-
if (!basis->interp && basis->tensor_basis) {
1395-
// Allocate
1396-
int ierr;
1397-
ierr = CeedMalloc(basis->Q*basis->P, &basis->interp); CeedChk(ierr);
1398-
1399-
// Initialize
1400-
for (CeedInt i=0; i<basis->Q*basis->P; i++)
1401-
basis->interp[i] = 1.0;
1402-
1403-
// Calculate
1404-
for (CeedInt d=0; d<basis->dim; d++)
1405-
for (CeedInt qpt=0; qpt<basis->Q; qpt++)
1406-
for (CeedInt node=0; node<basis->P; node++) {
1407-
CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d;
1408-
CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d;
1409-
basis->interp[qpt*(basis->P)+node] *= basis->interp_1d[q*basis->P_1d+p];
1410-
}
1411-
}
1412-
*interp = basis->interp;
1413-
break;
1414-
case 2: // H(div) discretization
1415-
*interp = basis->interp;
1416-
break;
1417-
case 0: // L2 discretization
1418-
*interp = basis->interp;
1419-
break;
1392+
if (!basis->interp && basis->tensor_basis) {
1393+
// Allocate
1394+
int ierr;
1395+
ierr = CeedMalloc(basis->Q*basis->P, &basis->interp); CeedChk(ierr);
1396+
1397+
// Initialize
1398+
for (CeedInt i=0; i<basis->Q*basis->P; i++)
1399+
basis->interp[i] = 1.0;
1400+
1401+
// Calculate
1402+
for (CeedInt d=0; d<basis->dim; d++)
1403+
for (CeedInt qpt=0; qpt<basis->Q; qpt++)
1404+
for (CeedInt node=0; node<basis->P; node++) {
1405+
CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d;
1406+
CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d;
1407+
basis->interp[qpt*(basis->P)+node] *= basis->interp_1d[q*basis->P_1d+p];
1408+
}
14201409
}
1410+
*interp = basis->interp;
14211411
return CEED_ERROR_SUCCESS;
14221412
}
14231413

‎interface/ceed-operator.c‎

Lines changed: 2 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -88,20 +88,16 @@ static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field,
8888
// LCOV_EXCL_STOP
8989

9090
}
91-
bool is_oriented;
92-
ierr = CeedElemRestrictionIsOriented(r, &is_oriented); CeedChk(ierr);
93-
CeedInt scale_r;
94-
scale_r = is_oriented ? qf_field->size : 1;
9591
// Field size
9692
switch(eval_mode) {
9793
case CEED_EVAL_NONE:
98-
if (size != restr_num_comp*scale_r)
94+
if (size != restr_num_comp)
9995
// LCOV_EXCL_START
10096
return CeedError(ceed, CEED_ERROR_DIMENSION,
10197
"Field '%s' of size %" CeedInt_FMT " and EvalMode %s: ElemRestriction has "
10298
CeedInt_FMT " components",
10399
qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode],
104-
restr_num_comp*scale_r);
100+
restr_num_comp);
105101
// LCOV_EXCL_STOP
106102
break;
107103
case CEED_EVAL_INTERP:

0 commit comments

Comments
 (0)