sandbox/prouvost/AMR_examples/no_hmin/sharp.c
We search an adapted mesh containing N_{obj} elements which minimizes the total error (error between the numerical solution of the Basilisk Poisson-Helmholtz solver) and the (known) analytical solution.
The solution is an boundary layer-like function.
#include "poisson.h" // solver
#include "./boundarylayer.h" // ''problem''
#include "./../../AMR_tools/amr.h" // AMR
#include "./../../utils/gauss_quadrature.h" // compute total error
#include "utils.h"
double TOL = 1.e-7;
const face vector alp[] = {D,D};
const scalar lam[] =s;
scalar psi[]; // the numerical solution
psi[left] = dirichlet(exact(x,y,0.));
psi[right] = dirichlet(exact(x,y,0.));
psi[top] = dirichlet(exact(x,y,0.));
psi[bottom] = dirichlet(exact(x,y,0.));
scalar psi_exact[]; // the analytical solution
psi_exact[left] = dirichlet(exact(x,y,0.));
psi_exact[right] = dirichlet(exact(x,y,0.));
psi_exact[top] = dirichlet(exact(x,y,0.));
psi_exact[bottom] = dirichlet(exact(x,y,0.));
int main() {
FILE * fpglobal = fopen("error","w");
L0=1.;
int mylev=6;
int Nobj = pow(2,mylev*2); // initial objective number of element
init_grid (1 << (mylev-1));
scalar rhs[]; // RHS term for Poisson-Helmholtz equation
We do a loop to obtain several adapted meshes having N_{obj} elements.
For each N_{obj}, we do a loop to obtain the objective number of elements.
for (int j=0; Nobj <= pow(2,2*8)/1.5; j++){
for (int ki=0; ( fabs((double)(grid->tn - Nobj)/Nobj)) > 0.03 ; ki++){
foreach ()
rhs[] = src(x,y,0);
poisson (psi, rhs, alp, lam, tolerance = TOL); // poisson solver
We update an epsilon criterion to obtain the objective number of elements
// AMR criterion
if (ki < 5 && j == 0) { // estimate AMReps the first run with uniform refinement
struct PreFactorData cd = compute_prefactors (2, {psi} );
AMReps = cd.cuniform*pow(Nobj,-1)/pow(Nobj,1./2.);
} else {
AMReps = update_epsilon_control(Nobj);
}
We adapt the mesh.
adapt_metric( {psi} );
/* astats st = adapt_metric( {psi} ); */
}
We compute total error and interpolation error.
foreach ()
rhs[] = src(x,y,0);
poisson (psi, rhs, alp, lam, tolerance = TOL);
// TOTAL ERROR
double errtot = norm_gauss_5p (psi, user_norm, exact); // compute total error : ||u_num - u_exact||
// INTERPOLATION ERROR estimate with metric-based formulae
double Interr = Interpolation_error(user_norm);
// INTERPOLATION ERROR exact
foreach ()
psi_exact[] = exact(x,y,0);
double Int_err_exact = norm_gauss_5p (psi_exact, user_norm, exact); // compute exact interpolation error : ||u_interp - u_exact||
// theoretical optimal and uniform errors
struct PreFactorData cd;
cd = compute_prefactors (2, {psi}, NULL );
fprintf(fpglobal,"%ld %g %.10g %.10g %.10g %g\n", grid->tn, errtot, Int_err_exact, Interr, cd.copt*pow(grid->tn,-1), cd.cuniform*pow(grid->tn,-1));
We plot the total error and the interpolation error. We see that the measured total error and the interpolation error are close to the optimal. Thus, we show that, in this case, the obtained meshes minimizes the interpolation error and the total error.
reset
set term pngcairo enhanced size 500,500
set output 'error.png'
set logscale
set xtics (16,32,64,128,256)
set format y "10^{%T}"
set xrange [64:190]
p "error" u (sqrt($1)):6 w l t "uniform",\
"error" u (sqrt($1)):5 w l t "optimal",\
"error" u (sqrt($1)):2 w p t "total error",\
"error" u (sqrt($1)):3 w p t "interp error"
We verify that we obtained the expected result
assert ( fabs((double)(grid->tn - Nobj)/Nobj) <= 0.03);
if (sqrt(Nobj)>150)
assert(errtot<3.*cd.copt*pow(grid->tn,-1) );
We ouput the mesh to compare with the computation with constraint
scalar lev[];
foreach()
lev[]=level;
char fname[80];
sprintf(fname,"mesh_%i",Nobj);
FILE * fp = fopen(fname,"w");
output_field ({lev}, fp, n = 1 << 8, linear=false);
fclose(fp);
loop update
Nobj *= 1.5;
}
free_grid();
fclose(fpglobal);
}
reset
set term pngcairo enhanced size 350,350 font ",8"
set output 'mesh.png'
# palette turbo
set palette defined (\
0.0 0.18995 0.07176 0.23217,\
0.00390625 0.19483 0.08339 0.26149,\
0.0078125 0.19956 0.09498 0.29024,\
0.01171875 0.20415 0.10652 0.31844,\
0.015625 0.2086 0.11802 0.34607,\
0.01953125 0.21291 0.12947 0.37314,\
0.0234375 0.21708 0.14087 0.39964,\
0.02734375 0.22111 0.15223 0.42558,\
0.03125 0.225 0.16354 0.45096,\
0.03515625 0.22875 0.17481 0.47578,\
0.0390625 0.23236 0.18603 0.50004,\
0.04296875 0.23582 0.1972 0.52373,\
0.046875 0.23915 0.20833 0.54686,\
0.05078125 0.24234 0.21941 0.56942,\
0.0546875 0.24539 0.23044 0.59142,\
0.05859375 0.2483 0.24143 0.61286,\
0.0625 0.25107 0.25237 0.63374,\
0.06640625 0.25369 0.26327 0.65406,\
0.0703125 0.25618 0.27412 0.67381,\
0.07421875 0.25853 0.28492 0.693,\
0.078125 0.26074 0.29568 0.71162,\
0.08203125 0.2628 0.30639 0.72968,\
0.0859375 0.26473 0.31706 0.74718,\
0.08984375 0.26652 0.32768 0.76412,\
0.09375 0.26816 0.33825 0.7805,\
0.09765625 0.26967 0.34878 0.79631,\
0.1015625 0.27103 0.35926 0.81156,\
0.10546875 0.27226 0.3697 0.82624,\
0.109375 0.27334 0.38008 0.84037,\
0.11328125 0.27429 0.39043 0.85393,\
0.1171875 0.27509 0.40072 0.86692,\
0.12109375 0.27576 0.41097 0.87936,\
0.125 0.27628 0.42118 0.89123,\
0.12890625 0.27667 0.43134 0.90254,\
0.1328125 0.27691 0.44145 0.91328,\
0.13671875 0.27701 0.45152 0.92347,\
0.140625 0.27698 0.46153 0.93309,\
0.14453125 0.2768 0.47151 0.94214,\
0.1484375 0.27648 0.48144 0.95064,\
0.15234375 0.27603 0.49132 0.95857,\
0.15625 0.27543 0.50115 0.96594,\
0.16015625 0.27469 0.51094 0.97275,\
0.1640625 0.27381 0.52069 0.97899,\
0.16796875 0.27273 0.5304 0.98461,\
0.171875 0.27106 0.54015 0.9893,\
0.17578125 0.26878 0.54995 0.99303,\
0.1796875 0.26592 0.55979 0.99583,\
0.18359375 0.26252 0.56967 0.99773,\
0.1875 0.25862 0.57958 0.99876,\
0.19140625 0.25425 0.5895 0.99896,\
0.1953125 0.24946 0.59943 0.99835,\
0.19921875 0.24427 0.60937 0.99697,\
0.203125 0.23874 0.61931 0.99485,\
0.20703125 0.23288 0.62923 0.99202,\
0.2109375 0.22676 0.63913 0.98851,\
0.21484375 0.22039 0.64901 0.98436,\
0.21875 0.21382 0.65886 0.97959,\
0.22265625 0.20708 0.66866 0.97423,\
0.2265625 0.20021 0.67842 0.96833,\
0.23046875 0.19326 0.68812 0.9619,\
0.234375 0.18625 0.69775 0.95498,\
0.23828125 0.17923 0.70732 0.94761,\
0.2421875 0.17223 0.7168 0.93981,\
0.24609375 0.16529 0.7262 0.93161,\
0.25 0.15844 0.73551 0.92305,\
0.25390625 0.15173 0.74472 0.91416,\
0.2578125 0.14519 0.75381 0.90496,\
0.26171875 0.13886 0.76279 0.8955,\
0.265625 0.13278 0.77165 0.8858,\
0.26953125 0.12698 0.78037 0.8759,\
0.2734375 0.12151 0.78896 0.86581,\
0.27734375 0.11639 0.7974 0.85559,\
0.28125 0.11167 0.80569 0.84525,\
0.28515625 0.10738 0.81381 0.83484,\
0.2890625 0.10357 0.82177 0.82437,\
0.29296875 0.10026 0.82955 0.81389,\
0.296875 0.0975 0.83714 0.80342,\
0.30078125 0.09532 0.84455 0.79299,\
0.3046875 0.09377 0.85175 0.78264,\
0.30859375 0.09287 0.85875 0.7724,\
0.3125 0.09267 0.86554 0.7623,\
0.31640625 0.0932 0.87211 0.75237,\
0.3203125 0.09451 0.87844 0.74265,\
0.32421875 0.09662 0.88454 0.73316,\
0.328125 0.09958 0.8904 0.72393,\
0.33203125 0.10342 0.896 0.715,\
0.3359375 0.10815 0.90142 0.70599,\
0.33984375 0.11374 0.90673 0.69651,\
0.34375 0.12014 0.91193 0.6866,\
0.34765625 0.12733 0.91701 0.67627,\
0.3515625 0.13526 0.92197 0.66556,\
0.35546875 0.14391 0.9268 0.65448,\
0.359375 0.15323 0.93151 0.64308,\
0.36328125 0.16319 0.93609 0.63137,\
0.3671875 0.17377 0.94053 0.61938,\
0.37109375 0.18491 0.94484 0.60713,\
0.375 0.19659 0.94901 0.59466,\
0.37890625 0.20877 0.95304 0.58199,\
0.3828125 0.22142 0.95692 0.56914,\
0.38671875 0.23449 0.96065 0.55614,\
0.390625 0.24797 0.96423 0.54303,\
0.39453125 0.2618 0.96765 0.52981,\
0.3984375 0.27597 0.97092 0.51653,\
0.40234375 0.29042 0.97403 0.50321,\
0.40625 0.30513 0.97697 0.48987,\
0.41015625 0.32006 0.97974 0.47654,\
0.4140625 0.33517 0.98234 0.46325,\
0.41796875 0.35043 0.98477 0.45002,\
0.421875 0.36581 0.98702 0.43688,\
0.42578125 0.38127 0.98909 0.42386,\
0.4296875 0.39678 0.99098 0.41098,\
0.43359375 0.41229 0.99268 0.39826,\
0.4375 0.42778 0.99419 0.38575,\
0.44140625 0.44321 0.99551 0.37345,\
0.4453125 0.45854 0.99663 0.3614,\
0.44921875 0.47375 0.99755 0.34963,\
0.453125 0.48879 0.99828 0.33816,\
0.45703125 0.50362 0.99879 0.32701,\
0.4609375 0.51822 0.9991 0.31622,\
0.46484375 0.53255 0.99919 0.30581,\
0.46875 0.54658 0.99907 0.29581,\
0.47265625 0.56026 0.99873 0.28623,\
0.4765625 0.57357 0.99817 0.27712,\
0.48046875 0.58646 0.99739 0.26849,\
0.484375 0.59891 0.99638 0.26038,\
0.48828125 0.61088 0.99514 0.2528,\
0.4921875 0.62233 0.99366 0.24579,\
0.49609375 0.63323 0.99195 0.23937,\
0.5 0.64362 0.98999 0.23356,\
0.50390625 0.65394 0.98775 0.22835,\
0.5078125 0.66428 0.98524 0.2237,\
0.51171875 0.67462 0.98246 0.2196,\
0.515625 0.68494 0.97941 0.21602,\
0.51953125 0.69525 0.9761 0.21294,\
0.5234375 0.70553 0.97255 0.21032,\
0.52734375 0.71577 0.96875 0.20815,\
0.53125 0.72596 0.9647 0.2064,\
0.53515625 0.7361 0.96043 0.20504,\
0.5390625 0.74617 0.95593 0.20406,\
0.54296875 0.75617 0.95121 0.20343,\
0.546875 0.76608 0.94627 0.20311,\
0.55078125 0.77591 0.94113 0.2031,\
0.5546875 0.78563 0.93579 0.20336,\
0.55859375 0.79524 0.93025 0.20386,\
0.5625 0.80473 0.92452 0.20459,\
0.56640625 0.8141 0.91861 0.20552,\
0.5703125 0.82333 0.91253 0.20663,\
0.57421875 0.83241 0.90627 0.20788,\
0.578125 0.84133 0.89986 0.20926,\
0.58203125 0.8501 0.89328 0.21074,\
0.5859375 0.85868 0.88655 0.2123,\
0.58984375 0.86709 0.87968 0.21391,\
0.59375 0.8753 0.87267 0.21555,\
0.59765625 0.88331 0.86553 0.21719,\
0.6015625 0.89112 0.85826 0.2188,\
0.60546875 0.8987 0.85087 0.22038,\
0.609375 0.90605 0.84337 0.22188,\
0.61328125 0.91317 0.83576 0.22328,\
0.6171875 0.92004 0.82806 0.22456,\
0.62109375 0.92666 0.82025 0.2257,\
0.625 0.93301 0.81236 0.22667,\
0.62890625 0.93909 0.80439 0.22744,\
0.6328125 0.94489 0.79634 0.228,\
0.63671875 0.95039 0.78823 0.22831,\
0.640625 0.9556 0.78005 0.22836,\
0.64453125 0.96049 0.77181 0.22811,\
0.6484375 0.96507 0.76352 0.22754,\
0.65234375 0.96931 0.75519 0.22663,\
0.65625 0.97323 0.74682 0.22536,\
0.66015625 0.97679 0.73842 0.22369,\
0.6640625 0.98 0.73 0.22161,\
0.66796875 0.98289 0.7214 0.21918,\
0.671875 0.98549 0.7125 0.2165,\
0.67578125 0.98781 0.7033 0.21358,\
0.6796875 0.98986 0.69382 0.21043,\
0.68359375 0.99163 0.68408 0.20706,\
0.6875 0.99314 0.67408 0.20348,\
0.69140625 0.99438 0.66386 0.19971,\
0.6953125 0.99535 0.65341 0.19577,\
0.69921875 0.99607 0.64277 0.19165,\
0.703125 0.99654 0.63193 0.18738,\
0.70703125 0.99675 0.62093 0.18297,\
0.7109375 0.99672 0.60977 0.17842,\
0.71484375 0.99644 0.59846 0.17376,\
0.71875 0.99593 0.58703 0.16899,\
0.72265625 0.99517 0.57549 0.16412,\
0.7265625 0.99419 0.56386 0.15918,\
0.73046875 0.99297 0.55214 0.15417,\
0.734375 0.99153 0.54036 0.1491,\
0.73828125 0.98987 0.52854 0.14398,\
0.7421875 0.98799 0.51667 0.13883,\
0.74609375 0.9859 0.50479 0.13367,\
0.75 0.9836 0.49291 0.12849,\
0.75390625 0.98108 0.48104 0.12332,\
0.7578125 0.97837 0.4692 0.11817,\
0.76171875 0.97545 0.4574 0.11305,\
0.765625 0.97234 0.44565 0.10797,\
0.76953125 0.96904 0.43399 0.10294,\
0.7734375 0.96555 0.42241 0.09798,\
0.77734375 0.96187 0.41093 0.0931,\
0.78125 0.95801 0.39958 0.08831,\
0.78515625 0.95398 0.38836 0.08362,\
0.7890625 0.94977 0.37729 0.07905,\
0.79296875 0.94538 0.36638 0.07461,\
0.796875 0.94084 0.35566 0.07031,\
0.80078125 0.93612 0.34513 0.06616,\
0.8046875 0.93125 0.33482 0.06218,\
0.80859375 0.92623 0.32473 0.05837,\
0.8125 0.92105 0.31489 0.05475,\
0.81640625 0.91572 0.3053 0.05134,\
0.8203125 0.91024 0.29599 0.04814,\
0.82421875 0.90463 0.28696 0.04516,\
0.828125 0.89888 0.27824 0.04243,\
0.83203125 0.89298 0.26981 0.03993,\
0.8359375 0.88691 0.26152 0.03753,\
0.83984375 0.88066 0.25334 0.03521,\
0.84375 0.87422 0.24526 0.03297,\
0.84765625 0.8676 0.2373 0.03082,\
0.8515625 0.86079 0.22945 0.02875,\
0.85546875 0.8538 0.2217 0.02677,\
0.859375 0.84662 0.21407 0.02487,\
0.86328125 0.83926 0.20654 0.02305,\
0.8671875 0.83172 0.19912 0.02131,\
0.87109375 0.82399 0.19182 0.01966,\
0.875 0.81608 0.18462 0.01809,\
0.87890625 0.80799 0.17753 0.0166,\
0.8828125 0.79971 0.17055 0.0152,\
0.88671875 0.79125 0.16368 0.01387,\
0.890625 0.7826 0.15693 0.01264,\
0.89453125 0.77377 0.15028 0.01148,\
0.8984375 0.76476 0.14374 0.01041,\
0.90234375 0.75556 0.13731 0.00942,\
0.90625 0.74617 0.13098 0.00851,\
0.91015625 0.73661 0.12477 0.00769,\
0.9140625 0.72686 0.11867 0.00695,\
0.91796875 0.71692 0.11268 0.00629,\
0.921875 0.7068 0.1068 0.00571,\
0.92578125 0.6965 0.10102 0.00522,\
0.9296875 0.68602 0.09536 0.00481,\
0.93359375 0.67535 0.0898 0.00449,\
0.9375 0.66449 0.08436 0.00424,\
0.94140625 0.65345 0.07902 0.00408,\
0.9453125 0.64223 0.0738 0.00401,\
0.94921875 0.63082 0.06868 0.00401,\
0.953125 0.61923 0.06367 0.0041,\
0.95703125 0.60746 0.05878 0.00427,\
0.9609375 0.5955 0.05399 0.00453,\
0.96484375 0.58336 0.04931 0.00486,\
0.96875 0.57103 0.04474 0.00529,\
0.97265625 0.55852 0.04028 0.00579,\
0.9765625 0.54583 0.03593 0.00638,\
0.98046875 0.53295 0.03169 0.00705,\
0.984375 0.51989 0.02756 0.0078,\
0.98828125 0.50664 0.02354 0.00863,\
0.9921875 0.49321 0.01963 0.00955,\
0.99609375 0.4796 0.01583 0.01055)
unset key
set size ratio -1
set pm3d map
sp "mesh_13824" u 1:2:3