8 patches for repository http://basilisk.fr/basilisk: patch 0e7260828903298ca329a9fe301bb844f9ea7fad Author: Jose M. Lopez-Herrera Date: Mon Dec 21 16:20:27 CET 2020 * VOF and embed solids are made compatible. As suggested by Stephane the variable *fms* collects both metric and solid factors. This patch focus on the necessary changes to have the N-S centered solver operative. patch 76f40cac47a756c555eaa29dbf47751a58452cce Author: Jose M. Lopez-Herrera Date: Mon Dec 21 16:26:25 CET 2020 * ref files are also updated because fms. patch 118252e55b3f81609455bd0052973d9fceb2df0a Author: Jose M. Lopez-Herrera Date: Mon Dec 21 16:41:55 CET 2020 * Since now fm is fms, update is needed for several solvers. In addition some metric factors lacking in all-mach & momentum.h are completed. patch d51ae4de0b8ec5bfe7465f6ca7be153c032396e9 Author: Jose M. Lopez-Herrera Date: Mon Dec 21 17:18:34 CET 2020 * test cases and examples updated to the new *fms* variable. patch 763229fe28669f63ed0bdd79ee11b062ae537f7a Author: Jose M. Lopez-Herrera Date: Mon Dec 21 17:39:41 CET 2020 * A figure to explain better the small cell problem in vof variables patch 02938c01b965c53f4ca3f0b0b64e6ecf54177895 Author: Jose M. Lopez-Herrera Date: Sat Jan 16 17:38:08 CET 2021 * changes in "axi.h" to make compatible "embed.h" patch 8d0a95e57c508b0d46d9e7803e9fbdaf345f7402 Author: Jose M. Lopez-Herrera Date: Mon Feb 1 17:49:42 CET 2021 * additional changes to embed.h + axi.h: poisson solver works patch 2a2ec127af5bafed67c090287e2c2e38e18df3a4 Author: Jose M. Lopez-Herrera Date: Wed Feb 24 13:06:19 CET 2021 * Fixed minimal conflicts in karman.c New patches: [VOF and embed solids are made compatible. Jose M. Lopez-Herrera **20201221152027 Ignore-this: abb3af4da69f5e527ea1bdaccc067ac9 As suggested by Stephane the variable *fms* collects both metric and solid factors. This patch focus on the necessary changes to have the N-S centered solver operative. ] hunk ./src/README 219 -![Metric scale factors](figures/metric.svg) - -The [face vector](/Basilisk C#face-and-vertex-fields) *fm* is the +![Metric scale and factors](figures/metric.svg) + +The [face vector](/Basilisk C#face-and-vertex-fields) *fms* is the hunk ./src/README 223 -$fm\Delta$ and the scalar field *cm* is the scale factor for the area -of the cell i.e. the physical area is $cm\Delta^2$. By default, these +$fms\Delta$ and the scalar field *cms* is the scale factor for the area +of the cell i.e. the physical area is $cms\Delta^2$. By default, these hunk ./src/README 227 +Note, however, that if an embed solid is present, *fms* and*cms* will +include also the solid factors for sake of brevity and simplicity. + hunk ./src/bcg.h 32 - strictly be `dt*(uf.x[i] + uf.x[i+1])/((fm.x[] + - fm.x[i+1])*Delta)` but this causes trouble with boundary + strictly be `dt*(uf.x[i] + uf.x[i+1])/((fms.x[] + + fms.x[i+1])*Delta)` but this causes trouble with boundary hunk ./src/bcg.h 36 - double un = dt*uf.x[]/(fm.x[]*Delta + SEPS), s = sign(un); + double un = dt*uf.x[]/(fms.x[]*Delta + SEPS), s = sign(un); hunk ./src/bcg.h 44 - if (fm.y[i] && fm.y[i,1]) { - double vn = (uf.y[i] + uf.y[i,1])/(fm.y[i] + fm.y[i,1]); + if (fms.y[i] && fms.y[i,1]) { + double vn = (uf.y[i] + uf.y[i,1])/(fms.y[i] + fms.y[i,1]); hunk ./src/bcg.h 51 - if (fm.z[i] && fm.z[i,0,1]) { - double wn = (uf.z[i] + uf.z[i,0,1])/(fm.z[i] + fm.z[i,0,1]); + if (fms.z[i] && fms.z[i,0,1]) { + double wn = (uf.z[i] + uf.z[i,0,1])/(fms.z[i] + fms.z[i,0,1]); hunk ./src/bcg.h 101 - f[] += p.dt*(flux.x[] - flux.x[1])/(Delta*cm[]); + f[] += p.dt*(flux.x[] - flux.x[1])/(Delta*cms[]); hunk ./src/common.h 886 -# define dv() (Delta*cm[]) +# define dv() (Delta*cms[]) hunk ./src/common.h 889 -# define dv() (sq(Delta)*cm[]) +# define dv() (sq(Delta)*cms[]) hunk ./src/common.h 892 -# define dv() (cube(Delta)*cm[]) +# define dv() (cube(Delta)*cms[]) hunk ./src/common.h 1201 -// Metric - -(const) face vector fm = unityf; -(const) scalar cm = unity; +// Metric & Solid + +(const) face vector fms = unityf; +(const) scalar cms = unity; hunk ./src/common.h 1209 +(const) face vector fs = unityf; +(const) scalar cs = unity; + hunk ./src/embed.h 19 -scalar cs[]; -face vector fs[]; +extern scalar cs; +extern face vector fs; hunk ./src/embed.h 828 - already stores the product $f_su$. */ + already stores the product $u_f \, f_{ms}$. */ hunk ./src/embed.h 836 - double dtmax = Delta*cs[]/(umax + SEPS); + double dtmax = Delta*cms[]/(umax + SEPS); hunk ./src/embed.h 844 - F /= Delta*cs[]; + F /= Delta*cms[]; hunk ./src/embed.h 866 - scs += sq(cs[]); - e[] = (dt - dtmax)*F*cs[]/scs; + scs += sq(cms[]); //Correct? + e[] = (dt - dtmax)*F*cms[]/scs; hunk ./src/embed.h 886 +### The advection of VOF-like variables + +In the explicit advection of VOF-like variables, the timestep +could be also limited by the CFL condition when applied to cells of small +volume. The function below performs the merging technique for cells +susceptible of "overflow". This function is an adaption of [sweep_x +()](vof.h#one-dimensional-advection) where additional comments can +be found.*/ + +#define FILL 0 +foreach_dimension() +void vof_update_x (scalar c, scalar cc, scalar flux, + scalar * tracers, scalar * tcl, scalar * tfluxl, + face vector uf, double dt) +{ + /** + Only for cells containing fluid, i.e. $c_s \neq 0$, the vof-like + variables are updated.*/ + + foreach() { + if (cs[] > 0.) { + + /** + Full cells can not overflow. The VOF-like variables are advanced + in time.*/ + + if (cs[] >= 1.) { + c[] += dt*(flux[] - flux[1] + + cc[]*(uf.x[1] - uf.x[]))/(cms[]*Delta + SEPS); + scalar t, tc, tflux; + for (t, tc, tflux in tracers, tcl, tfluxl) + t[] += dt*(tflux[] - tflux[1] + + tc[]*(uf.x[1] - uf.x[]))/(cms[]*Delta + SEPS); + } + + /** + Only in mixed cell (cells with fluid and solid) has sense to + determine *dtmax* and check if the overflow exists. */ + + else { + double umax = 0.; + for (int i = 0; i <= 1; i++) { + if (fabs(uf.x[i]) > umax) + umax = fabs(uf.x[i]); + } + double dtmax = Delta*cms[]/(umax + SEPS); + + if ( dt <= dtmax) { + c[] += dt*(flux[] - flux[1] + + cc[]*(uf.x[1] - uf.x[]))/(cms[]*Delta + SEPS); + + /** + In mixed cells, if the room available for the fluid left by + the solid is smaller that the fraction of fluid in that + cell, we set $c = 1$ since the mixed cell must be full of + fluid. */ + +#if FILL + c[] = (cs[] < c[] ? 1.0 : c[]); +#endif + scalar t, tc, tflux; + for (t, tc, tflux in tracers, tcl, tfluxl) + t[] += dt*(tflux[] - tflux[1] + + tc[]*(uf.x[1] - uf.x[]))/(cms[]*Delta + SEPS); + } + + /** + If an overflow exists, we will consider the merging of the + small cell with the contiguous one to which the relevant + normal component points out. + + ![Merging cells](figures/merged.svg) + + Assuming that *cc*[merged] is equal to *cc*[1], the excess + that goes to the adjacent cell is equal to: + + $$ + E = (\Delta t^{max} - \Delta t) (F_x[1] - F_x[]) + (u_f.x[] + -u_f.x[1])(cc[] \Delta t^{max} -cc[1] \Delta t) + $$ + */ + + else { + fprintf (ferr, "WARNING: Small cell found (dt/dtmax = %g)\n", + dt/dtmax), fflush (ferr); + int i = (fs.x[] >= fs.x[1] ? -1 : 1); + c[] += dtmax*(flux[] - flux[1] + + cc[]*(uf.x[1] - uf.x[]))/(cms[]*Delta + SEPS); + c[i] += ((dt -dtmax)*(flux[] - flux[i]) + + (uf.x[1] - uf.x[])*(cc[]*dtmax -cc[i]*dt))/(cms[i]*Delta + SEPS); +#if FILL + c[] = (cs[] < c[] ? 1.0 : c[]); + c[i] = (cs[i] < c[i] ? 1.0 : c[i]); +#endif + scalar t, tc, tflux; + for (t, tc, tflux in tracers, tcl, tfluxl) { + t[] += dtmax*(tflux[] - tflux[1] + + tc[]*(uf.x[1] - uf.x[]))/(cms[]*Delta + SEPS); + t[i] += ((dt-dtmax)*(tflux[] - tflux[i]) + + (uf.x[1] - uf.x[])*(tc[]*dtmax -tc[i]*dt))/(cms[i]*Delta + SEPS); + } + } + } + } + } + } + +/** hunk ./src/embed.h 997 -the metric using the embedded fractions. */ +the metric using the embedded fractions. The variables are not longer +constant fields and they must be allocated. Also, the solid fields +are moved to the beginning of the list of variables. This ensures that +solid fields are the first to be updated. */ hunk ./src/embed.h 1004 + if (is_constant(fs.x)) { + scalar * l = list_copy (all); + fs = new face vector; + free (all); + all = list_concat ((scalar *){fs}, l); + free (l); + } + + face vector fsv = fs; + foreach_face() + fsv.x[] = 1.; + + if (is_constant(cs)) { + scalar * l = list_copy (all); + cs = new scalar; + free (all); + all = list_concat ({cs}, l); + free (l); + } + + scalar csv = cs; hunk ./src/embed.h 1026 - cs[] = 1.; - foreach_face() - fs.x[] = 1.; + csv[] = 1.; + hunk ./src/embed.h 1050 - // fixme: embedded boundaries cannot be combined with (another) metric yet - assert (is_constant (cm) || cm.i == cs.i); + /** + The variables *cms* amd *fms* include both the metric and solid + factors. These variables are intialized to the solid factors + here. The metric factors would be added elsewhere (see for example + [axi.h](axi.h#mg_solve)) */ hunk ./src/embed.h 1056 - cm = cs; - fm = fs; + cms = cs; + fms = fs; + hunk ./src/grid/multigrid-common.h 38 - sum += cm[]*s[]; - s[] = sum/(1 << dimension)/(cm[] + 1e-30); + sum += cms[]*s[]; + s[] = sum/(1 << dimension)/(cms[] + 1e-30); hunk ./src/grid/multigrid-common.h 203 - double sc = s[], cmc = 4.*cm[], sum = cm[]*(1 << dimension); + double sc = s[], cmc = 4.*cms[], sum = cms[]*(1 << dimension); hunk ./src/grid/multigrid-common.h 207 - s[] += child.x*g.x*cm[-child.x]/cmc; - sum -= cm[]; + s[] += child.x*g.x*cms[-child.x]/cmc; + sum -= cms[]; hunk ./src/grid/tree-common.h 170 - if (is_constant(cm)) + if (is_constant(cms)) hunk ./src/grid/tree-common.h 173 - scalar * listr = list_concat ({cm}, p.slist); + scalar * listr = (cms.i == cs.i) ? list_concat ({cms}, p.slist) : + list_concat ({cs, cms}, p.slist); hunk ./src/grid/tree-mpi.h 1251 - scalar * listm = is_constant(cm) ? NULL : (scalar *){fm}; + scalar * listms = NULL; + if (!is_constant (cms)) + listms = (cms.i != cs.i) ? (scalar *) {fs, fms} : (scalar *) {fms}; hunk ./src/grid/tree-mpi.h 1280 - refine_cell (point, listm, 0, NULL); - } + refine_cell (point, listms, 0, NULL); + } hunk ./src/grid/tree-mpi.h 1321 - refine_cell (point, listm, 0, NULL); + refine_cell (point, listms, 0, NULL); hunk ./src/iforce.h 108 - ia.x[] += alpha.x[]/fm.x[]*phif*(f[] - f[-1])/Delta; + ia.x[] += alpha.x[]/(fms.x[] + SEPS)*phif*(f[] - f[-1])/Delta; hunk ./src/input.h 269 - scalar * listm = {cm,fm}; - scalar * listr = !is_constant(cm) ? listm : NULL; + scalar * listm = {cms,fms}; + scalar * listr = !is_constant(cms) ? listm : NULL; hunk ./src/navier-stokes/centered.h 85 -# define neumann_pressure(i) (alpha.n[i] ? a.n[i]*fm.n[i]/alpha.n[i] : \ - a.n[i]*rho[]/(cm[] + SEPS)) +# define neumann_pressure(i) (alpha.n[i] ? a.n[i]*fms.n[i]/alpha.n[i] : \ + a.n[i]*rho[]/(cms[] + SEPS)) hunk ./src/navier-stokes/centered.h 88 -# define neumann_pressure(i) (a.n[i]*fm.n[i]/alpha.n[i]) +# define neumann_pressure(i) (a.n[i]*fms.n[i]/alpha.n[i]) hunk ./src/navier-stokes/centered.h 119 - g->x = rho[]/(cm[] + SEPS)*(a.x[] + a.x[1])/2.; + g->x = rho[]/(cms[] + SEPS)*(a.x[] + a.x[1])/2.; hunk ./src/navier-stokes/centered.h 137 - The default density field is set to unity (times the metric). */ + The default density field is set to unity (times the metric and the solid factors). */ hunk ./src/navier-stokes/centered.h 140 - alpha = fm; - rho = cm; + alpha = fms; + rho = cms; hunk ./src/navier-stokes/centered.h 146 - alphav.x[] = fm.x[]; + alphav.x[] = fms.x[]; hunk ./src/navier-stokes/centered.h 193 - uf.x[] = fm.x[]*face_value (u.x, 0); + uf.x[] = fms.x[]*face_value (u.x, 0); hunk ./src/navier-stokes/centered.h 287 - if (fm.y[i,0] && fm.y[i,1]) { + if (fms.y[i,0] && fms.y[i,1]) { hunk ./src/navier-stokes/centered.h 293 - if (fm.z[i,0,0] && fm.z[i,0,1]) { + if (fms.z[i,0,0] && fms.z[i,0,1]) { hunk ./src/navier-stokes/centered.h 298 - uf.x[] *= fm.x[]; + uf.x[] *= fms.x[]; hunk ./src/navier-stokes/centered.h 384 - uf.x[] = fm.x[]*(face_value (u.x, 0) + dt*a.x[]); + uf.x[] = fms.x[]*(face_value (u.x, 0) + dt*a.x[]); hunk ./src/navier-stokes/centered.h 404 - gf.x[] = fm.x[]*a.x[] - alpha.x[]*(p[] - p[-1])/Delta; + gf.x[] = fms.x[]*a.x[] - alpha.x[]*(p[] - p[-1])/Delta; hunk ./src/navier-stokes/centered.h 414 - g.x[] = (gf.x[] + gf.x[1])/(fm.x[] + fm.x[1] + SEPS); + g.x[] = (gf.x[] + gf.x[1])/(fms.x[] + fms.x[1] + SEPS); hunk ./src/navier-stokes/conserving.h 18 - rhou += cm[]*rho(f[])*u[]; - double du = u[] - rhou/((1 << dimension)*cm[]*rho(f[])); + rhou += cms[]*rho(f[])*u[]; + double du = u[] - rhou/((1 << dimension)*(cms[] + SEPS)*rho(f[])); hunk ./src/navier-stokes/conserving.h 28 - rhou += cm[]*rho(f[])*u[]; - u[] = rhou/((1 << dimension)*cm[]*rho(f[])); + rhou += cms[]*rho(f[])*u[]; + u[] = rhou/((1 << dimension)*(cms[] + SEPS)*rho(f[])); hunk ./src/navier-stokes/double-projection.h 116 - Af.x[] = - fm.x[]*face_value (u.x, 0); + Af.x[] = - fms.x[]*face_value (u.x, 0); hunk ./src/navier-stokes/double-projection.h 134 - Af.x[] += fm.x[]*(face_value (u.x, 0) + dt*a.x[]); + Af.x[] += fms.x[]*(face_value (u.x, 0) + dt*a.x[]); hunk ./src/output.h 1030 - scalar * list = is_constant(cm) ? NULL : list_concat ({cm}, NULL); + scalar * list = NULL; + if (!is_constant (cms)) + list = (cms.i != cs.i) ? list_concat ({cs, cms}, NULL) : list_concat ({cms}, NULL); hunk ./src/output.h 1034 - if (!s.face && !s.nodump && s.i != cm.i) + if (!s.face && !s.nodump && s.i != cms.i && s.i != cs.i) hunk ./src/output.h 1307 + + scalar * listm = NULL; + if (!is_constant (cms)) + listm = (cms.i != cs.i) ? (scalar *) {fs, fms} : (scalar *) {fms}; hunk ./src/output.h 1312 - scalar * listm = is_constant(cm) ? NULL : (scalar *){fm}; hunk ./src/tension.h 39 - We first compute the minimum and maximum values of $\alpha/f_m = + We first compute the minimum and maximum values of $\alpha/(f_{ms}) = hunk ./src/tension.h 43 - foreach_face (reduction(min:amin) reduction(max:amax) reduction(min:dmin)) { - if (alpha.x[]/fm.x[] > amax) amax = alpha.x[]/fm.x[]; - if (alpha.x[]/fm.x[] < amin) amin = alpha.x[]/fm.x[]; - if (Delta < dmin) dmin = Delta; - } + foreach_face (reduction(min:amin) reduction(max:amax) reduction(min:dmin)) + if (fms.x[] > 0.) { + if (alpha.x[]/fms.x[] > amax) amax = alpha.x[]/fms.x[]; + if (alpha.x[]/fms.x[] < amin) amin = alpha.x[]/fms.x[]; + if (Delta < dmin) dmin = Delta; + } hunk ./src/timestep.h 1 -// note: u is weighted by fm +// note: u is weighted by fms hunk ./src/timestep.h 10 - assert (fm.x[]); - dt *= fm.x[]; + assert (fms.x[]); + dt *= fms.x[]; hunk ./src/timestep.h 13 - dt *= cm[]; + dt *= cms[]; hunk ./src/two-phase.h 99 - alphav.x[] = fm.x[]/rho(ff); + alphav.x[] = fms.x[]/rho(ff); hunk ./src/two-phase.h 102 - muv.x[] = fm.x[]*mu(ff); + muv.x[] = fms.x[]*mu(ff); hunk ./src/two-phase.h 106 - rhov[] = cm[]*rho(sf[]); + rhov[] = cms[]*rho(sf[]); hunk ./src/vof.h 61 + + /** + In case of mixed cells, to calculate the concentration value $f_j$ + we have to divide $t_j$ with the effective volume fraction occupied + by the fluid, cs[]. */ + +#if EMBED + cl = cl*cs[-1]; + cc = cc*cs[]; + cr = cr*cs[1]; +#endif + hunk ./src/vof.h 105 - double sc = s.inverse ? s[]/(1. - f[]) : s[]/f[], cmc = 4.*cm[]; + double sc = s.inverse ? s[]/(1. - f[]) : s[]/f[], cmc = 4.*cms[]; hunk ./src/vof.h 109 - s[] += child.x*g.x*cm[-child.x]/cmc; + s[] += child.x*g.x*cms[-child.x]/cmc; hunk ./src/vof.h 201 - double un = uf.x[]*dt/(Delta*fm.x[] + SEPS), s = sign(un); + double un = uf.x[]*dt/(Delta*fms.x[] + SEPS), s = sign(un); hunk ./src/vof.h 207 - if (un*fm.x[]*s/(cm[] + SEPS) > cfl) - cfl = un*fm.x[]*s/(cm[] + SEPS); + if (un*fms.x[]*s/(cms[] + SEPS) > cfl) + cfl = un*fms.x[]*s/(cms[] + SEPS); hunk ./src/vof.h 242 +#if EMBED + ci *= cs[i]; +#endif hunk ./src/vof.h 320 +#if !EMBED hunk ./src/vof.h 322 - c[] += dt*(flux[] - flux[1] + cc[]*(uf.x[1] - uf.x[]))/(cm[]*Delta + SEPS); + c[] += dt*(flux[] - flux[1] + cc[]*(uf.x[1] - uf.x[]))/(cms[]*Delta + SEPS); hunk ./src/vof.h 326 - (cm[]*Delta + SEPS); + (cms[]*Delta + SEPS); hunk ./src/vof.h 328 +#else + vof_update_x (c, cc, flux, tracers, tcl, tfluxl, uf, dt); +#endif [ref files are also updated because fms. Jose M. Lopez-Herrera **20201221152625 Ignore-this: 1bfa6fbdd456661377684baccceaa683 ] hunk ./src/test/axiadvection.ref 54 -0.410526 0.009348046441 0.009348046428 0.0000000009 0.0000000005 -0.421053 0.009348046441 0.009348046425 0.0000000009 0.0000000008 -0.431579 0.009348046441 0.009348046421 0.0000000009 0.0000000012 +0.410526 0.009348046441 0.009348046428 0.0000000008 0.0000000005 +0.421053 0.009348046441 0.009348046425 0.0000000008 0.0000000008 +0.431579 0.009348046441 0.009348046421 0.0000000008 0.0000000012 hunk ./src/test/axiadvection.ref 58 -0.452632 0.009348046446 0.009348046406 0.0000000014 0.0000000029 -0.463158 0.009348046447 0.009348046392 0.0000000015 0.0000000044 -0.473684 0.009348046447 0.009348046372 0.0000000015 0.0000000065 +0.452632 0.009348046445 0.009348046406 0.0000000013 0.0000000029 +0.463158 0.009348046446 0.009348046392 0.0000000014 0.0000000043 +0.473684 0.009348046446 0.009348046372 0.0000000014 0.0000000065 hunk ./src/test/axiadvection.ref 62 -0.494737 0.009348046452 0.009348046301 0.0000000021 0.0000000141 -0.505263 0.009348046463 0.009348046242 0.0000000032 0.0000000204 -0.515789 0.009348046463 0.009348046159 0.0000000032 0.0000000293 -0.526316 0.009348046471 0.009348046044 0.0000000040 0.0000000417 -0.536842 0.009348046471 0.009348045885 0.0000000040 0.0000000586 -0.547368 0.009348046471 0.009348045670 0.0000000040 0.0000000817 -0.557895 0.009348046471 0.009348045379 0.0000000041 0.0000001128 -0.568421 0.009348046471 0.009348044990 0.0000000041 0.0000001544 -0.578947 0.009348046476 0.009348044474 0.0000000046 0.0000002096 -0.589474 0.009348046476 0.009348043793 0.0000000046 0.0000002824 -0.6 0.009348046479 0.009348042902 0.0000000049 0.0000003777 -0.611765 0.009348046479 0.009348041607 0.0000000049 0.0000005163 -0.623529 0.009348046479 0.009348039871 0.0000000049 0.0000007019 -0.635294 0.009348046479 0.009348037560 0.0000000049 0.0000009492 -0.647059 0.009348046484 0.009348034494 0.0000000054 0.0000012772 -0.658824 0.009348046484 0.009348030439 0.0000000054 0.0000017109 -0.670588 0.009348046483 0.009348025081 0.0000000054 0.0000022841 -0.681373 0.009348046482 0.009348018582 0.0000000054 0.0000029793 -0.692157 0.009348046492 0.009348010209 0.0000000063 0.0000038750 -0.702941 0.009348046515 0.009347999423 0.0000000088 0.0000050289 -0.713725 0.009348046529 0.009347985525 0.0000000103 0.0000065155 -0.72451 0.009348046529 0.009347967645 0.0000000103 0.0000084284 -0.735294 0.009348046529 0.009347944701 0.0000000103 0.0000108828 -0.744538 0.009348046608 0.009347919508 0.0000000187 0.0000135778 -0.753782 0.009348046608 0.009347888610 0.0000000187 0.0000168831 -0.763025 0.009348046611 0.009347850753 0.0000000190 0.0000209329 -0.772269 0.009348046664 0.009347804460 0.0000000248 0.0000258853 -0.781513 0.009348046707 0.009347747968 0.0000000293 0.0000319285 -0.790756 0.009348046706 0.009347679160 0.0000000293 0.0000392895 -0.8 0.009348046706 0.009347595371 0.0000000293 0.0000482531 +0.494737 0.009348046451 0.009348046302 0.0000000020 0.0000000140 +0.505263 0.009348046462 0.009348046243 0.0000000031 0.0000000203 +0.515789 0.009348046462 0.009348046160 0.0000000031 0.0000000292 +0.526316 0.009348046469 0.009348046046 0.0000000039 0.0000000414 +0.536842 0.009348046469 0.009348045888 0.0000000039 0.0000000583 +0.547368 0.009348046469 0.009348045673 0.0000000039 0.0000000812 +0.557895 0.009348046469 0.009348045384 0.0000000039 0.0000001122 +0.568421 0.009348046469 0.009348044997 0.0000000039 0.0000001536 +0.578947 0.009348046474 0.009348044483 0.0000000044 0.0000002085 +0.589474 0.009348046474 0.009348043806 0.0000000044 0.0000002810 +0.6 0.009348046477 0.009348042919 0.0000000047 0.0000003759 +0.611765 0.009348046477 0.009348041630 0.0000000047 0.0000005138 +0.623529 0.009348046477 0.009348039903 0.0000000047 0.0000006986 +0.635294 0.009348046477 0.009348037601 0.0000000047 0.0000009448 +0.647059 0.009348046481 0.009348034548 0.0000000052 0.0000012714 +0.658824 0.009348046481 0.009348030510 0.0000000052 0.0000017034 +0.670588 0.009348046481 0.009348025173 0.0000000052 0.0000022743 +0.681373 0.009348046480 0.009348018697 0.0000000052 0.0000029670 +0.692157 0.009348046489 0.009348010354 0.0000000060 0.0000038596 +0.702941 0.009348046510 0.009347999602 0.0000000083 0.0000050097 +0.713725 0.009348046522 0.009347985747 0.0000000096 0.0000064918 +0.72451 0.009348046522 0.009347967918 0.0000000096 0.0000083991 +0.735294 0.009348046522 0.009347945036 0.0000000096 0.0000108469 +0.744538 0.009348046592 0.009347919909 0.0000000170 0.0000135349 +0.753782 0.009348046592 0.009347889091 0.0000000170 0.0000168316 +0.763025 0.009348046594 0.009347851330 0.0000000172 0.0000208712 +0.772269 0.009348046640 0.009347805151 0.0000000222 0.0000258113 +0.781513 0.009348046677 0.009347748794 0.0000000261 0.0000318402 +0.790756 0.009348046677 0.009347680144 0.0000000261 0.0000391842 +0.8 0.009348046677 0.009347596505 0.0000000261 0.0000481318 hunk ./src/test/dirichlet3D.ref 2 -32 3.43e-06 0.000317 7.99e-06 0.000317 2.53e-06 2.12e-05 36 31 +32 6.09e-06 0.000317 1.43e-05 0.000317 2.53e-06 2.12e-05 36 31 hunk ./src/test/dirichlet3D.ref 4 -64 4.13e-07 7.55e-05 1.07e-06 7.55e-05 3.52e-07 3.38e-06 33 23 +64 6.47e-07 7.55e-05 2.14e-06 7.55e-05 3.52e-07 3.38e-06 33 23 hunk ./src/test/dirichlet3D.ref 6 -128 5.63e-08 5.34e-06 1.47e-07 5.34e-06 5.21e-08 4.68e-07 22 4 +128 7.28e-08 5.34e-06 2.94e-07 5.34e-06 5.21e-08 4.68e-07 22 4 hunk ./src/test/dirichlet3D.ref 8 -256 8.29e-09 1.54e-06 2.12e-08 1.54e-06 8e-09 6.16e-08 20 4 +256 1.05e-08 1.54e-06 6.43e-08 1.54e-06 8e-09 6.16e-08 20 4 hunk ./src/test/fall.ref 6 -0.181654 0.999922 +0.181654 0.999924 hunk ./src/test/fall.ref 8 -0.261653 1.00033 +0.261653 1.00034 hunk ./src/test/fall.ref 11 -0.381653 1.00098 -0.421653 1.00118 -0.461653 1.00136 -0.501653 1.00151 -0.541653 1.00161 -0.581653 1.00167 -0.621653 1.00188 -0.661653 1.00211 -0.701653 1.00228 -0.741653 1.00238 -0.781653 1.00239 -0.821653 1.00251 -0.861653 1.00275 -0.901653 1.00286 -0.941653 1.00285 -0.981653 1.00302 -1.02165 1.00324 -1.06165 1.00329 -1.10165 1.00339 -1.14165 1.00367 -1.18165 1.00373 -1.22165 1.00397 -1.26165 1.00425 -1.30165 1.00431 -1.34165 1.00486 -1.37254 1.00537 -1.40767 1.00741 -1.44688 1.01258 -1.48676 1.02356 -1.52674 1.08482 -1.56674 1.19664 -1.60674 1.29492 -1.64674 1.37725 -1.68674 1.44427 -1.72674 1.49826 -1.76674 1.54668 -1.80674 1.58644 -1.84674 1.62194 -1.88674 1.65071 -1.92674 1.67345 -1.96674 1.68978 -2.00674 1.70309 -2.04674 1.71178 -2.08674 1.71656 -2.12674 1.71786 -2.16674 1.71618 -2.20674 1.71198 -2.24674 1.70564 -2.28674 1.69678 -2.32674 1.68681 -2.36674 1.67609 -2.40674 1.66479 -2.44674 1.65229 -2.48674 1.63983 -2.52674 1.62739 -2.56674 1.6145 -2.60674 1.60173 -2.64674 1.5893 -2.68674 1.57683 -2.72674 1.56464 -2.76674 1.55306 -2.80674 1.5419 -2.84674 1.53097 -2.88674 1.52067 -2.92674 1.51104 -2.96674 1.50195 -3.00674 1.4936 -3.04674 1.4865 -3.08674 1.48027 -3.12674 1.47486 -3.16674 1.47027 -3.20674 1.46649 -3.24674 1.46349 -3.28674 1.46115 -3.32674 1.45953 -3.36674 1.45877 -3.40674 1.45876 -3.44674 1.45968 -3.48674 1.46158 -3.52674 1.46406 -3.56674 1.46696 -3.60674 1.47027 -3.64674 1.4742 -3.68674 1.47865 -3.72674 1.48342 -3.76674 1.48845 -3.80674 1.49371 -3.84674 1.49912 -3.88674 1.50449 -3.92674 1.50964 -3.96674 1.51475 -4.00674 1.51975 -4.04674 1.52458 -4.08674 1.52919 -4.12674 1.53358 -4.16674 1.53772 -4.20674 1.54194 -4.24674 1.54568 -4.28674 1.54889 -4.32674 1.55169 -4.36674 1.55407 -4.40674 1.55606 -4.44674 1.55769 -4.48674 1.55898 -4.52674 1.55996 -4.56674 1.56067 -4.60674 1.56114 -4.64674 1.5614 -4.68674 1.56148 -4.72674 1.56143 -4.76674 1.56126 -4.80674 1.561 -4.84674 1.56069 -4.88674 1.56036 -4.92674 1.56001 -4.96674 1.55968 +0.381653 1.00099 +0.421653 1.00119 +0.461653 1.00137 +0.501653 1.00152 +0.541653 1.00162 +0.581653 1.00168 +0.621653 1.00189 +0.661653 1.00213 +0.701653 1.0023 +0.741653 1.00239 +0.781653 1.0024 +0.821653 1.00252 +0.861653 1.00276 +0.901653 1.00288 +0.941653 1.00286 +0.981653 1.00303 +1.02165 1.00325 +1.06165 1.0033 +1.10165 1.0034 +1.14165 1.00368 +1.18165 1.00374 +1.22165 1.00398 +1.26165 1.00426 +1.30165 1.00433 +1.34165 1.00488 +1.37256 1.00539 +1.40765 1.00743 +1.44685 1.01259 +1.48673 1.02357 +1.52672 1.08478 +1.56671 1.19662 +1.60671 1.29491 +1.64671 1.37726 +1.68671 1.44429 +1.72671 1.49828 +1.76671 1.54671 +1.80671 1.58646 +1.84671 1.62198 +1.88671 1.65074 +1.92671 1.67349 +1.96671 1.68981 +2.00671 1.70312 +2.04671 1.7118 +2.08671 1.71657 +2.12671 1.71787 +2.16671 1.71618 +2.20671 1.71197 +2.24671 1.70562 +2.28671 1.69674 +2.32671 1.68676 +2.36671 1.67602 +2.40671 1.66472 +2.44671 1.65221 +2.48671 1.63975 +2.52671 1.6273 +2.56671 1.61439 +2.60671 1.60163 +2.64671 1.5892 +2.68671 1.57673 +2.72671 1.56454 +2.76671 1.55296 +2.80671 1.5418 +2.84671 1.53088 +2.88671 1.52059 +2.92671 1.51098 +2.96671 1.50189 +3.00671 1.49359 +3.04671 1.48651 +3.08671 1.4803 +3.12671 1.47492 +3.16671 1.47035 +3.20671 1.4666 +3.24671 1.46363 +3.28671 1.46134 +3.32671 1.45974 +3.36671 1.45902 +3.40671 1.45906 +3.44671 1.45999 +3.48671 1.46185 +3.52671 1.46433 +3.56671 1.46723 +3.60671 1.47072 +3.64671 1.47475 +3.68671 1.47919 +3.72671 1.48397 +3.76671 1.48902 +3.80671 1.4943 +3.84671 1.49973 +3.88671 1.50512 +3.92671 1.51029 +3.96671 1.51543 +4.00671 1.52045 +4.04671 1.52532 +4.08671 1.52996 +4.12671 1.53439 +4.16671 1.53855 +4.20671 1.54286 +4.24671 1.54655 +4.28671 1.54976 +4.32671 1.55257 +4.36671 1.55496 +4.40671 1.55697 +4.44671 1.55861 +4.48671 1.5599 +4.52671 1.56088 +4.56671 1.56159 +4.60671 1.56206 +4.64671 1.5623 +4.68671 1.56237 +4.72671 1.56229 +4.76671 1.56211 +4.80671 1.56184 +4.84671 1.56149 +4.88671 1.56111 +4.92671 1.56072 +4.96671 1.56035 hunk ./src/test/neumann3D.ref 2 -32 5.3e-06 0.000255 9.91e-06 0.000255 4.39e-06 2.07e-05 35 31 +32 7.69e-06 0.000255 1.53e-05 0.000255 4.39e-06 2.07e-05 35 31 hunk ./src/test/neumann3D.ref 4 -64 1.21e-06 6.94e-05 2.31e-06 6.94e-05 1.11e-06 5.77e-06 32 21 +64 1.43e-06 6.94e-05 3.04e-06 6.94e-05 1.11e-06 5.77e-06 32 21 hunk ./src/test/neumann3D.ref 6 -128 2.92e-07 5.02e-06 5.69e-07 5.02e-06 2.79e-07 1.47e-06 20 4 +128 3.08e-07 5.02e-06 6.18e-07 5.02e-06 2.79e-07 1.47e-06 20 4 hunk ./src/test/neumann3D.ref 8 -256 7.14e-08 1.48e-06 1.41e-07 1.48e-06 6.99e-08 3.72e-07 19 4 +256 7.34e-08 1.48e-06 1.5e-07 1.48e-06 6.99e-08 3.72e-07 19 4 hunk ./src/test/oscillation-momentum.ref 1 -fit 6.4 0.000300 1.97 153 +fit 6.4 0.000298 1.92 153 hunk ./src/test/oscillation-momentum.ref 4 -fit 51.2 0.000275 0.02 154 +fit 51.2 0.000294 0.02 154 hunk ./src/test/poiseuille-axi.ref 1 -8 0.000976587 0.000976587 0.000976651 +8 0.000976587 0.000976587 0.000976652 hunk ./src/test/rising-axi.ref 1 -1.32386 0 -1.3239 0.0078125 - -1.3239 0.0078125 -1.324 0.015625 - -1.324 0.015625 -1.32417 0.0234375 - -1.32417 0.0234375 -1.3244 0.03125 - -1.3244 0.03125 -1.32471 0.0390625 - -1.32471 0.0390625 -1.32509 0.046875 - -1.32509 0.046875 -1.32554 0.0546875 - -1.32554 0.0546875 -1.32607 0.0625 - -1.32607 0.0625 -1.32667 0.0703125 - -1.32667 0.0703125 -1.32735 0.078125 - -1.32735 0.078125 -1.32812 0.0859375 - -1.32812 0.0859375 -1.32812 0.0859887 - -1.32812 0.0859755 -1.32898 0.09375 - -1.32898 0.09375 -1.32993 0.101562 - -1.32993 0.101562 -1.33097 0.109375 - -1.33097 0.109375 -1.33212 0.117188 - -1.33212 0.117188 -1.33338 0.125 - -1.33338 0.125 -1.33475 0.132812 - -1.33475 0.132812 -1.33594 0.139014 - -1.33594 0.138997 -1.33625 0.140625 - -1.33624 0.140625 -1.33788 0.148438 - -1.33788 0.148438 -1.33965 0.15625 - -1.33965 0.15625 -1.34159 0.164062 - -1.34158 0.164062 -1.34369 0.171875 - -1.34367 0.171875 -1.34375 0.172154 - -1.34375 0.172102 -1.34598 0.179688 - -1.34597 0.179688 -1.34847 0.1875 - -1.34846 0.1875 -1.35119 0.195312 - -1.35116 0.195312 -1.35156 0.196364 - -1.35156 0.196295 -1.35417 0.203125 - -1.35416 0.203125 -1.35743 0.210938 - -1.35741 0.210938 -1.35938 0.215208 - -1.35938 0.215159 -1.36103 0.21875 - -1.36101 0.21875 -1.36498 0.226562 - -1.36494 0.226562 -1.36719 0.230532 - -1.36719 0.230462 -1.3694 0.234375 - -1.36936 0.234375 -1.37431 0.242188 - -1.37424 0.242188 -1.375 0.243308 - -1.375 0.243184 -1.37988 0.25 - -1.37986 0.25 -1.38281 0.253777 - -1.38281 0.253628 -1.38623 0.257812 - -1.38623 0.257812 -1.39062 0.262616 - -1.39062 0.262561 -1.39368 0.265625 - -1.3936 0.265625 -1.39844 0.270057 - -1.39844 0.270061 -1.40267 0.273438 - -1.40247 0.273438 -1.40625 0.276314 - -1.40625 0.276277 -1.41373 0.28125 - -1.41349 0.28125 -1.41406 0.281601 - -1.41406 0.281478 -1.42188 0.285731 - -1.42188 0.285677 -1.42962 0.289062 - -1.42939 0.289062 -1.42969 0.289194 - -1.42969 0.289111 -1.4375 0.291764 - -1.4375 0.291731 -1.44531 0.293659 - -1.44531 0.293648 -1.45312 0.294942 - -1.45312 0.294936 -1.46094 0.295633 - -1.46094 0.295628 -1.46875 0.295747 - -1.46875 0.295743 -1.47656 0.295304 - -1.47656 0.2953 -1.48438 0.29432 - -1.48438 0.294315 -1.49219 0.292815 - -1.49219 0.292805 -1.5 0.290814 - -1.625 0.185492 -1.62378 0.1875 - -1.57812 0.243755 -1.57124 0.25 - -1.59375 0.227621 -1.58754 0.234375 - -1.58594 0.236034 -1.57974 0.242188 - -1.57812 0.243704 -1.57971 0.242188 - -1.58594 0.235988 -1.58749 0.234375 - -1.60156 0.218618 -1.60146 0.21875 - -1.60938 0.208496 -1.60758 0.210938 - -1.60156 0.218454 -1.60752 0.210938 - -1.61719 0.19751 -1.6133 0.203125 - -1.62374 0.1875 -1.61869 0.195312 - -1.61719 0.197481 -1.61865 0.195312 - -1.60938 0.208388 -1.61327 0.203125 - -1.60135 0.21875 -1.5947 0.226562 - -1.59375 0.227601 -1.59467 0.226562 +1.32404 0 +1.32408 0.0078125 + +1.32408 0.0078125 +1.32418 0.015625 + +1.32418 0.015625 +1.32435 0.0234375 + +1.32435 0.0234375 +1.32458 0.03125 + +1.32458 0.03125 +1.32489 0.0390625 + +1.32489 0.0390625 +1.32527 0.046875 + +1.32527 0.046875 +1.32572 0.0546875 + +1.32572 0.0546875 +1.32624 0.0625 + +1.32624 0.0625 +1.32685 0.0703125 + +1.32685 0.0703125 +1.32753 0.078125 + +1.32753 0.078125 +1.32812 0.0841714 + +1.32812 0.0841621 +1.3283 0.0859375 + +1.3283 0.0859375 +1.32915 0.09375 + +1.32915 0.09375 +1.3301 0.101562 + +1.3301 0.101562 +1.33115 0.109375 + +1.33115 0.109375 +1.3323 0.117188 + +1.3323 0.117188 +1.33355 0.125 + +1.33355 0.125 +1.33492 0.132812 + +1.33492 0.132812 +1.33594 0.138119 + +1.33594 0.138101 +1.33642 0.140625 + +1.33642 0.140625 +1.33805 0.148438 + +1.33805 0.148438 +1.33982 0.15625 + +1.33982 0.15625 +1.34175 0.164062 + +1.34175 0.164062 +1.34375 0.171491 + +1.34375 0.171478 +1.34386 0.171875 + +1.34385 0.171875 +1.34614 0.179688 + +1.34614 0.179688 +1.34864 0.1875 + +1.34863 0.1875 +1.35135 0.195312 + +1.35132 0.195312 +1.35156 0.195942 + +1.35156 0.19587 +1.35433 0.203125 + +1.35432 0.203125 +1.35759 0.210938 + +1.35756 0.210938 +1.35938 0.214868 + +1.35938 0.214813 +1.36119 0.21875 + +1.36116 0.21875 +1.36514 0.226562 + +1.3651 0.226562 +1.36719 0.230263 + +1.36719 0.230189 +1.36955 0.234375 + +1.36951 0.234375 +1.37446 0.242188 + +1.37439 0.242188 +1.375 0.243097 + +1.375 0.242977 +1.38002 0.25 + +1.38 0.25 +1.38281 0.253602 + +1.38281 0.253455 +1.38637 0.257812 + +1.38637 0.257812 +1.39062 0.262474 + +1.39062 0.262418 +1.39382 0.265625 + +1.39373 0.265625 +1.39844 0.269948 + +1.39844 0.26995 +1.40279 0.273438 + +1.40258 0.273438 +1.40625 0.276229 + +1.40625 0.276194 +1.41382 0.28125 + +1.41358 0.28125 +1.41406 0.281547 + +1.41406 0.28142 +1.42188 0.285696 + +1.42188 0.285641 +1.42965 0.289062 + +1.42942 0.289062 +1.42969 0.28918 + +1.42969 0.289095 +1.4375 0.291768 + +1.4375 0.291734 +1.44531 0.29368 + +1.44531 0.293668 +1.45312 0.29498 + +1.45312 0.294974 +1.46094 0.295686 + +1.46094 0.295681 +1.46875 0.295816 + +1.46875 0.295811 +1.47656 0.295386 + +1.47656 0.295382 +1.48438 0.294416 + +1.48438 0.29441 +1.49219 0.292921 + +1.49219 0.292911 +1.5 0.290929 + +1.625 0.185628 +1.62387 0.1875 + +1.57812 0.24391 +1.57141 0.25 + +1.59375 0.227767 +1.58768 0.234375 + +1.58594 0.236181 +1.5799 0.242188 + +1.57812 0.243858 +1.57987 0.242188 + +1.58594 0.236138 +1.58763 0.234375 + +1.60938 0.208644 +1.60769 0.210938 + +1.60156 0.218607 +1.60763 0.210938 + +1.61719 0.197644 +1.6134 0.203125 + +1.62382 0.1875 +1.61878 0.195312 + +1.61719 0.197613 +1.61874 0.195312 + +1.60938 0.20853 +1.61337 0.203125 + +1.60147 0.21875 +1.59483 0.226562 + +1.59375 0.227745 +1.5948 0.226562 hunk ./src/test/rising-axi.ref 224 -1.67592 0.0078125 - -1.67592 0.0078125 -1.67567 0.015625 - -1.67567 0.015625 -1.67526 0.0234375 - -1.67526 0.0234375 +1.67593 0.0078125 + +1.67593 0.0078125 +1.67568 0.015625 + +1.67568 0.015625 +1.67527 0.0234375 + +1.67527 0.0234375 hunk ./src/test/rising-axi.ref 235 -1.67188 0.0552192 -1.67071 0.0625 +1.67188 0.0552955 +1.67072 0.0625 hunk ./src/test/rising-axi.ref 239 -1.67394 0.0390625 - -1.67394 0.0390625 -1.67304 0.046875 - -1.67304 0.046875 -1.67196 0.0546875 - -1.67188 0.0551845 -1.67195 0.0546875 - -1.65625 0.119025 -1.65415 0.125 - -1.67071 0.0625 -1.66928 0.0703125 - -1.66928 0.0703125 -1.66767 0.078125 - -1.66406 0.0933002 -1.66395 0.09375 - -1.66767 0.078125 -1.66589 0.0859375 - -1.66406 0.0932153 -1.66589 0.0859375 - -1.66394 0.09375 -1.66178 0.101562 - -1.66178 0.101562 -1.65943 0.109375 - -1.65943 0.109375 -1.65689 0.117188 - -1.65625 0.118985 -1.65688 0.117188 - -1.64844 0.139642 -1.64804 0.140625 - -1.65414 0.125 -1.65118 0.132812 - -1.64844 0.139537 -1.65117 0.132812 - -1.64801 0.140625 -1.6446 0.148438 - -1.64459 0.148438 -1.64095 0.15625 - -1.64062 0.156925 -1.63706 0.164062 - -1.63706 0.164062 -1.63289 0.171875 - -1.63281 0.172047 -1.62847 0.179688 - -1.625 0.185375 -1.62845 0.179688 - -1.63281 0.172022 -1.63289 0.171875 - -1.64062 0.156898 -1.64095 0.15625 - -1.53125 0.278029 -1.52484 0.28125 - -1.50781 0.288355 -1.50557 0.289062 - -1.5 0.290816 -1.50555 0.289062 - -1.50781 0.288359 -1.51562 0.285387 - -1.51562 0.285361 -1.52344 0.28192 - -1.52344 0.281945 -1.52482 0.28125 - -1.53125 0.277995 -1.53906 0.273579 - -1.54688 0.268711 -1.53933 0.273438 - -1.53906 0.273604 -1.53934 0.273438 - -1.55469 0.263338 -1.5514 0.265625 - -1.5625 0.257473 -1.56207 0.257812 - -1.55469 0.263265 -1.56184 0.257812 - -1.54688 0.268669 -1.55126 0.265625 - -1.5625 0.257329 -1.57031 0.250828 - -1.57031 0.250801 -1.57123 0.25 +1.67395 0.0390625 + +1.67395 0.0390625 +1.67305 0.046875 + +1.67305 0.046875 +1.67197 0.0546875 + +1.67188 0.0552601 +1.67197 0.0546875 + +1.65625 0.119124 +1.65419 0.125 + +1.67072 0.0625 +1.66929 0.0703125 + +1.66929 0.0703125 +1.66769 0.078125 + +1.66406 0.0933954 +1.66397 0.09375 + +1.66769 0.078125 +1.66591 0.0859375 + +1.66406 0.0933092 +1.66591 0.0859375 + +1.66396 0.09375 +1.66181 0.101562 + +1.6618 0.101562 +1.65946 0.109375 + +1.65946 0.109375 +1.65692 0.117188 + +1.65625 0.119083 +1.65692 0.117188 + +1.64844 0.139757 +1.64809 0.140625 + +1.65418 0.125 +1.65122 0.132812 + +1.64844 0.139651 +1.65121 0.132812 + +1.64806 0.140625 +1.64465 0.148438 + +1.64465 0.148438 +1.641 0.15625 + +1.64062 0.15704 +1.63713 0.164062 + +1.63712 0.164062 +1.63296 0.171875 + +1.63281 0.172166 +1.62855 0.179688 + +1.625 0.185508 +1.62853 0.179688 + +1.63281 0.172148 +1.63296 0.171875 + +1.64062 0.157011 +1.641 0.15625 + +1.53125 0.278172 +1.52512 0.28125 + +1.50781 0.288478 +1.50596 0.289062 + +1.5 0.290929 +1.50593 0.289062 + +1.50781 0.288483 +1.51562 0.28552 + +1.51562 0.285493 +1.52344 0.282058 + +1.52344 0.282081 +1.52509 0.28125 + +1.53125 0.278136 +1.53906 0.273726 + +1.54688 0.268864 +1.53957 0.273438 + +1.53906 0.273749 +1.53958 0.273438 + +1.55469 0.263496 +1.55162 0.265625 + +1.5625 0.257639 +1.56228 0.257812 + +1.55469 0.263419 +1.56204 0.257812 + +1.54688 0.268821 +1.55147 0.265625 + +1.5625 0.257482 +1.57031 0.250983 + +1.57031 0.25095 +1.5714 0.25 [Since now fm is fms, update is needed for several solvers. Jose M. Lopez-Herrera **20201221154155 Ignore-this: 3f1e37d561b4e28a69e8a80eaaa27849 In addition some metric factors lacking in all-mach & momentum.h are completed. ] hunk ./src/all-mach.h 71 - alpha = fm; + alpha = fms; + + if (rho.i == unity.i) + rho = cms; hunk ./src/all-mach.h 148 + hunk ./src/all-mach.h 162 - q.x[] = (q.x[] + dt*g.x[])/rho[]; + q.x[] = (q.x[] + dt*g.x[])*cms[]/rho[]; hunk ./src/all-mach.h 167 - q.x[] = q.x[]*rho[] - dt*g.x[]; + q.x[] = q.x[]*rho[]/cms[] - dt*g.x[]; hunk ./src/all-mach.h 169 - } + } hunk ./src/all-mach.h 177 - uf.x[] = alpha.x[]*(q.x[] + q.x[-1])/2. + dt*fm.x[]*a.x[]; + uf.x[] = alpha.x[]*(q.x[] + q.x[-1])/2. + dt*fms.x[]*a.x[]; hunk ./src/all-mach.h 221 - lambda[] = - cm[]/(sq(dt)*rhoc2[]); + lambda[] = - cms[]/(sq(dt)*rhoc2[]); hunk ./src/all-mach.h 237 - tolerance](poisson.h#project) identical to that used for - incompressible flows. */ + tolerance](poisson.h#377) identical to that used for incompressible + flows. */ hunk ./src/all-mach.h 253 - gf.x[] = a.x[] - dp/fm.x[]; + gf.x[] = a.x[] - dp/fms.x[]; hunk ./src/all-mach.h 264 - g.x[] = rho[]*(gf.x[] + gf.x[1])/2.; + g.x[] = rho[]*(gf.x[] + gf.x[1])/(2.*cms[]); hunk ./src/all-mach.h 270 -/** -Some derived solvers need to hook themselves at the end of the -timestep. */ - -event end_timestep (i++, last); - hunk ./src/axi.h 20 -static void refine_cm_axi (Point point, scalar cm) +static void refine_cm_axi (Point point, scalar cms) hunk ./src/axi.h 22 - fine(cm,0,0) = fine(cm,1,0) = y - Delta/4.; - fine(cm,0,1) = fine(cm,1,1) = y + Delta/4.; + fine(cms,0,0) = fine(cms,1,0) = y - Delta/4.; + fine(cms,0,1) = fine(cms,1,1) = y + Delta/4.; hunk ./src/axi.h 26 -static void refine_face_x_axi (Point point, scalar fm) +static void refine_face_x_axi (Point point, scalar fms) hunk ./src/axi.h 29 - fine(fm,0,0) = y - Delta/4.; - fine(fm,0,1) = y + Delta/4.; + fine(fms,0,0) = y - Delta/4.; + fine(fms,0,1) = y + Delta/4.; hunk ./src/axi.h 33 - fine(fm,2,0) = y - Delta/4.; - fine(fm,2,1) = y + Delta/4.; + fine(fms,2,0) = y - Delta/4.; + fine(fms,2,1) = y + Delta/4.; hunk ./src/axi.h 36 - fine(fm,1,0) = y - Delta/4.; - fine(fm,1,1) = y + Delta/4.; + fine(fms,1,0) = y - Delta/4.; + fine(fms,1,1) = y + Delta/4.; hunk ./src/axi.h 40 -static void refine_face_y_axi (Point point, scalar fm) +static void refine_face_y_axi (Point point, scalar fms) hunk ./src/axi.h 43 - fine(fm,0,0) = fine(fm,1,0) = max(y - Delta/2., 1e-20); + fine(fms,0,0) = fine(fms,1,0) = max(y - Delta/2., 1e-20); hunk ./src/axi.h 45 - fine(fm,0,2) = fine(fm,1,2) = y + Delta/2.; - fine(fm,0,1) = fine(fm,1,1) = y; + fine(fms,0,2) = fine(fms,1,2) = y + Delta/2.; + fine(fms,0,1) = fine(fms,1,1) = y; hunk ./src/axi.h 53 - By default *cm* is a constant scalar field. To make it variable, we + By default *cms* is a constant scalar field. To make it variable, we hunk ./src/axi.h 58 - if (is_constant(cm)) { + if (is_constant (cms)) { hunk ./src/axi.h 60 - cm = new scalar; + cms = new scalar; hunk ./src/axi.h 62 - all = list_concat ({cm}, l); + all = list_concat ({cms}, l); hunk ./src/axi.h 67 + If embeded solid are presents we have to allocate a new field for + the combined metric and solid factors, *cms* beside to the already + defined *cs*. Also the first in the list of the variables should be + those corresponding to solid factors. This will be done a init event */ + + +#if EMBED + if (cms.i == cs.i) + cms = new scalar; +#endif + + /** hunk ./src/axi.h 83 - scalar cmv = cm; + scalar cmsv = cms; hunk ./src/axi.h 85 - cmv[] = y; - cm[top] = dirichlet(y); - cm[bottom] = dirichlet(y); + cmsv[] = y; + cms[top] = dirichlet(y); + cms[bottom] = dirichlet(y); hunk ./src/axi.h 95 - if (is_constant(fm.x)) { + if (is_constant(fms.x)) { hunk ./src/axi.h 97 - fm = new face vector; + fms = new face vector; hunk ./src/axi.h 99 - all = list_concat ((scalar *){fm}, l); + all = list_concat ((scalar *){fms}, l); hunk ./src/axi.h 102 - face vector fmv = fm; + + #if EMBED + if (fms.x.i == fs.x.i) + fms = new face vector; +#endif + + + face vector fmsv = fms; hunk ./src/axi.h 111 - fmv.x[] = max(y, 1e-20); - fm.t[top] = dirichlet(y); - fm.t[bottom] = dirichlet(y); + fmsv.x[] = max(y, 1e-20); + fms.t[top] = dirichlet(y); + fms.t[bottom] = dirichlet(y); hunk ./src/axi.h 119 - cm.refine = cm.prolongation = refine_cm_axi; - fm.x.prolongation = refine_face_x_axi; - fm.y.prolongation = refine_face_y_axi; + cms.refine = cms.prolongation = refine_cm_axi; + fms.x.prolongation = refine_face_x_axi; + fms.y.prolongation = refine_face_y_axi; hunk ./src/axi.h 124 - boundary ({cm, fm}); + boundary ({cms, fms}); hunk ./src/conservation.h 233 - dtmax = riemann (r, l, Delta*cm[]/fm.x[], f, len, dtmax); + dtmax = riemann (r, l, Delta*cms[]/fms.x[], f, len, dtmax); hunk ./src/conservation.h 236 - fs.x[] = fm.x[]*f[i++]; + fs.x[] = fms.x[]*f[i++]; hunk ./src/conservation.h 238 - fv.x.x[] = fm.x[]*f[i++]; + fv.x.x[] = fms.x[]*f[i++]; hunk ./src/conservation.h 240 - fv.y.x[] = fm.x[]*f[i++]; + fv.y.x[] = fms.x[]*f[i++]; hunk ./src/conservation.h 243 - fv.z.x[] = fm.x[]*f[i++]; + fv.z.x[] = fms.x[]*f[i++]; hunk ./src/conservation.h 261 - ds[] += (f.x[] - f.x[1])/(cm[]*Delta); + ds[] += (f.x[] - f.x[1])/(cms[]*Delta); hunk ./src/contact.h 12 + +#undef interface_normal hunk ./src/ehd/implicit.h 52 - epsilon.x[] = K.x[] = fm.x[]; + epsilon.x[] = K.x[] = fms.x[]; hunk ./src/ehd/implicit.h 64 - rhs[] = - rhoe[]*cm[]; + rhs[] = - rhoe[]*cms[]; hunk ./src/ehd/implicit.h 93 - rhoe[] /= cm[]*Delta; + rhoe[] /= cms[]*Delta; hunk ./src/ehd/implicit.h 101 - rhoe[] /= cm[]*sq(Delta); + rhoe[] /= cms[]*sq(Delta); hunk ./src/ehd/stress.h 49 - f.x[] = (Mx.x[1,0] - Mx.x[] + Mx.y[0,1] - Mx.y[])/(Delta*cm[]); + f.x[] = (Mx.x[1,0] - Mx.x[] + Mx.y[0,1] - Mx.y[])/(Delta*cms[]); hunk ./src/ehd/stress.h 59 - sq(phi[0,1] - phi[0,-1]))/(8.*cm[]*sq(Delta)) - *(epsilon.x[]/fm.x[] + epsilon.y[]/fm.y[] + - epsilon.x[1,0]/fm.x[1,0] + epsilon.y[0,1]/fm.y[0,1])/4.; + sq(phi[0,1] - phi[0,-1]))/(8.*cms[]*sq(Delta)) + *(epsilon.x[]/fms.x[] + epsilon.y[]/fms.y[] + + epsilon.x[1,0]/fms.x[1,0] + epsilon.y[0,1]/fms.y[0,1])/4.; hunk ./src/ehd/stress.h 71 - av.x[] += alpha.x[]/fm.x[]*(f.x[] + f.x[-1])/2.; + av.x[] += alpha.x[]/fms.x[]*(f.x[] + f.x[-1])/2.; hunk ./src/layered/hydro.h 164 - double ax = - fm.x[]*G*(eta[] - eta[-1])/((cm[] + cm[-1])*Delta/2.); + double ax = - fms.x[]*G*(eta[] - eta[-1])/((cms[] + cms[-1])*Delta/2.); hunk ./src/layered/hydro.h 196 - hf.x[] *= fm.x[]; + hf.x[] *= fms.x[]; hunk ./src/layered/hydro.h 256 - Hf /= fm.x[]; + Hf /= fms.x[]; hunk ./src/layered/hydro.h 262 - double dt = hf.x[]/fm.x[]*cm[]*Delta/c; + double dt = hf.x[]/fms.x[]*cms[]*Delta/c; hunk ./src/layered/hydro.h 328 - s[] += dt*(flux.x[] - flux.x[1])/(Delta*cm[]); + s[] += dt*(flux.x[] - flux.x[1])/(Delta*cms[]); hunk ./src/layered/hydro.h 344 - h[] += dt*(hu.x[] - hu.x[1])/(Delta*cm[]); + h[] += dt*(hu.x[] - hu.x[1])/(Delta*cms[]); hunk ./src/layered/hydro.h 411 - double dmdl = (fm.x[1,0] - fm.x[])/(cm[]*Delta); - double dmdt = (fm.y[0,1] - fm.y[])/(cm[]*Delta); + double dmdl = (fms.x[1,0] - fms.x[])/(cms[]*Delta); + double dmdt = (fms.y[0,1] - fms.y[])/(cms[]*Delta); hunk ./src/layered/hydro.h 585 - double dp = (p[1].x - p[0].x)*Delta/Delta_x*(fm.y[] + fm.y[0,1])/2.; + double dp = (p[1].x - p[0].x)*Delta/Delta_x*(fms.y[] + fms.y[0,1])/2.; hunk ./src/layered/implicit.h 135 - alpha.x[] *= 2.*fm.x[]/(cm[] + cm[-1]); + alpha.x[] *= 2.*fms.x[]/(cms[] + cms[-1]); hunk ./src/layered/implicit.h 188 - double ax = - fm.x[]*G*(etap[] - etap[-1])/((cm[] + cm[-1])*Delta/2.); + double ax = - fms.x[]*G*(etap[] - etap[-1])/((cms[] + cms[-1])*Delta/2.); hunk ./src/layered/nh-staggered.h 680 - double ax = - fm.x[]*(ac[l]*phib[-1] + bc[l][0]*phi[-1] + cc[l]*phit[-1] + + double ax = - fms.x[]*(ac[l]*phib[-1] + bc[l][0]*phi[-1] + cc[l]*phit[-1] + hunk ./src/layered/nh.h 251 - rhs[] += ((h[] + h[1])*fm.x[1]* + rhs[] += ((h[] + h[1])*fms.x[1]* hunk ./src/layered/nh.h 253 - (h[] + h[-1])*fm.x[]* + (h[] + h[-1])*fms.x[]* hunk ./src/layered/nh.h 298 - ax = fm.x[]*((h[]*phi[] - h[-1]*phi[-1])/Delta + + ax = fms.x[]*((h[]*phi[] - h[-1]*phi[-1])/Delta + hunk ./src/layered/nh.h 302 - ax = fm.x[]* + ax = fms.x[]* hunk ./src/layered/nh.h 313 - ax = fm.x[]*((h[]*phi[] - h[-1]*phi[-1])/Delta + + ax = fms.x[]*((h[]*phi[] - h[-1]*phi[-1])/Delta + hunk ./src/layered/nh.h 316 - ax = fm.x[]*((h[]*(phi[] + phi[0,0,1]) - + ax = fms.x[]*((h[]*(phi[] + phi[0,0,1]) - hunk ./src/log-conform.h 541 - if (fm.x[] > 1e-20) { - double shear = (tau_p.x.y[0,1]*cm[0,1] + tau_p.x.y[-1,1]*cm[-1,1] - - tau_p.x.y[0,-1]*cm[0,-1] - tau_p.x.y[-1,-1]*cm[-1,-1])/4.; - av.x[] += (shear + cm[]*tau_p.x.x[] - cm[-1]*tau_p.x.x[-1])* - alpha.x[]/(sq(fm.x[])*Delta); + if (fms.x[] > 1e-20) { + double shear = (tau_p.x.y[0,1]*cms[0,1] + tau_p.x.y[-1,1]*cms[-1,1] - + tau_p.x.y[0,-1]*cms[0,-1] - tau_p.x.y[-1,-1]*cms[-1,-1])/4.; + av.x[] += (shear + cms[]*tau_p.x.x[] - cms[-1]*tau_p.x.x[-1])* + alpha.x[]/(sq(fms.x[])*Delta); hunk ./src/momentum.h 44 - rhov[] = rho(f[]); + rhov[] = rho(f[])*cms[]; hunk ./src/momentum.h 47 - alphav.x[] = 2./(rhov[] + rhov[-1]); hunk ./src/momentum.h 48 - muv.x[] = fm.x[]*mu(ff); + alphav.x[] = fms.x[]/rho(ff); + muv.x[] = fms.x[]*mu(ff); hunk ./src/radial.h 30 - if (is_constant(cm)) { + if (is_constant(cms)) { hunk ./src/radial.h 32 - cm = new scalar; + cms = new scalar; hunk ./src/radial.h 34 - all = list_concat ({cm}, l); + all = list_concat ({cms}, l); hunk ./src/radial.h 37 - if (is_constant(fm.x)) { + if (is_constant(fms.x)) { hunk ./src/radial.h 41 - all = list_concat ((scalar *){fm}, l); + all = list_concat ((scalar *){fms}, l); hunk ./src/radial.h 57 - scalar cmv = cm; + scalar cmsv = cms; hunk ./src/radial.h 59 - cmv[] = r*dtheta/L0; + cmsv[] = r*dtheta/L0; hunk ./src/radial.h 65 - cm[left] = dirichlet (r*dtheta/L0); - cm[right] = dirichlet (r*dtheta/L0); + cms[left] = dirichlet (r*dtheta/L0); + cms[right] = dirichlet (r*dtheta/L0); hunk ./src/radial.h 69 - The (length) metric factor `fm` is the ratio of the mapped length of + The (length) metric factor `fms` is the ratio of the mapped length of hunk ./src/radial.h 76 - face vector fmv = fm; + face vector fmsv = fms; hunk ./src/radial.h 78 - fmv.x[] = 1.; + fmsv.x[] = 1.; hunk ./src/radial.h 80 - fmv.x[] = max(r*dtheta/L0, 1e-20); + fmsv.x[] = max(r*dtheta/L0, 1e-20); hunk ./src/radial.h 82 - boundary ({cm, fm}); + boundary ({cms, fms}); hunk ./src/saint-venant.h 245 - kurganov (hm, hp, um, up, Delta*cm[]/fm.x[], &fh, &fu, &dtmax); + kurganov (hm, hp, um, up, Delta*cms[]/fms.x[], &fh, &fu, &dtmax); hunk ./src/saint-venant.h 271 - Fh.x[] = fm.x[]*fh; - Fq.x.x[] = fm.x[]*(fu - sl); - S.x[] = fm.x[]*(fu - sr); - Fq.y.x[] = fm.x[]*fv; + Fh.x[] = fms.x[]*fh; + Fq.x.x[] = fms.x[]*(fu - sl); + S.x[] = fms.x[]*(fu - sr); + Fq.y.x[] = fms.x[]*fv; hunk ./src/saint-venant.h 291 - layer[l]*(Fh.x[1,0] - Fh.x[] + Fh.y[0,1] - Fh.y[])/(cm[]*Delta); + layer[l]*(Fh.x[1,0] - Fh.x[] + Fh.y[0,1] - Fh.y[])/(cms[]*Delta); hunk ./src/saint-venant.h 294 - dhu.x[] = (Fq.x.x[] + Fq.x.y[] - S.x[1,0] - Fq.x.y[0,1])/(cm[]*Delta); + dhu.x[] = (Fq.x.x[] + Fq.x.y[] - S.x[1,0] - Fq.x.y[0,1])/(cms[]*Delta); hunk ./src/saint-venant.h 321 - double dmdl = (fm.x[1,0] - fm.x[])/(cm[]*Delta); - double dmdt = (fm.y[0,1] - fm.y[])/(cm[]*Delta); + double dmdl = (fms.x[1,0] - fms.x[])/(cms[]*Delta); + double dmdt = (fms.y[0,1] - fms.y[])/(cms[]*Delta); hunk ./src/spherical.h 29 -cm = fm.y = \cos(\phi) +cms = fms.y = \cos(\phi) hunk ./src/spherical.h 32 -fm.x = 1 +fms.x = 1 hunk ./src/spherical.h 35 -the length scale factors *fm*. - -For the volume scale factor *cm*, more care needs to be taken to +the length scale factors *fms*. + +For the volume scale factor *cms*, more care needs to be taken to hunk ./src/spherical.h 40 -\sum_{children}cm_{child} = 4cm_{parent} +\sum_{children}cms_{child} = 4cms_{parent} hunk ./src/spherical.h 45 -\Delta^2 (\sin(\phi + d\phi/2) - \sin(\phi - d\phi/2)) = \Delta^2 cm +\Delta^2 (\sin(\phi + d\phi/2) - \sin(\phi - d\phi/2)) = \Delta^2 cms hunk ./src/spherical.h 50 -static void refine_cm_lonlat (Point point, scalar cm) +static void refine_cm_lonlat (Point point, scalar cms) hunk ./src/spherical.h 54 - fine(cm,0,0) = fine(cm,1,0) = (sinphi - sin(phi - dphi))/dphi; - fine(cm,0,1) = fine(cm,1,1) = (sin(phi + dphi) - sinphi)/dphi; + fine(cms,0,0) = fine(cms,1,0) = (sinphi - sin(phi - dphi))/dphi; + fine(cms,0,1) = fine(cms,1,1) = (sin(phi + dphi) - sinphi)/dphi; hunk ./src/spherical.h 58 -static void refine_face_x_lonlat (Point point, scalar fm) +static void refine_face_x_lonlat (Point point, scalar fms) hunk ./src/spherical.h 61 - fine(fm,0,0) = fine(fm,0,1) = 1.; + fine(fms,0,0) = fine(fms,0,1) = 1.; hunk ./src/spherical.h 63 - fine(fm,2,0) = fine(fm,2,1) = 1.; - fine(fm,1,0) = fine(fm,1,1) = 1.; + fine(fms,2,0) = fine(fms,2,1) = 1.; + fine(fms,1,0) = fine(fms,1,1) = 1.; hunk ./src/spherical.h 67 -static void refine_face_y_lonlat (Point point, scalar fm) +static void refine_face_y_lonlat (Point point, scalar fms) hunk ./src/spherical.h 71 - fine(fm,0,0) = fine(fm,1,0) = cos(phi - dphi); + fine(fms,0,0) = fine(fms,1,0) = cos(phi - dphi); hunk ./src/spherical.h 73 - fine(fm,0,2) = fine(fm,1,2) = cos(phi + dphi); - fine(fm,0,1) = fine(fm,1,1) = cos(phi); + fine(fms,0,2) = fine(fms,1,2) = cos(phi + dphi); + fine(fms,0,1) = fine(fms,1,1) = cos(phi); hunk ./src/spherical.h 84 - if (is_constant(cm)) { + if (is_constant(cms)) { hunk ./src/spherical.h 86 - cm = new scalar; + cms = new scalar; hunk ./src/spherical.h 88 - all = list_concat ({cm}, l); + all = list_concat ({cms}, l); hunk ./src/spherical.h 91 - scalar cmv = cm; + scalar cmsv = cms; hunk ./src/spherical.h 94 - cmv[] = (sin(phi + dphi) - sin(phi - dphi))/(2.*dphi); + cmsv[] = (sin(phi + dphi) - sin(phi - dphi))/(2.*dphi); hunk ./src/spherical.h 97 - if (is_constant(fm.x)) { + if (is_constant(fms.x)) { hunk ./src/spherical.h 99 - fm = new face vector; + fms = new face vector; hunk ./src/spherical.h 101 - all = list_concat ((scalar *){fm}, l); + all = list_concat ((scalar *){fms}, l); hunk ./src/spherical.h 104 - face vector fmv = fm; + face vector fmsv = fms; hunk ./src/spherical.h 106 - fmv.x[] = 1.; + fmsv.x[] = 1.; hunk ./src/spherical.h 108 - fmv.y[] = cos(y*pi/180.); + fmsv.y[] = cos(y*pi/180.); hunk ./src/spherical.h 114 - cm.refine = cm.prolongation = refine_cm_lonlat; - fm.x.prolongation = refine_face_x_lonlat; - fm.y.prolongation = refine_face_y_lonlat; + cms.refine = cms.prolongation = refine_cm_lonlat; + fms.x.prolongation = refine_face_x_lonlat; + fms.y.prolongation = refine_face_y_lonlat; hunk ./src/spherical.h 119 - boundary ({cm, fm}); + boundary ({cms, fms}); [test cases and examples updated to the new *fms* variable. Jose M. Lopez-Herrera **20201221161834 Ignore-this: b87a5798a52459417e522d7db40b2210 ] conflictor [ hunk ./src/examples/karman.c 53 - muv.x[] = fm.x[]*0.125/160.; + muv.x[] = fm.x[]*0.125/Reynolds; ] : hunk ./src/examples/karman.c 53 - muv.x[] = fm.x[]*0.125/160.; + muv.x[] = fms.x[]*0.125/160.; hunk ./src/test/Makefile 152 +rising-axi-momentum.c: rising.c + ln -s rising.c rising-axi-momentum.c +rising-axi-momentum.tst: rising.s +rising-axi-momentum.tst: CFLAGS += -DAXIS=1 -DMOMENTUM=1 + hunk ./src/test/axi.c 113 - mg = viscosity (u, fm, cm, dt); + mg = viscosity (u, fms, cms, dt); hunk ./src/test/couette.c 5 -between two rotating cylinders. */ +between two rotating cylinders. If SLIP is defined, the outer cylinder +is free of viscous stresses. */ hunk ./src/test/cyl_axi.c 90 - epsilon.x[] = (ff*beta + (1. - ff))*fm.x[]; - K.x[] = cond*s.x[]*fm.x[]; + epsilon.x[] = (ff*beta + (1. - ff))*fms.x[]; + K.x[] = cond*s.x[]*fms.x[]; hunk ./src/test/cyl_planar.c 57 - epsilon.x[] = (ff*beta + (1. - ff))*fm.x[]; - K.x[] = cond*s.x[]*fm.x[]; + epsilon.x[] = (ff*beta + (1. - ff))*fms.x[]; + K.x[] = cond*s.x[]*fms.x[]; hunk ./src/test/cylinders.c 102 - mu = fm; + mu = fms; hunk ./src/test/ehd_axi_stress.c 95 - rhs[] = -rhoe[]*cm[]; + rhs[] = -rhoe[]*cms[]; hunk ./src/test/ehd_axi_stress.c 100 - epsilon.x[] = fm.x[]; - alpha.x[] = fm.x[]; + epsilon.x[] = fms.x[]; + alpha.x[] = fms.x[]; hunk ./src/test/hydrostatic2.c 76 - alpha = fm; + alpha = fs; hunk ./src/test/hydrostatic2.c 92 - gf.x[] = fm.x[] ? fm.x[]*a.x[] - alpha.x[]*(p[] - p[-1])/Delta : 0.; + gf.x[] = fms.x[] ? fms.x[]*a.x[] - alpha.x[]*(p[] - p[-1])/Delta : 0.; hunk ./src/test/hydrostatic2.c 98 - g.x[] = (gf.x[] + gf.x[1])/(fm.x[] + fm.x[1] + SEPS); + g.x[] = (gf.x[] + gf.x[1])/(fms.x[] + fms.x[1] + SEPS); hunk ./src/test/hydrostatic3.c 83 - alpha = fm; + alpha = fs; hunk ./src/test/hydrostatic3.c 104 - gf.x[] = fm.x[] ? fm.x[]*a.x[] - alpha.x[]*(p[] - p[-1])/Delta : 0.; + gf.x[] = fs.x[] ? fs.x[]*a.x[] - alpha.x[]*(p[] - p[-1])/Delta : 0.; hunk ./src/test/hydrostatic3.c 110 - g.x[] = (gf.x[] + gf.x[1])/(fm.x[] + fm.x[1] + SEPS); + g.x[] = (gf.x[] + gf.x[1])/(fs.x[] + fs.x[1] + SEPS); hunk ./src/test/lonlat.c 122 - restriction ({cm,zb,h}); // fixme: for restriction on eta + restriction ({cms,zb,h}); // fixme: for restriction on eta hunk ./src/test/neumann.c 13 +scalar csv[]; +face vector fsv[]; + hunk ./src/test/neumann.c 49 + cs = csv; + fs = fsv; hunk ./src/test/neumann.c 77 - fractions (phi, cs, fs); + fractions (phi, cs, fs); hunk ./src/test/neumann.c 84 - cm = cs; - fm = fs; + cms = cs; + fms = fs; hunk ./src/test/neumann.c 221 - draw_vof ("cs", "fs"); + draw_vof ("csv", "fsv"); hunk ./src/test/neumann.c 226 - draw_vof ("cs", "fs"); + draw_vof ("csv", "fsv"); hunk ./src/test/neumann3D.c 13 +scalar csv[]; +face vector fsv[]; + hunk ./src/test/neumann3D.c 38 + cs = csv; + fs = fsv; hunk ./src/test/neumann3D.c 57 - - cm = cs; - fm = fs; hunk ./src/test/neumann3D.c 191 - draw_vof("cs", "fs", color = "e"); + draw_vof("csv", "fsv", color = "e"); hunk ./src/test/poiseuille-axi.c 25 - mu = fm; + mu = fms; hunk ./src/test/poiseuille45.c 37 - mu = fm; + mu = fms; hunk ./src/test/porous.c 100 - mu = fm; + mu = fms; hunk ./src/test/porous1.c 125 - mu = fm; + mu = fs; hunk ./src/test/rising.c 18 +#if MOMENTUM +#include "momentum.h" +#else hunk ./src/test/rising.c 23 +#endif hunk ./src/test/rising.c 37 +#if MOMENTUM +q.t[right] = dirichlet(0); +q.t[left] = dirichlet(0); +#else hunk ./src/test/rising.c 43 +#endif hunk ./src/test/rising.c 130 - xb += x*dv; +#if MOMENTUM + vb += q.x[]*dv/rho(f[]); +#else hunk ./src/test/rising.c 134 +#endif + xb += x*dv; hunk ./src/test/rising.c 140 +#if !MOMENTUM hunk ./src/test/rising.c 144 +#endif hunk ./src/test/rising.c 192 - '../rising-axi/log' u 1:2 w l t 'Basilisk (axisymmetric)' + '../rising-axi/log' u 1:2 w l t 'Basilisk (axisymmetric)', \ + '../rising-axi-momentum/log' u 1:2 w l t 'Basilisk (momentum)' hunk ./src/test/rising.c 214 - '../rising-axi/out' u 1:5 w l t 'Basilisk (axisymmetric)' + '../rising-axi/out' u 1:5 w l t 'Basilisk (axisymmetric)', \ + '../rising-axi-momentum/out' u 1:5 w l t 'Basilisk (momentum)' addfile ./src/test/slot.c hunk ./src/test/slot.c 1 +/** +# Interface through a converging-diverging 2D slot + +The test consists in the time evolution of an interface forced to flow +through a converging-diverging nozzle thanks to an inward velocity +$U$. + */ + +//#include "grid/multigrid.h" +#include "embed.h" +#include "navier-stokes/centered.h" +#include "vof.h" + +scalar f[], * interfaces = {f}; + +/** +Flow velocity */ + +#define U 0.01 + +/** +Geometrical parameters of the slot. */ + +#define A 0.4 +#define B 0.1 + +u.n[top] = neumann(0.); +p[top] = dirichlet(0.); + +f[bottom] = 1.; +u.n[bottom] = dirichlet(U); +u.t[bottom] = dirichlet(0.); + +u.n[embed] = dirichlet(0.); +u.t[embed] = dirichlet(0.); + +int levelmax; + +int main() +{ + origin(-0.5,0.); + //DT = 0.005; + init_grid (16); + for (levelmax=5; levelmax <= 7; levelmax++) + run(); + + /** + Let's clear the log file of *WARNINGS* */ + + system("sed -i '/^W/d' log"); +} + +event init (i = 0) +{ + vertex scalar phi[]; + foreach_vertex() { + double cone = x*x-(y-A)*(y-A); + double slot = x*x-B*B; + phi[] = -intersection(cone, slot); + } + boundary ({phi}); + refine ((phi[0] + phi[1] + phi[0,1] + phi[1,1]) > 0 && level < levelmax); + fractions (phi, cs, fs); + if (levelmax == 7){ + FILE * fp = fopen ("solid", "w"); + output_facets (cs, fp, fs); + fclose (fp); + } + + /** + Initially, the interface was at height 0.1 */ + + fraction (f, 0.1 - y); + +} + +#if 0 +event adapt (i++) +{ + double uemax = 0.005; + adapt_wavelet ({f,u}, (double[]){0.001,uemax,uemax}, levelmax, 5); +} +#endif + +/** +## Outputs + +In the *update_profund* event we save the average depth of the +interface and the total volume of the denser fluid. We do not use the +height variable to calculate this depth. */ + +event update_vol (i += 5 ; t < 20) +{ + double vol = 0.; + + foreach (reduction (+:vol)) + vol += f[]*cms[]*sq(Delta); + fprintf (stderr,"%g %g \n", t, vol); + fflush(stderr); +} + +/** + +~~~gnuplot Time evolution of the volume. Continuity is preserved. +plot 'log' u 1:2 pt 5, 0.0698437 + 0.008*x lt 5 lw 2 +~~~ +*/ + +/** +We save in files the interface position at t = 10, 15 and 20. */ + +event save_depth (t = {10, 15, 20}) { + char name[50]; + sprintf (name, "depth-%g-%d", t, levelmax); + FILE *fp = fopen (name, "w"); + output_facets (f, fp); + fclose (fp); +} + +/** + +~~~gnuplot Interface position at instants t = 10, 15 and 20 for the different levels +set xrange [-0.5:0.5] +set yrange [0:1] +set size ratio -1 +set key off +plot 'solid' w l lw 2, for[i = 10:20:5] 'depth-'.i.'-5' lt 7 t 'level = 5',\ + for[i = 10:20:5] 'depth-'.i.'-6' pt 4 lc rgb "green" t 'level = 6',\ + for[i = 10:20:5] 'depth-'.i.'-7' w l lt 8 lw 2 t 'level = 7' +~~~ + +Eventually, we could save the last iteration time and/or make a movie. */ + +#if 0 +event save_dump (t = end) { + dump("dump"); +} + +#include "view.h" +event snapshot (i += 10) { + char str[99]; + sprintf (str,"In mixed cells f[] = 1 if cs[] < f[] "); + int w = 1280; + view (fov = 20.7, width = w, height = w, bg = {0.95, 0.95, 0.95}, ty = -0.5); + draw_vof("cs", lw = 3); + draw_vof("f", lc = {1,0,1}, lw = 3); + scalar ff[]; + foreach() + ff[] = cs[] > 0. ? f[] : nodata; + squares("ff", spread = -1); + cells (n = {0,1,0}); + draw_string (str, pos = 2, lw = 3); + save ("movie.mp4"); +} +#endif addfile ./src/test/slot.ref hunk ./src/test/slot.ref 1 +0 0.0696061 + res: 0.610633 sum: -167.778 nrelax: 100 +1.07188 0.0770462 +2.35755 0.0872868 +3.74502 0.0983596 +5.1666 0.109695 +6.58819 0.121034 +8.00978 0.132377 +9.43137 0.143721 +10.8824 0.155314 +12.3039 0.166662 +13.652 0.177416 +14.8652 0.18709 +16.0526 0.196568 +17.2862 0.206402 +18.5197 0.216239 +19.7533 0.226073 +0 0.0698684 + res: 0.257955 sum: -241.846 nrelax: 100 +0.724815 0.0731849 +1.31923 0.0779081 +1.96579 0.083074 +2.64236 0.0884682 +3.33699 0.0940126 +4.06196 0.0997971 +4.81625 0.105814 +5.57857 0.111894 +6.34088 0.117974 +7.1032 0.124057 +7.86552 0.130142 +8.62783 0.136223 +9.39015 0.142302 +10.1562 0.148153 +10.9375 0.153903 +11.7188 0.15958 +12.5 0.165288 +13.2812 0.171092 +13.9974 0.176398 +14.7135 0.181598 +15.4412 0.186826 +16.1535 0.191955 +16.8494 0.196968 +17.5343 0.201899 +18.2192 0.206865 +18.8552 0.211521 +19.4912 0.216176 +0 0.0699789 + res: 0.48036 sum: -1751.65 nrelax: 100 +0.342969 0.0714826 +0.600851 0.0735289 +0.870357 0.0756828 +1.14405 0.0778646 +1.42187 0.0800814 +1.70599 0.0823476 +2.0002 0.0846946 +2.30333 0.0871126 +2.61681 0.0896133 +2.94071 0.0921972 +3.2744 0.094859 +3.61834 0.0976027 +3.96907 0.100401 +4.32558 0.103245 +4.68472 0.10611 +5.04386 0.108974 +5.403 0.111839 +5.76214 0.114704 +6.12128 0.117569 +6.48042 0.120434 +6.83957 0.123299 +7.19871 0.126164 +7.55785 0.129028 +7.91699 0.131894 +8.27613 0.134758 +8.63527 0.137623 +8.98992 0.140452 +9.32661 0.143138 +9.66331 0.145824 +10 0.14851 +10.3521 0.151319 +10.7042 0.154108 +11.0517 0.156876 +11.3981 0.159648 +11.7444 0.162426 +12.0907 0.16509 +12.4371 0.167796 +12.7743 0.170451 +13.1115 0.173118 +13.4488 0.175774 +13.786 0.178407 +14.1232 0.181055 +14.4604 0.183668 +14.7977 0.186256 +15.1351 0.188805 +15.473 0.191316 +15.8068 0.193827 +16.1396 0.196321 +16.4667 0.198781 +16.7867 0.201192 +17.0992 0.203566 +17.4037 0.205881 +17.7 0.208146 +17.9875 0.210347 +18.2699 0.212518 +18.549 0.214664 +18.828 0.216814 +19.1071 0.218966 +19.3818 0.221087 +19.6394 0.223142 +19.897 0.22521 hunk ./src/test/spheres.c 110 - mu = fm; + mu = fs; hunk ./src/test/starting.c 87 - muv.x[] = fm.x[]/Re; + muv.x[] = fms.x[]/Re; hunk ./src/test/uf.c 41 - mu = fm; + mu = fms; hunk ./src/test/wannier.c 88 - mu = fm; + mu = fs; [A figure to explain better the small cell problem in vof variables Jose M. Lopez-Herrera **20201221163941 Ignore-this: c6f0f077a708e6bdfab85e30ca9b9ed3 ] addfile ./src/figures/merged.svg hunk ./src/figures/merged.svg 1 + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + 0 + 1 + 2 + + + + n + 0 + 1 + + + nx + + [changes in "axi.h" to make compatible "embed.h" Jose M. Lopez-Herrera **20210116163808 Ignore-this: dcb870645c97b02f67aeeb931501067c ] hunk ./src/axi.h 22 +#if !EMBED hunk ./src/axi.h 25 +#else + foreach_child() { + if(cs[] > 0. && cs[] < 1.) { + coord p, n = facet_normal (point, cs, fs); + double alpha = plane_alpha (cs[], n); + plane_center (n, alpha, cs[], &p); + cms[] = (y + Delta*p.y)*cs[]; + } + else + cms[] = y*cs[]; + } +#endif hunk ./src/axi.h 41 +#if !EMBED hunk ./src/axi.h 52 +#else + double sig = 0.; + double ff = 0.; + if(cs[] > 0. && cs[] < 1.) { + coord n = facet_normal (point, cs, fs); + sig = sign(n.y); + } + if (!is_refined(neighbor(-1))) { + ff = fine(fs.x,0,0); + fine(fms,0,0) = (y - Delta/4. + sig*(1.-ff))*ff; + ff = fine(fs.x,0,1); + fine(fms,0,1) = (y + Delta/4. + sig*(1.-ff))*ff; + } + if (!is_refined(neighbor(1)) && neighbor(1).neighbors) { + ff = fine(fs.x,2,0); + fine(fms,2,0) = (y - Delta/4. + sig*(1.-ff))*ff; + ff = fine(fs.x,2,1); + fine(fms,2,1) = (y + Delta/4. + sig*(1.-ff))*ff; + } + ff = fine(fs.x,1,0); + fine(fms,1,0) = (y - Delta/4. + sig*(1.-ff))*ff; + ff = fine(fs.x,1,1); + fine(fms,1,1) = (y + Delta/4. + sig*(1.-ff))*ff; +#endif hunk ./src/axi.h 80 +#if !EMBED hunk ./src/axi.h 86 +#else + if (!is_refined(neighbor(0,-1))) { + fine(fms,0,0) = (max(y - Delta/2., 1e-20))*fine(fs.y,0,0) ; + fine(fms,1,0) = (max(y - Delta/2., 1e-20))*fine(fs.y,1,0); + } + if (!is_refined(neighbor(0,1)) && neighbor(0,1).neighbors) { + fine(fms,0,2) = (y + Delta/2.)*fine(fs.y,0,2); + fine(fms,1,2) = (y + Delta/2.)*fine(fs.y,1,2); + } + fine(fms,0,1) = y*fine(fs.y,0,1); + fine(fms,1,1) = y*fine(fs.y,1,1); +#endif hunk ./src/axi.h 101 +#if EMBED +double axi_factor (Point point, coord p) { + return (y + p.y*Delta); +} +#endif + + hunk ./src/axi.h 125 - If embeded solid are presents we have to allocate a new field for - the combined metric and solid factors, *cms* beside to the already + If embeded solids are presents we have to allocate a new field for + the combined metric and solid factors, *cms*, beside to the already hunk ./src/axi.h 128 - those corresponding to solid factors. This will be done a init event */ - - -#if EMBED - if (cms.i == cs.i) + those corresponding to solid factors. */ + + +#if EMBED + metric_embed_factor = axi_factor; + if (cms.i == cs.i) { + scalar * l = list_copy (all); hunk ./src/axi.h 136 + free (all); + all = list_concat ({cs, cms}, NULL); + for (scalar s in l) + if (s.i != cms.i && s.i != cs.i) + all = list_add (all, s); + free (l); + } hunk ./src/axi.h 151 - foreach() + foreach() { +#if EMBED + if(cs[] > 0. && cs[] < 1.) { + coord p, n = facet_normal (point, cs, fs); + double alpha = plane_alpha (cs[], n); + plane_center (n, alpha, cs[], &p); + cmsv[] = (y + Delta*p.y)*cs[]; + } + else + cmsv[] = y*cs[]; +#else hunk ./src/axi.h 163 +#endif + } hunk ./src/axi.h 183 - if (fms.x.i == fs.x.i) + if (fms.x.i == fs.x.i) { + scalar * l = list_copy (all); hunk ./src/axi.h 186 + free (all); + all = list_concat ((scalar *){fs, fms}, NULL); + for (scalar s in l) + foreach_dimension() + if (s.i != fms.x.i && s.i != fs.x.i) + all = list_add (all, s); + free (l); + } hunk ./src/axi.h 198 - foreach_face() +#if EMBED + foreach_face(x) { + double sig = 0.; + if(cs[] > 0. && cs[] < 1.) { + coord n = facet_normal (point, cs, fs); + sig = sign(n.y); + } + fmsv.x[] = (y + sig*(1.-fs.x[])*Delta/2.)*fs.x[]; + } + foreach_face(y) + fmsv.y[] = max(y, 1e-20)*fs.y[]; +#else + foreach_face() { hunk ./src/axi.h 212 +#endif hunk ./src/embed.h 22 + +double cartesian_factor (Point point, coord p) { + return 1.; +} + +double (*metric_embed_factor) (Point, coord) = cartesian_factor; + hunk ./src/embed.h 691 + area *= metric_embed_factor(point, p); hunk ./src/embed.h 1010 -event metric (i = 0) +event solid (i = 0) [additional changes to embed.h + axi.h: poisson solver works Jose M. Lopez-Herrera **20210201164942 Ignore-this: 8449d77c58d843c862cb819b1d75d637 ] hunk ./src/axi.h 26 - foreach_child() { - if(cs[] > 0. && cs[] < 1.) { - coord p, n = facet_normal (point, cs, fs); - double alpha = plane_alpha (cs[], n); - plane_center (n, alpha, cs[], &p); - cms[] = (y + Delta*p.y)*cs[]; + if(cs[] > 0. && cs[] < 1.) { + coord n = interface_normal (point, cs); + // coord n = facet_normal (point, cs, fs); Better? involve fs (troubles w prolongation) + foreach_child() { + if(cs[] > 0. && cs[] < 1.) { + coord p; + double alpha = plane_alpha (cs[], n); + plane_center (n, alpha, cs[], &p); + cms[] = (y + Delta*p.y)*cs[]; + } + else + cms[] = y*cs[]; hunk ./src/axi.h 39 - else - cms[] = y*cs[]; hunk ./src/axi.h 40 + else + foreach_child() + cms[] = y*cs[]; hunk ./src/axi.h 64 - sig = sign(n.y); + sig = sign(n.y)*Delta/4.; hunk ./src/axi.h 68 - fine(fms,0,0) = (y - Delta/4. + sig*(1.-ff))*ff; + fine(fms,0,0) = (y - Delta/4. - sig*(1.-ff))*ff; hunk ./src/axi.h 70 - fine(fms,0,1) = (y + Delta/4. + sig*(1.-ff))*ff; + fine(fms,0,1) = (y + Delta/4. - sig*(1.-ff))*ff; hunk ./src/axi.h 74 - fine(fms,2,0) = (y - Delta/4. + sig*(1.-ff))*ff; + fine(fms,2,0) = (y - Delta/4. - sig*(1.-ff))*ff; hunk ./src/axi.h 76 - fine(fms,2,1) = (y + Delta/4. + sig*(1.-ff))*ff; + fine(fms,2,1) = (y + Delta/4. - sig*(1.-ff))*ff; hunk ./src/axi.h 79 - fine(fms,1,0) = (y - Delta/4. + sig*(1.-ff))*ff; + fine(fms,1,0) = (y - Delta/4. - sig*(1.-ff))*ff; hunk ./src/axi.h 81 - fine(fms,1,1) = (y + Delta/4. + sig*(1.-ff))*ff; + fine(fms,1,1) = (y + Delta/4. - sig*(1.-ff))*ff; hunk ./src/axi.h 108 +/** +If embeded solids are presents, *cms*, *fms* and fluxes need to be +update accordingly to the combined effect of axisymmetric +cylindrical coordinates and solid factors. */ + + hunk ./src/axi.h 118 + +void cms_update (scalar cms, scalar cs, face vector fs) +{ + foreach() { + if(cs[] > 0. && cs[] < 1.) { + coord p, n = facet_normal (point, cs, fs); + double alpha = plane_alpha (cs[], n); + plane_center (n, alpha, cs[], &p); + cms[] = (y + Delta*p.y)*cs[]; + } + else + cms[] = y*cs[]; + } + cms[top] = dirichlet(y*cs[]); + cms[bottom] = dirichlet(y*cs[]); +} + +void fms_update (face vector fms, scalar cs, face vector fs) +{ + foreach_face(x) { + double sig = 0.; + if(cs[] > 0. && cs[] < 1.) { + coord n = facet_normal (point, cs, fs); + sig = sign(n.y)*Delta/2.; + } + fms.x[] = (y - sig*(1.-fs.x[]))*fs.x[]; + } + foreach_face(y) + fms.y[] = max(y, 1e-20)*fs.y[]; + fms.t[top] = dirichlet(y*fs.t[]); + fms.t[bottom] = dirichlet(y*fs.t[]); +} hunk ./src/axi.h 170 - If embeded solids are presents we have to allocate a new field for - the combined metric and solid factors, *cms*, beside to the already - defined *cs*. Also the first in the list of the variables should be - those corresponding to solid factors. */ - + In case of embedded solids fluxes through them are affected by + metric factors. */ hunk ./src/axi.h 175 - if (cms.i == cs.i) { - scalar * l = list_copy (all); - cms = new scalar; - free (all); - all = list_concat ({cs, cms}, NULL); - for (scalar s in l) - if (s.i != cms.i && s.i != cs.i) - all = list_add (all, s); - free (l); - } hunk ./src/axi.h 179 - to set boundary conditions at the top and bottom so that *cm* is + to set boundary conditions at the top and bottom so that *cms* is hunk ./src/axi.h 183 - foreach() { -#if EMBED - if(cs[] > 0. && cs[] < 1.) { - coord p, n = facet_normal (point, cs, fs); - double alpha = plane_alpha (cs[], n); - plane_center (n, alpha, cs[], &p); - cmsv[] = (y + Delta*p.y)*cs[]; - } - else - cmsv[] = y*cs[]; -#else + foreach() hunk ./src/axi.h 185 -#endif - } + hunk ./src/axi.h 203 - #if EMBED - if (fms.x.i == fs.x.i) { - scalar * l = list_copy (all); - fms = new face vector; - free (all); - all = list_concat ((scalar *){fs, fms}, NULL); - for (scalar s in l) - foreach_dimension() - if (s.i != fms.x.i && s.i != fs.x.i) - all = list_add (all, s); - free (l); - } -#endif - - hunk ./src/axi.h 204 -#if EMBED - foreach_face(x) { - double sig = 0.; - if(cs[] > 0. && cs[] < 1.) { - coord n = facet_normal (point, cs, fs); - sig = sign(n.y); - } - fmsv.x[] = (y + sig*(1.-fs.x[])*Delta/2.)*fs.x[]; - } - foreach_face(y) - fmsv.y[] = max(y, 1e-20)*fs.y[]; -#else - foreach_face() { + foreach_face() hunk ./src/axi.h 206 -#endif hunk ./src/embed.h 546 - fa += fs.x[] + fs.x[1]; + fa += fms.x[] + fms.x[1]; hunk ./src/embed.h 709 - fa += fs.x[] + fs.x[1]; + fa += fms.x[] + fms.x[1]; hunk ./src/embed.h 1005 -the metric using the embedded fractions. The variables are not longer -constant fields and they must be allocated. Also, the solid fields -are moved to the beginning of the list of variables. This ensures that -solid fields are the first to be updated. */ - -event solid (i = 0) +the embed solid using the embedded fractions. The variables are not longer +constant fields and they must be allocated. + +The variables *cms* amd *fms* include both the metric and embed solid +factors. These variables are the solid factors *cs* amd *fs* if the +coordinate system is cartesian (absence of any metric factor). On the +contrary, with non cartesian coordinate systems *cs* amd *fs* must be +different from *cms* amd *fms*. Also, the solid fields are moved to +the beginning of the list of variables. This ensures that solid fields +are the first to be updated. */ + +event metric (i = 0) hunk ./src/embed.h 1022 - all = list_concat ((scalar *){fs}, l); + if (is_constant(fms.x)) { + all = list_concat ((scalar *){fs}, l); + fms = fs; + } + else { + all = list_concat ((scalar *){fs, fms}, NULL); + for (scalar s in l) + foreach_dimension() + if (s.i != fms.x.i && s.i != fs.x.i) + all = list_add (all, s); + } hunk ./src/embed.h 1036 - face vector fsv = fs; hunk ./src/embed.h 1037 - fsv.x[] = 1.; + fs.x[] = 1.; hunk ./src/embed.h 1043 - all = list_concat ({cs}, l); + if (is_constant(cms)) { + all = list_concat ({cs}, l); + cms = cs; + } + else { + all = list_concat ({cs, cms}, NULL); + for (scalar s in l) + if (s.i != cms.i && s.i != cs.i) + all = list_add (all, s); + } hunk ./src/embed.h 1056 - scalar csv = cs; hunk ./src/embed.h 1057 - csv[] = 1.; + cs[] = 1.; hunk ./src/embed.h 1080 - - /** - The variables *cms* amd *fms* include both the metric and solid - factors. These variables are intialized to the solid factors - here. The metric factors would be added elsewhere (see for example - [axi.h](axi.h#mg_solve)) */ - - cms = cs; - fms = fs; - hunk ./src/test/Makefile 434 +dirichlet-axi.tst: CFLAGS += -DDIRICHLET=1 +dirichlet-axi.c: neumann-axi.c + ln -s neumann-axi.c dirichlet-axi.c + hunk ./src/test/Makefile 444 -embed-tests: neumann.tst dirichlet.tst poiseuille45.tst \ +embed-tests: neumann.tst dirichlet.tst neumann-axi.tst dirichlet-axi.tst \ + poiseuille45.tst \ addfile ./src/test/neumann-axi.c hunk ./src/test/neumann-axi.c 1 +/** +# Axisymmetric Poisson equation on complex domains + +The axisymmetric version of [neumann.c](src/neumann.c). */ + +#include "embed.h" +#include "axi.h" +#include "poisson.h" +#include "utils.h" +#include "view.h" + +scalar csv[], cmsv[]; +face vector fsv[], fmsv[]; + +/** +The exact solution is given by +$$ +\phi(r,\theta) = r^4\cos 3\theta +$$ +*/ + +static double exact (double x, double y) { + double theta = atan2 (y, x), r2 = x*x + y*y; + return sq(r2)*cos (3.*theta); +} + +double exact_gradient (Point point, double theta, double r) +{ + coord n = facet_normal (point, cs, fs); + normalize (&n); + double dsdtheta = - 3.*cube(r)*sin (3.*theta); + double dsdr = 4.*cube(r)*cos (3.*theta); + return (n.x*(dsdr*cos(theta) - dsdtheta*sin(theta)) + + n.y*(dsdr*sin(theta) + dsdtheta*cos(theta))); +} + +/** +We also need a function for homogeneous Dirichlet condition. */ + +static +double dirichlet_homogeneous_bc (Point point, Point neighbor, + scalar s, void * data) { + return 0.; +} + +int main() +{ + cms = cmsv; + fms = fmsv; + cs = csv; + fs = fsv; + for (N = 32; N <= 512; N *= 2) { + origin (-L0/2., 0.); + init_grid (N); + + /** + The shape of the domain is given by + $$ + \Omega = {(r,\theta): r \leq 0.30 + 0.15\cos 6\theta} + $$ + for Problem 1 and + $$ + \Omega = {(r,\theta): r \geq 0.25 + 0.05\cos 6\theta} + $$ + for Problem 2. + */ + + vertex scalar phi[]; + foreach_vertex() { + double theta = atan2(y, x), r = sqrt (sq(x) + sq(y)); +#if DIRICHLET + phi[] = 0.30 + 0.15*cos(6.*theta) - r; +#else + phi[] = r - (0.25 + 0.05*cos(6.*theta)); +#endif + } + boundary ({phi}); + fractions (phi, cs, fs); + cms_update (cms, cs, fs); + fms_update (fms, cs, fs); +#if TREE + cms.refine = cms.prolongation = refine_cm_axi; + cs.refine = cs.prolongation = fraction_refine; + fms.x.refine = refine_face_x_axi; + fms.y.refine = refine_face_y_axi; + metric_embed_factor = axi_factor; +#endif + boundary ({cs, fs, cms, fms}); + restriction ({cs, fs, cms, fms}); + + /** + Conditions on the box boundaries are set (only relevant for Problem 3). */ + + scalar a[], b[]; + + a[left] = exact (x - Delta/2., y); + a.boundary_homogeneous[left] = dirichlet_homogeneous_bc; + a[right] = exact (x + Delta/2., y); + a.boundary_homogeneous[right] = dirichlet_homogeneous_bc; + a[top] = exact (x, y + Delta/2.); + a.boundary_homogeneous[top] = dirichlet_homogeneous_bc; + a[bottom] = exact (x, y - Delta/2.); + a.boundary_homogeneous[bottom] = dirichlet_homogeneous_bc; + + /** + The boundary conditions on the embedded boundary are Dirichlet and + Neumann for Problem 1 and 3, respectively. */ + +#if DIRICHLET + a[embed] = dirichlet (exact (x,y)); +#else + a[embed] = neumann (exact_gradient (point, atan2(y, x), sqrt(x*x + y*y))); +#endif + + /** + We use "third-order" [face flux interpolation](/src/embed.h). */ + + a.third = true; + + /** + The right-hand-side + $$ + \Delta\phi = 2r^2 (-3 \cos \theta + 4 \cos 3\theta) + $$ + is defined using the coordinates of the + barycenter of the cut cell (xc,yc), which is calculated from the + cell and surface fractions. */ + + foreach() { + a[] = cs[] > 0. ? exact (x, y) : nodata; + + double xc = x, yc = y; + if (cs[] > 0. && cs[] < 1.) { + coord n = facet_normal (point, cs, fs), p; + double alpha = plane_alpha (cs[], n); + line_center (n, alpha, cs[], &p); + xc += p.x*Delta, yc += p.y*Delta; + } + // fprintf (stderr, "xc %g %g\n", xc, yc); + double theta = atan2(yc, xc), r2 = sq(xc) + sq(yc); + b[] = 2.*r2*(-3* cos(theta) + 4.* cos (3.*theta))*cms[]; + } + boundary ({a,b}); + +#if 0 + // output_cells (stdout); + // output_facets (cs, stdout, fs); + + scalar e[]; + foreach() { + if (cs[] > 0. && cs[] < 1.) { + scalar s = a; + coord n = facet_normal (point, cs, fs), p; + double alpha = plane_alpha (cs[], n); + double length = line_length_center (n, alpha, &p); + x += p.x*Delta, y += p.y*Delta; + double theta = atan2(y,x), r = sqrt(x*x + y*y); + + double dsdtheta = - 3.*cube(r)*sin (3.*theta); + double dsdr = 4.*cube(r)*cos (3.*theta); + double nn = sqrt (sq(n.x) + sq(n.y)); + n.x /= nn, n.y /= nn; + double dsdn = (n.x*(dsdr*cos(theta) - dsdtheta*sin(theta)) + + n.y*(dsdr*sin(theta) + dsdtheta*cos(theta))); + + e[] = dsdn - dirichlet_gradient (point, s, cs, n, p, exact (x, y)); +#if 1 + fprintf (stderr, "g %g %g %g %g\n", + x, y, dsdn, + dirichlet_gradient (point, s, cs, n, p, exact (x, y))); +#endif + } + else + e[] = nodata; + } + + norm n = normf (e); + fprintf (stderr, "%d %g %g\n", + N, n.rms, n.max); +#else + + /** + The Poisson equation is solved. */ + + struct Poisson p; + p.alpha = fms; + p.lambda = zeroc; + p.embed_flux = embed_flux; + scalar res[]; + double maxp = residual ({a}, {b}, {res}, &p), maxf = 0.; + foreach() + if (cs[] == 1. && fabs(res[]) > maxf) + maxf = fabs(res[]); + fprintf (stderr, "maxres %d %.3g %.3g\n", N, maxf, maxp); + + // FIXME: need to set minlevel to 4 + timer t = timer_start(); + mgstats s = poisson (a, b, alpha = fms, tolerance = 1e-6, minlevel = 4); + double dt = timer_elapsed (t); + printf ("%d %g %d %d\n", N, dt, s.i, s.nrelax); + + /** + The total (*e*), partial cells (*ep*) and full cells (*ef*) errors + fields and their norms are computed. */ + + scalar e[], ep[], ef[]; + foreach() { + if (cs[] == 0.) + ep[] = ef[] = e[] = nodata; + else { + e[] = a[] - exact (x, y); + ep[] = cs[] < 1. ? e[] : nodata; + ef[] = cs[] >= 1. ? e[] : nodata; + } + } + norm n = normf (e), np = normf (ep), nf = normf (ef); + fprintf (stderr, "%d %.3g %.3g %.3g %.3g %.3g %.3g %d %d\n", + N, n.avg, n.max, np.avg, np.max, nf.avg, nf.max, s.i, s.nrelax); + + /** + The solution is displayed using bview. */ + + view (fov = 18, ty = -0.45); + squares ("a", spread = -1); + draw_vof ("csv", "fsv"); + save ("a.png"); + + clear(); + squares ("e", spread = -1); + draw_vof ("csv", "fsv"); + save ("e.png"); +#endif + + dump ("dump"); + } +} + +/** +## Results + +### Problem 1 + +![Solution on the finest mesh](dirichlet-axi/a.png) + +![Error on the finest mesh](dirichlet-axi/e.png) + +For Dirichlet boundary conditions, the residual converges at first order +in partial cells. + +~~~gnuplot Maximum residual convergence +set xrange [*:*] +ftitle(a,b) = sprintf("%.3f/x^{%4.2f}", exp(a), -b) +f(x) = a + b*x +fit f(x) '< grep maxres ../dirichlet-axi/log' u (log($2)):(log($3)) via a,b +f2(x) = a2 + b2*x +fit f2(x) '' u (log($2)):(log($4)) via a2,b2 +set xlabel 'Resolution' +set logscale +set cbrange [1:2] +set xtics 16,2,1024 +set grid ytics +set ytics format "% .0e" +set xrange [16:1024] +set ylabel 'Maximum residual' +set yrange [1e-6:] +set key bottom left +plot '' u 2:3 t 'full cells', exp(f(log(x))) t ftitle(a,b), \ + '' u 2:4 t 'partial cells', exp(f2(log(x))) t ftitle(a2,b2), \ + 'neumann.table1' u 1:4 w lp t 'full cells (J and C, 1998)', \ + 'neumann.table1' u 1:2 w lp t 'partial cells (J and C, 1998)' +~~~ + +This leads to third-order overall convergence. + +~~~gnuplot Error convergence +set xrange [*:*] +fit f(x) '../dirichlet-axi/log' u (log($1)):(log($3)) via a,b +fit f2(x) '' u (log($1)):(log($4)) via a2,b2 +f3(x) = a3 + b3*x +fit f3(x) '' u (log($1)):(log($6)) via a3,b3 +set xrange [16:1024] +set ylabel 'Error' +set yrange [*:*] +plot '' u 1:3 pt 6 t 'max all cells', exp(f(log(x))) t ftitle(a,b), \ + '' u 1:4 t 'avg partial cells', exp(f2(log(x))) t ftitle(a2,b2), \ + '' u 1:6 t 'avg full cells', exp(f3(log(x))) t ftitle(a3,b3), \ + 'neumann.table2' u 1:2 w lp t 'max all cells (J and C, 1998)',\ + '' u 1:3 w lp t 'avg partial cells (J and C, 1998)', \ + '' u 1:4 w lp t 'avg full cells (J and C, 1998)' +~~~ + +### Problem 3 + +![Solution on the finest mesh](neumann-axi/a.png) + +![Error on the finest mesh](neumann-axi/e.png) + +For Neumann boundary conditions, the residual also converges at first order +in partial cells. + +~~~gnuplot Maximum residual convergence +set xrange [*:*] +fit f(x) '< grep maxres log' u (log($2)):(log($3)) via a,b +fit f2(x) '' u (log($2)):(log($4)) via a2,b2 +set ylabel 'Maximum residual' +set xrange [16:1024] +set yrange [1e-6:] +set key bottom left +plot '' u 2:3 t 'full cells', exp(f(log(x))) t ftitle(a,b), \ + '' u 2:4 t 'partial cells', exp(f2(log(x))) t ftitle(a2,b2), \ + 'neumann.table5' u 1:2 w lp t 'full cells (J and C, 1998)', \ + 'neumann.table5' u 1:3 w lp t 'partial cells (J and C, 1998)' +~~~ + +But this now leads to second-order overall convergence. + +~~~gnuplot Maximum error convergence +set xrange [*:*] +fit f(x) 'log' u (log($1)):(log($3)) via a,b +fit f2(x) '' u (log($1)):(log($5)) via a2,b2 +set ylabel 'Maximum error' +set xrange [16:1024] +set yrange [*:*] +plot '' u 1:3 pt 6 t 'all cells', exp(f(log(x))) t ftitle(a,b), \ + '' u 1:5 t 'partial cells', exp(f2(log(x))) t ftitle(a2,b2), \ + 'neumann.table5' u 1:5 w lp t 'all cells (J and C, 1998)', \ + 'neumann.table5' u 1:4 w lp t 'partial cells (J and C, 1998)' +~~~ + +## References + +~~~bib +@article{johansen1998, + title={A Cartesian grid embedded boundary method for Poisson's + equation on irregular domains}, + author={Johansen, Hans and Colella, Phillip}, + journal={Journal of Computational Physics}, + volume={147}, + number={1}, + pages={60--85}, + year={1998}, + publisher={Elsevier}, + url={https://pdfs.semanticscholar.org/17cd/babecd054d58da05c2ba009cccb3c687f58f.pdf} +} +~~~ +*/ [Fixed minimal conflicts in karman.c Jose M. Lopez-Herrera **20210224120619 Ignore-this: 6436c656ad02c0c06b5c2fb08be3a54c ] hunk ./src/examples/karman.c 53 - muv.x[] = fm.x[]*0.125/160.; + muv.x[] = fms.x[]*0.125/Reynolds; Context: [Default display for various solvers Stephane Popinet **20210131171057 Ignore-this: 589f235246e148b4ba8c6ceec6e07fc6 ] [TAG release 20-11-13 Stephane Popinet **20210201152605 Ignore-this: 8d63a681436472b586c079ca23ca60a2 ] Patch bundle hash: 8c4d15cb863819d97aa4223ae661a8060543f816