sandbox/Antoonvh/ABL/example_rough.c
Flow over a rough embedded topography
~ {.bash}
set size ratio -1
plot 'ustar50' matrix with image
#include "grid/octree.h"
#include "embed.h"
#include "navier-stokes/centered.h"
#include "view.h"
#include "lambda2.h"
#include "tracer.h"
#include "diffusion.h"
scalar s[], * tracers = {s};
Including rough-embed.h
switches the embedded boundary to a rough surface whose roughness length is defined via Z_0
.
double z0 = 0.01;
#define Z_0 (z0)
#include "rough-embed.h"
face vector muc[];
double ue = 0.05, cse = 0.01; //Refinement criteria
int maxlevel = 6;
int main() {
N = 32;
L0 = 2*pi;
X0 = Z0 = -L0/2.;
periodic (left);
periodic (front);
mu = muc;
run();
}
event init (t = 0) {
add_rough_scalar (s, s_val = 1); //<- set boundary conditions for s
vertex scalar phi[];
refine (y < 2.1 && y > 0.1 && level < maxlevel);
foreach_vertex()
phi[] = -1.1 + 0.5*sin(x) * cos(z*2.) + y;
boundary ({phi});
fractions (phi, cs, fs);
foreach()
u.x[] = cs[] > 0;
boundary ({u.x, cs, fs});
}
event defaults (i++) {
foreach_face()
muc.x[] = fm.x[]/10000.;
boundary ((scalar*){muc});
}
event tracer_diffusion (i++) {
diffusion (s, dt, muc);
}
event adapt (i++) {
cs.prolongation = refine_injection;
boundary ({cs});
adapt_wavelet ({cs, u}, (double[]){0.001, ue, ue, ue}, maxlevel, 1.);
cs.prolongation = fraction_refine;
boundary ({cs});
}
#define INDEX(i,j) (nx*i + j) // Cartesian index for nx*nx array
event output_flux_plane (t += 10) {
int nx = 1 << maxlevel;
double ustar[nx*nx*3]; //quasi 3D array initialized with high values
for (int i = 0; i < nx; i++)
for (int j = 0; j < nx; j++)
ustar[INDEX(i,j)] = nodata;
foreach() {
if (cs[] > 0. && cs[] < 1. && level == maxlevel) {
coord n = facet_normal (point, cs, fs), p;
double alpha = plane_alpha (cs[], n);
plane_area_center (n, alpha, &p);
int xi = point.i - GHOSTS, zj = point.k - GHOSTS;
double v[dimension];
ustar[INDEX(xi, zj)] = u_tau (point, cs, n ,p, (scalar*){u}, v);
}
}
#if _MPI //We should reduce all copies of ustar to a single one
MPI_Barrier (MPI_COMM_WORLD);
if (pid() == 0)
MPI_Reduce (MPI_IN_PLACE, ustar, nx*nx, MPI_DOUBLE, MPI_MIN, 0,
MPI_COMM_WORLD);
else
MPI_Reduce (ustar, NULL, nx*nx, MPI_DOUBLE, MPI_MIN, 0,
MPI_COMM_WORLD);
#endif
//Write it to some space-delimited asci file
if (pid() == 0) {
char fname[99];
sprintf (fname, "ustar%g", t);
FILE * fp = fopen (fname, "w");
for (int i = 0; i < nx; i++) {
for (int j = 0; j < nx; j++)
fprintf (fp, "%g ", ustar[INDEX(i,j)]);
fputc ('\n', fp);
}
fclose (fp);
}
}
event movie (t += 0.2 ; t < 50) {
scalar l2[];
lambda2 (u, l2);
boundary ({l2});
foreach()
l2[] *= cs[];
view (ty = -0.2, theta = 0.3, phi = 0.4);
draw_vof ("cs", "fs", fc = {0.1,0.7,0.2});
squares ("s", n = {0, 0 ,1}, alpha = -L0/2., min = -0.1, max = 1.1);
cells(n = {0, 0 ,1}, alpha = -L0/2.);
isosurface ("l2", -0.05);
save("mov.mp4");
}
We make profiles of the velocity components as a function of distance to the surface.
#include "../profile6.h"
#include "../frac-dist.h"
event profile_maker (t += 10) {
char file[99];
sprintf (file, "prof%g", t);
vertex scalar phi[];
distance_to_surface (cs, fs, phi = phi);
profiles (rough_list, phi, rf = 0.5,
fname = file, min = 0.05, max = 2.6);
}
set yr [0 : 2.3]
set xr [0 : 1.1]
set key box top left
set ylabel 'Distance to surface'
set xlabel 'Velocity: x-component'
plot 'prof20' u 2:1 w l lw 2, 'prof50' u 2:1 w l lw 2 , \
'prof90' u 2:1 w l lw 2