src/test/neumann-axi.c
Axisymmetric Poisson equation on complex domains
The axisymmetric version of neumann.c.
The exact solution is given by \displaystyle \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.;
}
scalar cmv[];
face vector fmv[];
int main()
{
cm = cmv, fm = fmv;
for (N = 32; N <= 512; N *= 2) {
size (1. [0]); // dimensionless space
origin (-L0/2., 0.);
init_grid (N);The shape of the domain is given by \displaystyle \Omega = {(r,\theta): r \leq 0.30 + 0.15\cos 6\theta} for Problem 1 and \displaystyle \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
}
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
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)));
#endifWe use “third-order” face flux interpolation.
a.third = true;The right-hand-side \displaystyle \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[];
}
#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);
#elseThe 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 ("cs", "fs");
save ("a.png");
clear();
squares ("e", spread = -1);
draw_vof ("cs", "fs");
save ("e.png");
#endif
dump ("dump");
}
}Results
Problem 1
Solution on the finest mesh
Error on the finest mesh
For Dirichlet boundary conditions, the residual converges at first order in partial cells.
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)'
Maximum residual convergence (script)
This leads to third-order overall 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)'
Error convergence (script)
Problem 3
Solution on the finest mesh
Error on the finest mesh
For Neumann boundary conditions, the residual also converges at first order in partial cells.
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)'
Maximum residual convergence (script)
But this now leads to second-order overall 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)'
Maximum error convergence (script)
References
| [johansen1998] |
Hans Johansen and Phillip Colella. A cartesian grid embedded boundary method for poisson’s equation on irregular domains. Journal of Computational Physics, 147(1):60–85, 1998. [ .pdf ] |
