/** 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](./../metric_amr_no_constraint/exponentiel.h) is an exponential function. */ #include "poisson.h" // solver #include "./../no_hmin/exponentiel.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 with an additionnal constraint on the minimal size computed from an estimation of the optimal conmpression ratio. */ /* restriction on the minimum grid size */ double etaopt = estimate_eta_opt(2, {psi}); maxlevel = 0.5*log(Nobj/etaopt)/log(2.); 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 and we compare with the case without constraint. We plot the error in function of the square root of the number of element (equivalent to mean cell size in 2D). We see that the total error is lower with the constraint than without constraint, but the interpolation error is slightly higher with the constraint. ~~~gnuplot error reset set term pngcairo enhanced size 700,400 set output 'error.png' set logscale set xtics (16,32,64,128,256) set format y "10^{%T}" set xrange [50:191] set key below set multiplot layout 1,2 set title "total error" 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 with constraint",\ "./../../metric_amr_no_constraint/convergence_exp/error" u (sqrt($1)):2 w p t "total error no constraint" set title "interpolation error" p "error" u (sqrt($1)):6 w l t "uniform",\ "error" u (sqrt($1)):5 w l t "optimal",\ "error" u (sqrt($1)):3 w p t "interp error with constraint",\ "./../../metric_amr_no_constraint/convergence_exp/error" u (sqrt($1)):3 w p t "interp error no constraint" unset multiplot ~~~ */ /** We verify that we obtained the expected result */ assert ( fabs((double)(grid->tn - Nobj)/Nobj) <= 0.03 && errtot<5.*cd.copt*pow(grid->tn,-(2./3.)) ); /** 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); } /** We compare the meshes with and without constraint. They are different. ~~~gnuplot Example of mesh which minimizes the interpolation error reset set term pngcairo enhanced size 700,700 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) set cbrange [3:10] unset key set size ratio -1 set pm3d map set multiplot layout 1,2 set title "With constraint" sp "mesh_20736" u 1:2:3 set title "No constraint" sp "./../../no_hmin/exp/mesh_20736" u 1:2:3 unset multiplot ~~~ */