Using Basilisk for Atmospheric Boundary Layer Flows

[Antoon van Hooft]

Hello all,

I have a few questions regarding the usage of Basilisk for a numerical study on atmospheric boundary layer flows.

In this study I would like to investigate stably stratified turbulence in a simplified atmosphere modelled with a cooled (no-slip) bottom. The effect of thermal stratification is modelled according to the Boussinesq-approximation (i.e. changes in density (temperature) only enter the momentum equation in the gravity term) such that the flow remains incompressible. For this I need to solve apart from the NS-equations an advection-diffusion equation for temperature and couple these two. Furtheremore, I want to employ different LES sub-grid-scale (SGS) models to investigate their performance by comparing them to benchmark results. I the future i want to incorporate more realstic features to the model and study their effects in the dynamics (i.e. A realistic modelling of the soil as a reservoir for heat). The questions that arise are:

0) Would you advice me to use Basilisk for this study?

1) In the 'solvers and functions' section of the website I could not find an appropriate solver for scalar advection-diffusion. Can the standard scalar-Advection scheme be easily extended to incorporate diffusion? 

2) I Noticed LES models are not present in Basilisk. Since I have case specific wishes for my SGS-models I want to define them myself. Is this simply done by defining an extra force term 'a' as a function of the resolved flow field? Is this as simple as adding a few lines in the source code calculating the SGS stresses from the resolved fields or am i missing something? 

Gr

Antoon

[Stephane Popinet]

Dear Anton,

> 0) Would you advice me to use Basilisk for this study?

Why not!

> 1) In the 'solvers and functions' section of the website I could not
> find an appropriate solver for scalar advection-diffusion. Can the
> standard scalar-Advection scheme be easily extended to incorporate
> diffusion?

You can combine the advection solver with the reaction-diffusion solver:

http://basilisk.fr/src/diffusion.h

> 2) I Noticed LES models are not present in Basilisk. Since I have case
> specific wishes for my SGS-models I want to define them myself. Is this

> simply done by defining an extra force term *'a'* as a function of the

> resolved flow field? Is this as simple as adding a few lines in the
> source code calculating the SGS stresses from the resolved fields or am
> i missing something?

There are various kinds of LES models which can be classified according
to the "order" of the closure. The simplest models are "zero-equation"
models such as Smagorinsky eddy viscosity which can indeed be computed
just as a local function of the local large-scale stresses. This could
indeed be done quite simply by defining the viscosity (e.g. by
overloading the "properties" event of the Navier-Stokes "centered.h"
solver), using the Smagorinsky formula.

Other models (i.e. k-epsilon) can involve additional advection/diffusion
equations, but this can be easily implemented using a combination of the
existing solvers.

cheers

Stephane

[Antoon van Hooft]

Dear Stephane,

Thank you for your answer. It has motivated me to start using Basilisk.

Now I have a more basic question: How can I add a volume force (e.g. gravity)?
In the example of Von Karman vortex street with tracer field i like to ad a force as if gravity is acting on the tracer-field f.

I tried the following:

1. event gravity (i++) {
2.    foreach()
3.    a.y[] = -f[]; 
4. }

but this returned an error on compiling. error: assignment of member ‘y’ in read-only object   
Is the acceleration term 'a' read-only?

gr

Antoon

[jmf]

Dear Antoon, 

tale a look to the rising.c file

http://basilisk.fr/src/test/rising.c

where the updtate of the a vector is explained 

event acceleration (i++) {
  face vector av = a;
  foreach_face(x)
    av.x[] -= 0.98;
}

vector a is a read only objet in Basilisk. 

Regards

josé

--
You received this message because you are subscribed to the Google Groups "basilisk-fr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to basilisk-fr...@googlegroups.com.
To post to this group, send email to basil...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/basilisk-fr/3ac0c2ea-d6d1-4bd8-9fed-b95c33e3ff65%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Antoon van Hooft]

Hello José

Thank you for the reference. Apparently i am still missing something

i tried adding:
...........

event acceleration (i++) {
  face  vector av = a;

  foreach_face(y)
    av.y[] -= f[];
}
..........
to the karman.c example file. However it returns "Floating point exception (core dumped)"  after the initialization step.

Do you know what is wrong? I can image i am making a simple mistake but i am new to this type of simulations



Op dinsdag 20 oktober 2015 23:10:20 UTC+2 schreef jmfullana:

[Andrés Castillo Castellanos]

Dear Antoon, greetings

I think you can find an example of what you are trying to here : http://basilisk.fr/sandbox/acastillo/rb.c

const face vector delta_ij[] = {0.,1.};
face vector a = new face vector;
foreach_face()
  a.x[] = T[]*delta_ij.x[];

or silmply

face vector a = new face vector;
foreach_face(x)
  a.x[] = T[];

Best regards,

[Antoon van Hooft]

Hello Andrés

Thank you for your answer. I copy/pasted your Rayleigh-Bernard code. (without the Auxiliaries) but it still does not work (i.e. the flow remains stationary). I have the feeling the tracer field T[] is handled correctly as it diffuses a bit in the course of the simulation. In order to stimulate the instability I increased the Rayleigh number and grid-resolution and added a sinusoidal perturbation on the temperature-field.
   
How does Basilisk know that this new face factor a should be implemented as the acceleration term "a" in navierstokes/centered.h code? Is it because it is called "a" of because it is defined in the event "acceleration" or are both required ?    

There are now two events acceleration. One in the rb.c code and one in the navierstokes/centered.h. is this the problem?

[Andrés Castillo Castellanos]

Try to compile the code and add the compiler flag -events, this is useful to understand in which order the different parts of the solver are invoked.
In general, when you define an event with the name acceleration, and this name is already used by some other file, Basilisk will join them together
and execute them in the order they were included into the file.
 
Greetings,

[Antoon van Hooft]

Hello, I have not been able to solve the problem. When including the "-events" tag I get (typically) something like this:

events (i = 188, t = 15)
  logfile                   karman2.c:55
188 15 1 1
  movies                    karman2.c:58
  adapt                     karman2.c:81
  set_dtmax                 src/navier-stokes/centered.h:158
  stability                 src/navier-stokes/centered.h:160
  vof                       src/navier-stokes/centered.h:170
  tracer_advection          src/tracer.h:38
  tracer_advection          src/navier-stokes/centered.h:171
  properties                src/navier-stokes/centered.h:178
  advection_term            src/navier-stokes/centered.h:240
  viscous_term              src/navier-stokes/centered.h:270
  acceleration              karman2.c:17
  acceleration              src/navier-stokes/centered.h:305
  projection                src/navier-stokes

Apperently both acceleration terms in my karman2.c file and the centered.h are executed. This seems ok right, or is it in the wrong order ? The karman2.c file should produce a gif file of a gravity current.

1. #include "navier-stokes/centered.h"
2. #include "tracer.h"
3. scalar f[];
4. scalar * tracers = {f};
5. int main() {
6.   DT = 0.1;
7.   L0 = 8.;
8.   origin (-0.5, -L0/2.);
9.   N = 512;
10.   run();
11. }
12. event acceleration (i++)
13. {

1.   const face vector delta_ij[] = {0.,1.};
2.   face vector a = new face vector;
3.   foreach_face()

1.     a.x[] = 100*f[]*delta_ij.x[];
2. }
3. event init (t = 0) {
4.      const face vector muc[] = {0.00078125,0.00078125};
5.      mu = muc;
6.      foreach() {
7.        f[] = (x>3)+0.01*noise();
8.   }
9.       boundary (all);
10. }
11. event logfile (i++)
12.   fprintf (stderr, "%d %g %d %d\n", i, t, mgp.i, mgu.i);
13. event movies (t += 0.5; t <= 15.) {
14.   static FILE * fp = popen ("ppm2gif > vort.gif", "w");
15.   scalar vorticity[];
16.   foreach()
17.     vorticity[] = (u.x[0,1] - u.x[0,-1] - u.y[1,0] + u.y[-1,0])/(2.*Delta);
18.   boundary ({vorticity});
19.   output_ppm (vorticity, fp, box = {{-0.5,-0.5},{7.5,0.5}},
20.           min = -10, max = 10, linear = true);
21.   static FILE * fp1 = popen ("ppm2gif > fg.gif", "w");
22.   output_ppm (f, fp1, box = {{-0.5,-0.5},{7.5,0.5}},
23.           linear = true, min = 0, max = 1);
24. }
25. event adapt (i++) {
26.   adapt_wavelet ({u,f}, (double[]){3e-2,3e-2,3e-2}, 9, 4);
27. }

[Stephane Popinet]

Dear Anton,

There are several things which one need to know before implementing new
acceleration terms. Unfortunately they are not documented yet. I will
try to update the doc when things have stabilised. In the meantime, here
is a brief summary:

1) one can define several events with the same name. This allows to
"overload" events (such as "acceleration") in the various solvers. This
works like this: when several events with the same name appear in a
code, they are executed in their inverse order of appearance. i.e. in
your case, you "acceleration" event will be executed before that in
navier-stokes/centered.h.

2) fields can be "constant" i.e. in your case the acceleration of
gravity is constant (in space). One defines a constant field like this
for example:

event init() {
const face vector g[] = {0, -9.81};
a = g;
}

By default the acceleration 'a' in centered/navier-stokes.h is a
constant field with value zero. Here we replace it with a constant
(vector) field with components (0,-9.81).

Note that we dont't need to overload 'acceleration()' in this case.

cheers

Stephane

PS: the syntax for constant fields is not fully stabilised and may
change in the near future.

[jose.lopez...@gmail.com]

Hi Stephane,

1) one can define several events with the same name. This allows to
"overload" events (such as "acceleration") in the various solvers. This
works like this: when several events with the same name appear in a
code, they are executed in their inverse order of appearance. i.e. in
your case, you "acceleration" event will be executed before that in
navier-stokes/centered.h.

 
What happen in case you intend to add an additional transport mechanism  of a tracer as it is in EHD problems?  In a EHD problem
the volumetric charge is in the most general case advected, diffused and electrically conducted, 

The advection and the diffusion are easily accounted for using "tracer.h". i.e.

scalar q[];  /*the volumetric charge */

#include "navier-stokes/centered.h"
#include "tracer.h"

scalar * tracers = {q}; /* Advection of q is taken into account */

event tracer_diffusion (i++) {

....   /* overload procedure: write your own tracer_diffusion event that will substitute that of "tracer.h" */
} 

I wrote "ehd.h" which is essentially an event named

event ohmic_conduction (i++)

what performs the electric ohmic conduction of q. 

Shall I change the name  of this event from ohmic_conduction (i++) to tracer_diffusion ()?. 
So doing it, a general EHD problem would write:

scalar q[];  /*the volumetric charge */

#include "navier-stokes/centered.h"
#include "tracer.h"
#include "ehd.h" /* with this the ohmic conduction is calculated (which is computed the last given 
                    the name of the event is the same, i.e tracer_diffusion()*/

scalar * tracers = {q}; /* Advection of q is taken into account */
event tracer_diffusion (i++) {

....   /* we write it to add the diffusion of $q$ as in the case of 'acceleration' 
          with the rule of equal name this tracer_diffusion() will be computed penultimate*/
} 

Thanks in advance

Cheers,

Jose.

[Stephane Popinet]

Hi Jose,

> *Shall I change the name of this event from*||*|ohmic_conduction (i++)| totracer_diffusion ()?*.

I think you could. It is really a matter of thinking about how it fits
within the overall time integration scheme. Another way to see this
"overloading" mechanism is as a way to "plug" a particular term at a
specific point within the time integration scheme.

cheers

Stephane

[jose.lopez...@gmail.com]

I understand, Stephane.

This overloading mechanism allows the users a precise "plug" insertion with the minimal cost of a (not so significant as it could be) name of the event.

This reduce a lot the repetition of code!

Cheers

José