src/test/rotation.c
Pure rotation of a smooth tracer field
set output 'error.png'
set pm3d
set pm3d map interpolate 1,1
set palette defined ( 0 0 0 0.5647, 0.125 0 0.05882 1, 0.25 0 0.5647 1,\
0.375 0.05882 1 0.9333, 0.5 0.5647 1 0.4392, 0.625 1 0.9333 0, \
0.75 1 0.4392 0, 0.875 0.9333 0 0, 1 0.498 0 0 )
set size ratio -1
splot 'out' t ""
reset
ftitle(a,b) = sprintf("%.0f/x^{%4.2f}", exp(a), -b)
f(x)=a+b*x
fit f(x) 'log' u (log($1)):(log($4)) via a,b
f2(x)=a2+b2*x
fit f2(x) 'log' u (log($1)):(log($2)) via a2,b2
set xlabel 'Maximum resolution'
set ylabel 'Maximum error'
set logscale
set cbrange [1:2]
set xrange [32:512]
set xtics 32,2,512
set grid ytics
plot 'log' u 1:4 t 'max', 'log' u 1:2 t 'norm1', \
exp(f(log(x))) t ftitle(a,b), exp(f2(log(x))) t ftitle(a2,b2)
We use the advection solver with a single tracer f
.
#include "grid/cartesian.h"
#include "advection.h"
scalar f[];
scalar * tracers = {f};
We impose no normal flux on the box boundaries.
u.n[left] = 0.;
u.n[right] = 0.;
u.n[top] = 0.;
u.n[bottom] = 0.;
int main()
{
// coordinates of lower-left corner
origin (-0.5, -0.5);
// maximum timestep
DT = .1;
// CFL number
CFL = 0.8;
The spatial resolution varies from 64 to 256 points per box size to study spatial convergence.
for (N = 64; N <= 256; N *= 2)
run ();
}
The initial tracer field is a somewhat complicated but compact and smooth function.
double bump (double x, double y)
{
double k = 20., a = 20000.;
double r2 = x*x + y*y;
double coeff = 20. + a*r2*r2*r2*r2;
return (1. + cos(k*x)*cos(k*y))*exp(-coeff*r2)/2.;
}
event init (i = 0)
{
foreach()
f[] = bump(x,y);
}
The advection velocity field is a simple solid rotation.
event velocity (i++) {
trash ({u});
foreach_face(x) u.x[] = -8.*y;
foreach_face(y) u.y[] = 8.*x;
}
We log statistics on mass conservation and errors after one full rotation.
#define end 0.785398
event logfile (t = {0,end}) {
stats s = statsf (f);
fprintf (stderr, "# %f %.12f %g %g\n", t, s.sum, s.min, s.max);
}
event field (t = end) {
scalar e[];
foreach()
e[] = f[] - bump(x,y);
norm n = normf (e);
fprintf (stderr, "%d %g %g %g\n", N, n.avg, n.rms, n.max);
if (N == 256)
output_field ({e}, stdout, N);
}