/** 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](./boundarylayer.h) 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. ~~~gnuplot 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](./../metric_amr_with_constraint/convergence_exp_restricted1.c) */ 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); } /** ~~~gnuplot Example of mesh which minimizes the interpolation error 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 ~~~ */