Skip to content

Commit

Permalink
Update ellip test to solve p-Bratu instead of plain p-Laplacian.
Browse files Browse the repository at this point in the history
Make ellip tests use -snes_monitor_short.
  • Loading branch information
jedbrown committed Sep 17, 2009
1 parent a52d383 commit 878c4fd
Show file tree
Hide file tree
Showing 9 changed files with 359 additions and 266 deletions.
7 changes: 4 additions & 3 deletions src/fs/tests/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,9 @@ dohp_add_test (fs-mf-op-proj-1 1 fs-ex1 -const_BDeg 7 -snes_mf_operator -ksp_co

dohp_add_test (bu-0 1 bu -dmesh_in "${Dohp_DATA_DIR}/dblock211.h5m" -dfs_ordering_type natural)

dohp_add_test (ellip-e0-b2-p18 1 ellip -dmesh_in "${Dohp_DATA_DIR}/dblock4.h5m" -snes_monitor -ksp_monitor -pc_type lu -exact_a 1.6 -exact_b 1.7 -exact_c 1.8 -exact 0 -const_bdeg 2 -ellip_p 1.6 -dfs_ordering_type natural)
dohp_add_test (ellip-e0-b4-p18 1 ellip -dmesh_in "${Dohp_DATA_DIR}/dblock4.h5m" -snes_monitor -ksp_monitor -pc_type lu -exact_a 1.6 -exact_b 1.7 -exact_c 1.8 -exact 0 -const_bdeg 4 -ellip_p 1.6 -djac_tensor_no_unroll -dfs_ordering_type natural)
dohp_add_test (ellip-e0-b4-p17-sse 1 ellip -dmesh_in "${Dohp_DATA_DIR}/dblock4.h5m" -snes_monitor -ksp_monitor -pc_type lu -exact_a 1.6 -exact_b 1.7 -exact_c 1.8 -exact 0 -const_bdeg 4 -ellip_p 1.4 -dfs_ordering_type natural)
dohp_add_test (ellip-e0-b2-p16 1 ellip -dmesh_in "${Dohp_DATA_DIR}/dblock4.h5m" -snes_monitor_short -ksp_monitor_short -pc_type lu -exact_a 1.6 -exact_b 1.7 -exact_c 1.8 -exact 0 -const_bdeg 2 -ellip_p 1.6 -dfs_ordering_type natural)
dohp_add_test (ellip-e0-b4-p16 1 ellip -dmesh_in "${Dohp_DATA_DIR}/dblock4.h5m" -snes_monitor_short -ksp_monitor_short -pc_type lu -exact_a 1.6 -exact_b 1.7 -exact_c 1.8 -exact 0 -const_bdeg 4 -ellip_p 1.6 -djac_tensor_no_unroll -dfs_ordering_type natural)
dohp_add_test (ellip-e0-b4-p14-sse 1 ellip -dmesh_in "${Dohp_DATA_DIR}/dblock4.h5m" -snes_monitor_short -ksp_monitor_short -pc_type lu -exact_a 1.6 -exact_b 1.7 -exact_c 1.8 -exact 0 -const_bdeg 4 -ellip_p 1.4)
dohp_add_test (ellip-e0-b4-p17-lam1-sse 1 ellip -dmesh_in "${Dohp_DATA_DIR}/dblock4.h5m" -snes_monitor_short -ksp_monitor_short -pc_type lu -exact 0 -const_bdeg 4 -ellip_p 1.7 -ellip_lam 1)
dohp_add_test (ellip-morph 1 ellip -dmesh_in "${Dohp_DATA_DIR}/dblock4.h5m" -morph 1 -twist 0.4 -ksp_monitor_short -snes_max_it 1 -ksp_rtol 1e-8 -ksp_type cg -pc_type lu)
dohp_add_test (elast-e0-b4-g01 1 elast -dmesh_in "${Dohp_DATA_DIR}/dblock4.h5m" -snes_monitor_short -ksp_monitor_short -pc_type ilu -exact_a 0.98 -exact_b 1.01 -exact_c 1.01 -exact 0 -const_bdeg 4 -elast_gamma 0.1)
52 changes: 21 additions & 31 deletions src/fs/tests/ellip.c
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
static const char help[] = "Solve a scalar elliptic problem, a regularized p-Laplacian using dual order hp elements.\n"
static const char help[] = "Solve a scalar elliptic problem, a regularized p-Bratu using dual order hp elements.\n"
"The model problem is\n"
" \\int_\\Omega (\\eta Dv \\cdot Du - f v) - \\int_\\Gamma v (\\eta Du \\cdot n) = 0\n"
" \\int_\\Omega (\\eta Dv \\cdot Du - \\lambda e^u - f v) - \\int_\\Gamma v (\\eta Du \\cdot n) = 0\n"
"where\n"
" \\eta(u) = (\\epsilon + 1/2 Du . Du)^(p-2)\n"
" (\\eta Du \\cdot n) = known OR function of u OR self (\"No condition\" outflow)\n\n";
Expand All @@ -23,6 +23,7 @@ typedef struct {
struct EllipParam {
dReal epsilon;
dReal p;
dReal lambda;
dTruth onlyproject;
dTruth bdy100;
};
Expand Down Expand Up @@ -99,8 +100,9 @@ static void EllipExact_3_Forcing(const struct EllipExactCtx *ctx,const struct El
}

struct EllipStore {
dReal eta,deta;
dReal Du[3];
dReal eta;
dReal sqrt_mdeta_Du[3];
dReal lambda_exp_u;
};

struct EllipCtx {
Expand Down Expand Up @@ -136,6 +138,7 @@ static dErr EllipCreate(MPI_Comm comm,Ellip *ellip)
prm = &elp->param;
prm->p = 2.0; /* p in p-Laplacian */
prm->epsilon = 1.0;
prm->lambda = 0.0; /* Bratu nonlinearity */
prm->onlyproject = dFALSE;

*ellip = elp;
Expand Down Expand Up @@ -175,6 +178,7 @@ static dErr EllipSetFromOptions(Ellip elp)
err = PetscOptionsTruth("-eta_monitor","Monitor nonlinearity","",elp->eta_monitor,&elp->eta_monitor,NULL);dCHK(err);
err = PetscOptionsReal("-ellip_p","p in p-Laplacian","",prm->p,&prm->p,NULL);dCHK(err);
err = PetscOptionsReal("-ellip_eps","Regularization in p-Laplacian","",prm->epsilon,&prm->epsilon,NULL);dCHK(err);
err = PetscOptionsReal("-ellip_lam","Strength of Bratu nonlinearity","",prm->lambda,&prm->lambda,NULL);dCHK(err);
err = PetscOptionsTruth("-onlyproject","Actually just do a projection","",prm->onlyproject,&prm->onlyproject,NULL);dCHK(err);
err = PetscOptionsTruth("-bdy100","Only use boundary 100","",prm->bdy100,&prm->bdy100,NULL);dCHK(err);
err = PetscOptionsInt("-exact","Exact solution choice","",exact,&exact,NULL);dCHK(err);
Expand Down Expand Up @@ -285,12 +289,13 @@ static dErr EllipDestroy(Ellip elp)

static inline void EllipPointwiseComputeStore(struct EllipParam *prm,const dReal dUNUSED x[3],const dScalar dUNUSED u[1],const dScalar Du[3],struct EllipStore *st)
{
dReal gamma,espg,p=prm->p;
dReal gamma,espg,p=prm->p,sqrt_mdeta;
gamma = 0.5 * (dSqr(Du[0]) + dSqr(Du[1]) + dSqr(Du[2]));
espg = dSqr(prm->epsilon) + gamma;
st->eta = pow(espg,(p-2)/2);
st->deta = ((p-2)/2) * st->eta / espg;
st->Du[0] = Du[0]; st->Du[1] = Du[1]; st->Du[2] = Du[2];
sqrt_mdeta = sqrt(-((p-2)/2) * st->eta / espg);
for (dInt i=0; i<3; i++) st->sqrt_mdeta_Du[i] = sqrt_mdeta*Du[i];
st->lambda_exp_u = prm->lambda * exp(u[0]);
}

static inline void EllipPointwiseFunction(struct EllipParam *prm,struct EllipExact *exact,struct EllipExactCtx *exactctx,
Expand All @@ -305,26 +310,25 @@ static inline void EllipPointwiseFunction(struct EllipParam *prm,struct EllipExa
v[0] += weight * u[0];
Dv[0] = Dv[1] = Dv[2] = 0;
} else {
v[0] += -weight * st->lambda_exp_u;
for (dInt i=0; i<3; i++) {
Dv[i] = weight * st->eta * Du[i]; /* Coefficient of Dv in weak form */
}
}
}
}

static inline void EllipPointwiseJacobian(struct EllipParam dUNUSED *prm,const struct EllipStore *restrict st,dReal weight,
const dScalar dUNUSED u[restrict static 1],const dScalar Du[restrict static 3],
const dScalar u[restrict static 1],const dScalar Du[restrict static 3],
dScalar v[restrict static 1],dScalar Dv[restrict static 3])
{
const dScalar dot = dDotScalar3(st->Du,Du);
const dReal etaw = st->eta*weight,dotdetaw = dot*st->deta*weight;
const dScalar dotw = dDotScalar3(st->sqrt_mdeta_Du,Du)*weight;
const dReal etaw = st->eta*weight;
if (prm->onlyproject) {
v[0] = weight * u[0];
Dv[0] = Dv[1] = Dv[2] = 0;
} else {
v[0] = 0;
Dv[0] = etaw*Du[0] + dotdetaw*st->Du[0];
Dv[1] = etaw*Du[1] + dotdetaw*st->Du[1];
Dv[2] = etaw*Du[2] + dotdetaw*st->Du[2];
v[0] = -weight * st->lambda_exp_u * u[0];
for (dInt i=0; i<3; i++) Dv[i] = etaw*Du[i] - dotw*st->sqrt_mdeta_Du[i];
}
}

Expand Down Expand Up @@ -363,20 +367,6 @@ static dErr EllipFunction(SNES dUNUSED snes,Vec gx,Vec gy,void *ctx)
}
err = dFSRestoreWorkspace(fs,__func__,&q,&jinv,&jw,&u,&v,&du,&dv);dCHK(err);
err = dFSRestoreElements(fs,&n,&off,&rule,&efs,&geomoff,&geom);dCHK(err);
#if 0
if (0) {
dMeshESH *sets;
dInt nsets;
err = dFSGetBoundarySets(fs,100,&nsets,&sets);dCHK(err);
for (dInt s=0; s<nsets; s++) {
err = dFSGetElements(fs,sets[s],&n,&off,&rule,&efs,&geomoff,&geom);dCHK(err);
for (dInt e=0; e<n; e++) {
// integrate weak forms over faces
}
err = dFSRestoreElements(fs,sets[s],&n,&off,&rule,&efs,&geomoff,&geom);dCHK(err);
}
}
#endif
err = VecRestoreArray(elp->x,&x);dCHK(err);
err = VecRestoreArray(elp->y,&y);dCHK(err);
err = dFSExpandedToGlobal(fs,elp->y,gy,dFS_INHOMOGENEOUS,INSERT_VALUES);dCHK(err);
Expand Down Expand Up @@ -518,12 +508,12 @@ static dErr EllipJacobian(SNES dUNUSED snes,Vec gx,Mat *J,Mat *Jp,MatStructure *
dReal v[1],Dv[3];
EllipPointwiseJacobian(&elp->param,&st,jw[lq],u,Du,v,Dv);
# if Q1SCALE
K[ltest][lp] += //basis[lq][ltest] * v[0]
K[ltest][lp] += basis[lq][ltest] * v[0] +
lscale[ltest]*lscale[lp]*(+ deriv[lq][ltest][0] * Dv[0]
+ deriv[lq][ltest][1] * Dv[1]
+ deriv[lq][ltest][2] * Dv[2]);
# else
K[ltest][lp] += //basis[lq][ltest] * v[0]
K[ltest][lp] += basis[lq][ltest] * v[0] +
(+ deriv[lq][ltest][0] * Dv[0]
+ deriv[lq][ltest][1] * Dv[1]
+ deriv[lq][ltest][2] * Dv[2]);
Expand Down
17 changes: 17 additions & 0 deletions src/fs/tests/refout/ellip-e0-b2-p16.refout
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
0 SNES Function norm 3.17818
0 KSP Residual norm 3.75851
1 KSP Residual norm < 1.e-11
1 SNES Function norm 0.152157
0 KSP Residual norm 0.199466
1 KSP Residual norm < 1.e-11
2 SNES Function norm 0.000928224
0 KSP Residual norm 0.00116026
1 KSP Residual norm < 1.e-11
3 SNES Function norm 4.24861e-08
0 KSP Residual norm 5.05985e-08
1 KSP Residual norm < 1.e-11
4 SNES Function norm < 1.e-11
Algebraic residual |x|_1 9.35e-16 |x|_2 2.59e-16 |x|_inf 1.67e-16
Interpolation residual |x|_1 2.43e-01 |x|_2 6.74e-02 |x|_inf 3.37e-02
Pointwise solution error |x|_1 3.80e-01 |x|_2 2.13e-01 |x|_inf 3.03e-01
Pointwise gradient error |x|_1 8.22e+00 |x|_2 7.19e+00 |x|_inf 1.59e+01
17 changes: 0 additions & 17 deletions src/fs/tests/refout/ellip-e0-b2-p18.refout

This file was deleted.

109 changes: 109 additions & 0 deletions src/fs/tests/refout/ellip-e0-b4-p14-sse.refout
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
0 SNES Function norm 3.08243
0 KSP Residual norm 49.3455
1 KSP Residual norm 7.44349
2 KSP Residual norm 3.97315
3 KSP Residual norm 2.15976
4 KSP Residual norm 1.20451
5 KSP Residual norm 0.698012
6 KSP Residual norm 0.39404
7 KSP Residual norm 0.183079
8 KSP Residual norm 0.098371
9 KSP Residual norm 0.0557577
10 KSP Residual norm 0.0296621
11 KSP Residual norm 0.0162693
12 KSP Residual norm 0.00891964
13 KSP Residual norm 0.00537996
14 KSP Residual norm 0.00337431
15 KSP Residual norm 0.00210799
16 KSP Residual norm 0.00127631
17 KSP Residual norm 0.000769679
18 KSP Residual norm 0.000454002
1 SNES Function norm 2.36119
0 KSP Residual norm 42.3983
1 KSP Residual norm 13.982
2 KSP Residual norm 5.18871
3 KSP Residual norm 2.46515
4 KSP Residual norm 0.96147
5 KSP Residual norm 0.386276
6 KSP Residual norm 0.129491
7 KSP Residual norm 0.0630208
8 KSP Residual norm 0.0244245
9 KSP Residual norm 0.0103638
10 KSP Residual norm 0.00498037
11 KSP Residual norm 0.00255219
12 KSP Residual norm 0.00107896
13 KSP Residual norm 0.00052551
14 KSP Residual norm 0.000263445
2 SNES Function norm 0.57823
0 KSP Residual norm 4.07038
1 KSP Residual norm 1.79894
2 KSP Residual norm 0.744108
3 KSP Residual norm 0.310638
4 KSP Residual norm 0.129062
5 KSP Residual norm 0.0560168
6 KSP Residual norm 0.0283736
7 KSP Residual norm 0.0140375
8 KSP Residual norm 0.00759276
9 KSP Residual norm 0.0034378
10 KSP Residual norm 0.00116475
11 KSP Residual norm 0.000500336
12 KSP Residual norm 0.000237274
13 KSP Residual norm 0.000110793
14 KSP Residual norm 5.4996e-05
15 KSP Residual norm 2.93445e-05
3 SNES Function norm 0.0557488
0 KSP Residual norm 0.268575
1 KSP Residual norm 0.115653
2 KSP Residual norm 0.0447676
3 KSP Residual norm 0.0202612
4 KSP Residual norm 0.00819633
5 KSP Residual norm 0.00343241
6 KSP Residual norm 0.00155376
7 KSP Residual norm 0.000784417
8 KSP Residual norm 0.000423942
9 KSP Residual norm 0.00021626
10 KSP Residual norm 6.97143e-05
11 KSP Residual norm 3.04921e-05
12 KSP Residual norm 1.24375e-05
13 KSP Residual norm 5.88531e-06
14 KSP Residual norm 3.1101e-06
15 KSP Residual norm 1.65415e-06
4 SNES Function norm 0.000346418
0 KSP Residual norm 0.00186954
1 KSP Residual norm 0.000836923
2 KSP Residual norm 0.000322917
3 KSP Residual norm 0.000140581
4 KSP Residual norm 5.06626e-05
5 KSP Residual norm 1.9646e-05
6 KSP Residual norm 8.13219e-06
7 KSP Residual norm 4.27875e-06
8 KSP Residual norm 2.17908e-06
9 KSP Residual norm 1.13513e-06
10 KSP Residual norm 4.3265e-07
11 KSP Residual norm 2.06638e-07
12 KSP Residual norm 8.7329e-08
13 KSP Residual norm 3.74862e-08
14 KSP Residual norm 1.75962e-08
5 SNES Function norm 2.60663e-08
0 KSP Residual norm 1.37492e-07
1 KSP Residual norm 6.25593e-08
2 KSP Residual norm 2.55046e-08
3 KSP Residual norm 1.19925e-08
4 KSP Residual norm 4.81454e-09
5 KSP Residual norm 2.4683e-09
6 KSP Residual norm 1.37205e-09
7 KSP Residual norm 7.914e-10
8 KSP Residual norm 3.593e-10
9 KSP Residual norm 1.500e-10
10 KSP Residual norm 5.967e-11
11 KSP Residual norm 3.153e-11
12 KSP Residual norm 1.643e-11
13 KSP Residual norm < 1.e-11
14 KSP Residual norm < 1.e-11
15 KSP Residual norm < 1.e-11
16 KSP Residual norm < 1.e-11
6 SNES Function norm < 1.e-11
Algebraic residual |x|_1 4.85e-12 |x|_2 1.95e-13 |x|_inf 4.73e-14
Interpolation residual |x|_1 4.90e-02 |x|_2 2.56e-03 |x|_inf 4.11e-04
Pointwise solution error |x|_1 1.14e-03 |x|_2 7.55e-04 |x|_inf 2.68e-03
Pointwise gradient error |x|_1 4.85e+00 |x|_2 5.11e+00 |x|_inf 1.70e+01
106 changes: 106 additions & 0 deletions src/fs/tests/refout/ellip-e0-b4-p16.refout
Original file line number Diff line number Diff line change
@@ -0,0 +1,106 @@
0 SNES Function norm 4.93408
0 KSP Residual norm 51.7059
1 KSP Residual norm 6.79209
2 KSP Residual norm 2.7863
3 KSP Residual norm 1.50687
4 KSP Residual norm 0.747774
5 KSP Residual norm 0.374842
6 KSP Residual norm 0.174141
7 KSP Residual norm 0.0784086
8 KSP Residual norm 0.0389488
9 KSP Residual norm 0.0193889
10 KSP Residual norm 0.0100507
11 KSP Residual norm 0.00610852
12 KSP Residual norm 0.00343091
13 KSP Residual norm 0.00190266
14 KSP Residual norm 0.00102771
15 KSP Residual norm 0.000539639
16 KSP Residual norm 0.000292823
1 SNES Function norm 2.62495
0 KSP Residual norm 16.3713
1 KSP Residual norm 6.34626
2 KSP Residual norm 2.33583
3 KSP Residual norm 0.998847
4 KSP Residual norm 0.392279
5 KSP Residual norm 0.157428
6 KSP Residual norm 0.0496749
7 KSP Residual norm 0.0206766
8 KSP Residual norm 0.0114255
9 KSP Residual norm 0.00679726
10 KSP Residual norm 0.00294285
11 KSP Residual norm 0.0012573
12 KSP Residual norm 0.000521005
13 KSP Residual norm 0.000213325
14 KSP Residual norm 9.09341e-05
2 SNES Function norm 1.08788
0 KSP Residual norm 5.39009
1 KSP Residual norm 1.80249
2 KSP Residual norm 0.628344
3 KSP Residual norm 0.265499
4 KSP Residual norm 0.100997
5 KSP Residual norm 0.0357152
6 KSP Residual norm 0.0124749
7 KSP Residual norm 0.0047978
8 KSP Residual norm 0.0022914
9 KSP Residual norm 0.00125394
10 KSP Residual norm 0.000622427
11 KSP Residual norm 0.000285466
12 KSP Residual norm 0.000142561
13 KSP Residual norm 6.54096e-05
14 KSP Residual norm 2.795e-05
3 SNES Function norm 0.0756234
0 KSP Residual norm 0.461806
1 KSP Residual norm 0.225487
2 KSP Residual norm 0.0890375
3 KSP Residual norm 0.0359211
4 KSP Residual norm 0.0121965
5 KSP Residual norm 0.00429802
6 KSP Residual norm 0.00196438
7 KSP Residual norm 0.00102016
8 KSP Residual norm 0.000570907
9 KSP Residual norm 0.000282919
10 KSP Residual norm 8.64812e-05
11 KSP Residual norm 4.16991e-05
12 KSP Residual norm 2.14124e-05
13 KSP Residual norm 1.0569e-05
14 KSP Residual norm 5.43069e-06
15 KSP Residual norm 2.91421e-06
4 SNES Function norm 0.000462468
0 KSP Residual norm 0.00247013
1 KSP Residual norm 0.00124566
2 KSP Residual norm 0.000471243
3 KSP Residual norm 0.000174178
4 KSP Residual norm 6.61506e-05
5 KSP Residual norm 2.40872e-05
6 KSP Residual norm 1.01713e-05
7 KSP Residual norm 4.01465e-06
8 KSP Residual norm 1.86325e-06
9 KSP Residual norm 1.03734e-06
10 KSP Residual norm 5.39176e-07
11 KSP Residual norm 2.78381e-07
12 KSP Residual norm 1.11244e-07
13 KSP Residual norm 4.42854e-08
14 KSP Residual norm 2.04564e-08
5 SNES Function norm 3.33707e-08
0 KSP Residual norm 1.66742e-07
1 KSP Residual norm 7.62914e-08
2 KSP Residual norm 3.18992e-08
3 KSP Residual norm 1.33728e-08
4 KSP Residual norm 5.28412e-09
5 KSP Residual norm 2.70809e-09
6 KSP Residual norm 1.68188e-09
7 KSP Residual norm 9.450e-10
8 KSP Residual norm 4.250e-10
9 KSP Residual norm 1.778e-10
10 KSP Residual norm 7.143e-11
11 KSP Residual norm 4.171e-11
12 KSP Residual norm 2.115e-11
13 KSP Residual norm 1.118e-11
14 KSP Residual norm < 1.e-11
15 KSP Residual norm < 1.e-11
16 KSP Residual norm < 1.e-11
6 SNES Function norm < 1.e-11
Algebraic residual |x|_1 6.18e-12 |x|_2 2.59e-13 |x|_inf 4.08e-14
Interpolation residual |x|_1 3.13e-02 |x|_2 1.57e-03 |x|_inf 2.55e-04
Pointwise solution error |x|_1 1.13e-03 |x|_2 7.50e-04 |x|_inf 2.66e-03
Pointwise gradient error |x|_1 4.85e+00 |x|_2 5.11e+00 |x|_inf 1.70e+01
Loading

0 comments on commit 878c4fd

Please sign in to comment.