5 patches for repository http://basilisk.fr/basilisk: patch 542ddc02417c6fbee806dbd3e17af3a431ca4fe4 Author: Jose M. Lopez-Herrera Date: Tue Mar 2 14:13:53 CET 2021 * This figure explain better changes in vof.h patch 28b8442366b00333b1ef52d693a7697ca885fd98 Author: Jose M. Lopez-Herrera Date: Tue Mar 2 14:19:45 CET 2021 * Metric + solid + vof: Major changes in embed.h and axi.h patch 09e87ada0cb9345bba411ae3c14dfb3975ec33c2 Author: Jose M. Lopez-Herrera Date: Tue Mar 2 14:25:42 CET 2021 * Making embed + vof + metric works: changes in vof.h (and common files) patch 49abd1fb6292ce6525ef313d4779afebc87ec9eb Author: Jose M. Lopez-Herrera Date: Tue Mar 2 14:32:45 CET 2021 * Minor changes in tension.h, iforces.h and contact.h to avoid division by 0 In addition some metric factors lacking in all-mach & momentum.h are completed patch 5caa133aa624dabae066749db1b761402d519a25 Author: Jose M. Lopez-Herrera Date: Wed Mar 3 14:07:11 CET 2021 * test cases for recent changes: some just modifications & some new. New patches: [This figure explain better changes in vof.h Jose M. Lopez-Herrera **20210302131353 Ignore-this: dd13c196400747633c9944ed3ce23999 ] addfile ./src/figures/merged.svg hunk ./src/figures/merged.svg 1 + + + + [Metric + solid + vof: Major changes in embed.h and axi.h Jose M. Lopez-Herrera **20210302131945 Ignore-this: cebcb7904460de3e451f38183c50666b ] hunk ./src/axi.h 22 +#if !EMBED hunk ./src/axi.h 25 +#else + 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); + cm[] = (y + Delta*p.y)*cs[]; + } + else + cm[] = y*cs[]; + } + } + else + foreach_child() + cm[] = y*cs[]; +#endif hunk ./src/axi.h 48 +#if !EMBED hunk ./src/axi.h 59 +#else + double sig = 0.; + double ff = 0.; + if(cs[] > 0. && cs[] < 1.) { + coord n = facet_normal (point, cs, fs); + sig = sign(n.y)*Delta/4.; + } + if (!is_refined(neighbor(-1))) { + ff = fine(fs.x,0,0); + fine(fm,0,0) = (y - Delta/4. - sig*(1.-ff))*ff; + ff = fine(fs.x,0,1); + fine(fm,0,1) = (y + Delta/4. - sig*(1.-ff))*ff; + } + if (!is_refined(neighbor(1)) && neighbor(1).neighbors) { + ff = fine(fs.x,2,0); + fine(fm,2,0) = (y - Delta/4. - sig*(1.-ff))*ff; + ff = fine(fs.x,2,1); + fine(fm,2,1) = (y + Delta/4. - sig*(1.-ff))*ff; + } + ff = fine(fs.x,1,0); + fine(fm,1,0) = (y - Delta/4. - sig*(1.-ff))*ff; + ff = fine(fs.x,1,1); + fine(fm,1,1) = (y + Delta/4. - sig*(1.-ff))*ff; +#endif hunk ./src/axi.h 87 +#if !EMBED hunk ./src/axi.h 93 +#else + if (!is_refined(neighbor(0,-1))) { + fine(fm,0,0) = (max(y - Delta/2., 1e-20))*fine(fs.y,0,0) ; + fine(fm,1,0) = (max(y - Delta/2., 1e-20))*fine(fs.y,1,0); + } + if (!is_refined(neighbor(0,1)) && neighbor(0,1).neighbors) { + fine(fm,0,2) = (y + Delta/2.)*fine(fs.y,0,2); + fine(fm,1,2) = (y + Delta/2.)*fine(fs.y,1,2); + } + fine(fm,0,1) = y*fine(fs.y,0,1); + fine(fm,1,1) = y*fine(fs.y,1,1); +#endif hunk ./src/axi.h 108 +/** +If embeded solids are presents, *cm*, *fm* and fluxes need to be +update accordingly to the combined effect of axisymmetric +cylindrical coordinates and solid factors. */ + + +#if EMBED +double axi_factor (Point point, coord p) { + return (y + p.y*Delta); +} + +void cm_update (scalar cm, 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); + cm[] = (y + Delta*p.y)*cs[]; + } + else + cm[] = y*cs[]; + } + cm[top] = dirichlet(y*cs[]); + cm[bottom] = dirichlet(y*cs[]); +} + +void fm_update (face vector fm, 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.; + } + fm.x[] = (y - sig*(1.-fs.x[]))*fs.x[]; + } + foreach_face(y) + fm.y[] = max(y, 1e-20)*fs.y[]; + fm.t[top] = dirichlet(y*fs.t[]); + fm.t[bottom] = dirichlet(y*fs.t[]); +} +#endif + + hunk ./src/axi.h 170 + In case of embedded solids fluxes through them are affected by + metric factors. */ + +#if EMBED + metric_embed_factor = axi_factor; +#endif + + /** hunk ./src/axi.h 185 + hunk ./src/axi.h 202 + hunk ./src/embed.h 19 -scalar cs[]; -face vector fs[]; +extern scalar cs; +extern face vector fs; + + +double cartesian_factor (Point point, coord p) { + return 1.; +} + +double (*metric_embed_factor) (Point, coord) = cartesian_factor; hunk ./src/embed.h 546 - fa += fs.x[] + fs.x[1]; + fa += fm.x[] + fm.x[1]; hunk ./src/embed.h 691 + area *= metric_embed_factor(point, p); hunk ./src/embed.h 709 - fa += fs.x[] + fs.x[1]; + fa += fm.x[] + fm.x[1]; hunk ./src/embed.h 712 - return - mua/(fa + SEPS)*coef*area/Delta;; + return - mua/(fa + SEPS)*coef*area/Delta; hunk ./src/embed.h 836 - already stores the product $f_su$. */ + already stores the product $u_f \, f_m$. */ hunk ./src/embed.h 844 - double dtmax = Delta*cs[]/(umax + SEPS); + double dtmax = Delta*cm[]/(umax + SEPS); hunk ./src/embed.h 852 - F /= Delta*cs[]; + F /= Delta*cm[]; hunk ./src/embed.h 874 - scs += sq(cs[]); - e[] = (dt - dtmax)*F*cs[]/scs; + scs += sq(cm[]); //Correct? + e[] = (dt - dtmax)*F*cm[]/scs; hunk ./src/embed.h 894 +### 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[]))/(cm[]*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[]))/(cm[]*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*cm[]/(umax + SEPS); + + if ( dt <= dtmax) { + c[] += dt*(flux[] - flux[1] + + cc[]*(uf.x[1] - uf.x[]))/(cm[]*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[]))/(cm[]*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[]))/(cm[]*Delta + SEPS); + c[i] += ((dt -dtmax)*(flux[] - flux[i]) + + (uf.x[1] - uf.x[])*(cc[]*dtmax -cc[i]*dt))/(cm[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[]))/(cm[]*Delta + SEPS); + t[i] += ((dt-dtmax)*(tflux[] - tflux[i]) + + (uf.x[1] - uf.x[])*(tc[]*dtmax -tc[i]*dt))/(cm[i]*Delta + SEPS); + } + } + } + } + } + } + +/** hunk ./src/embed.h 1005 -the metric using the embedded fractions. */ +the embed solid using the embedded fractions. The variables are not longer +constant fields and they must be allocated. + +The variables *cm* amd *fm* 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 *cm* amd *fm*. 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 1018 + if (is_constant(fs.x)) { + scalar * l = list_copy (all); + fs = new face vector; + free (all); + if (is_constant(fm.x)) { + all = list_concat ((scalar *){fs}, l); + fm = fs; + } + else { + all = list_concat ((scalar *){fs, fm}, NULL); + for (scalar s in l) + foreach_dimension() + if (s.i != fm.x.i && s.i != fs.x.i) + all = list_add (all, s); + } + free (l); + } + + foreach_face() + fs.x[] = 1.; + + if (is_constant(cs)) { + scalar * l = list_copy (all); + cs = new scalar; + free (all); + if (is_constant(cm)) { + all = list_concat ({cs}, l); + cm = cs; + } + else { + all = list_concat ({cs, cm}, NULL); + for (scalar s in l) + if (s.i != cm.i && s.i != cs.i) + all = list_add (all, s); + } + free (l); + } + hunk ./src/embed.h 1058 - foreach_face() - fs.x[] = 1.; + hunk ./src/embed.h 1080 - - // fixme: embedded boundaries cannot be combined with (another) metric yet - assert (is_constant (cm) || cm.i == cs.i); - - cm = cs; - fm = fs; [Making embed + vof + metric works: changes in vof.h (and common files) Jose M. Lopez-Herrera **20210302132542 Ignore-this: 21caa15538fd10ee124ec74836b1a251 ] hunk ./src/common.h 1209 +(const) face vector fs = unityf; +(const) scalar cs = unity; + hunk ./src/grid/tree-common.h 173 - scalar * listr = list_concat ({cm}, p.slist); + scalar * listr = (cm.i == cs.i) ? list_concat ({cm}, p.slist) : + list_concat ({cs, cm}, p.slist); hunk ./src/grid/tree-mpi.h 1251 - scalar * listm = is_constant(cm) ? NULL : (scalar *){fm}; + scalar * listm = NULL; + if (!is_constant(cm)) + listm = (cm.i != cs.i) ? (scalar *) {fs, fm} : (scalar *) {fm}; hunk ./src/grid/tree-mpi.h 1281 - } + } 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/output.h 1030 - scalar * list = is_constant(cm) ? NULL : list_concat ({cm}, NULL); + scalar * list = NULL; + if (!is_constant (cm)) + list = (cm.i != cs.i) ? list_concat ({cs, cm}, NULL) : list_concat ({cm}, NULL); hunk ./src/output.h 1034 - if (!s.face && !s.nodump && s.i != cm.i) + if (!s.face && !s.nodump && s.i != cm.i && s.i != cs.i) hunk ./src/output.h 1307 + + scalar * listm = NULL; + if (!is_constant (cm)) + listm = (cm.i != cs.i) ? (scalar *) {fs, fm} : (scalar *) {fm}; hunk ./src/output.h 1312 - scalar * listm = is_constant(cm) ? NULL : (scalar *){fm}; 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 242 +#if EMBED + ci *= cs[i]; +#endif hunk ./src/vof.h 320 +#if !EMBED hunk ./src/vof.h 328 +#else + vof_update_x (c, cc, flux, tracers, tcl, tfluxl, uf, dt); +#endif [Minor changes in tension.h, iforces.h and contact.h to avoid division by 0 Jose M. Lopez-Herrera **20210302133245 Ignore-this: 66e91c4e1f75481ea7b77b6d0b830021 In addition some metric factors lacking in all-mach & momentum.h are completed ] hunk ./src/all-mach.h 72 + + if (rho.i == unity.i) + rho = cm; 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[])*cm[]/rho[]; hunk ./src/all-mach.h 167 - q.x[] = q.x[]*rho[] - dt*g.x[]; + q.x[] = q.x[]*rho[]/cm[] - dt*g.x[]; hunk ./src/all-mach.h 169 - } + } 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 264 - g.x[] = rho[]*(gf.x[] + gf.x[1])/2.; + g.x[] = rho[]*(gf.x[] + gf.x[1])/(2.*cm[]); hunk ./src/contact.h 12 + +#undef interface_normal hunk ./src/iforce.h 108 - ia.x[] += alpha.x[]/fm.x[]*phif*(f[] - f[-1])/Delta; + ia.x[] += alpha.x[]/(fm.x[] + SEPS)*phif*(f[] - f[-1])/Delta; hunk ./src/momentum.h 44 - rhov[] = rho(f[]); + rhov[] = rho(f[])*cm[]; hunk ./src/momentum.h 47 - alphav.x[] = 2./(rhov[] + rhov[-1]); hunk ./src/momentum.h 48 + alphav.x[] = fm.x[]/rho(ff); hunk ./src/navier-stokes/conserving.h 19 - double du = u[] - rhou/((1 << dimension)*cm[]*rho(f[])); + double du = u[] - rhou/((1 << dimension)*(cm[] + SEPS)*rho(f[])); hunk ./src/navier-stokes/conserving.h 29 - u[] = rhou/((1 << dimension)*cm[]*rho(f[])); + u[] = rhou/((1 << dimension)*(cm[] + SEPS)*rho(f[])); 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 (fm.x[] > 0.) { + 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; + } 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/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)' [test cases for recent changes: some just modifications & some new. Jose M. Lopez-Herrera **20210303130711 Ignore-this: 47892e60342dceb3edaf228e33e52a3a ] 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 + 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[], cmv[]; +face vector fsv[], fmv[]; + +/** +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() +{ + cm = cmv; + fm = fmv; + 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); + cm_update (cm, cs, fs); + fm_update (fm, cs, fs); +#if TREE + cm.refine = cm.prolongation = refine_cm_axi; + cs.refine = cs.prolongation = fraction_refine; + fm.x.refine = refine_face_x_axi; + fm.y.refine = refine_face_y_axi; + metric_embed_factor = axi_factor; +#endif + boundary ({cs, fs, cm, fm}); + restriction ({cs, fs, cm, fm}); + + /** + 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))*cm[]; + } + 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 = fm; + 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 = fm, 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} +} +~~~ +*/ 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 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 194 - draw_vof("cs", "fs", color = "e"); + draw_vof("csv", "fsv", color = "e"); 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[]*cm[]*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 Context: [TAG release 20-11-13 Stephane Popinet **20210201152605 Ignore-this: 8d63a681436472b586c079ca23ca60a2 ] Patch bundle hash: 59cab9876c84fed81548df29db66fcd81b6f0ade