dirichlet_fixed_timestep.c

Solidifying planar interface

Solidifying planar interface

Exact solution. Speed should be constant throughout the simulation. Temperature profile :

\displaystyle T(x,t) = \left\{\begin{array}{cc} -1 + e^{-V(x-Vt)}& , x > Vt\\ 0 & , x \leq Vt\\ \end{array}\right.

Plot of the L_1-error.

f(x) = x
set key left
plot 'log0' u 1:2 w l t '32x32   VOF height',\
     'log1' u 1:2 w l t '64x64   VOF height',\
     'log2' u 1:2 w l t '128x128 VOF height',\
     'log0' u 1:3 w p pt 7 lc 'blue' t '32x32   LS',\
     'log1' u 1:3 every 2 w p pt 7 lc 'green' t '64x64   LS',\
     'log2' u 1:3 every 4 w p pt 7 lc 'red' t '128x128 LS',\
     f(x)+1.e-6

Position of the interface (script)

unset xrange
unset yrange

ftitle(a,b) = sprintf("$%.1fx^\{%4.2f\}$", exp(a), -b)

f(x) = a + b*x
fit f(x) 'log' u (log($1)):(log($2)) via a,b

f2(x) = a2 + b2*x
fit f2(x) 'log' u (log($1)):(log($4)) via a2,b2

f3(x) = a3 + b3*x
fit f3(x) 'log' u (log($1)):(log($6)) via a3,b3

set ylabel 'Average error'
set xrange [16:512]
set yrange [*:*]
set xtics 16,2,512
set format y "%.1e"
set logscale xy
set key bottom left
plot '' u 1:2 pt 7 lc 'blue' t 'all cells', exp(f(log(x))) t ftitle(a,b) lc 'blue', \
     '' u 1:4 pt 7 lc 'green' t 'partial cells', exp(f2(log(x))) t ftitle(a2,b2) lc 'green', \
     '' u 1:6 pt 7 lc 'red' t 'full cells', exp(f3(log(x))) t ftitle(a3,b3) lc 'red'

Average error convergence CFL = 0.5 (script)

set term svg size 1000,1000
set key right
unset logscale
unset xrange
set ylabel 'Temperature profile' 
unset xrange
set xtics 0,0.1,0.5
plot 'out0' u 1:2 w p pt 7 lc 'blue' t 'Final temperature 16x16  ',\
     'out1' u 1:2 w p pt 7 lc 'green' t '32x32  ',\
     'out2' u 1:2 w p pt 7 lc 'orange' t '64x64  ',\
     'out3' u 1:2 w p pt 7 lc 'red' t '128x128', \
     '' u 1:3 w l lt -1 t 'Reference final',\
     'init0' u 1:2 w p pt 7 lc 'blue' t 'Initial temperature 16x16 ',\
     'init1' u 1:2 w p pt 7 lc 'green' t '32x32  ',\
     'init2' u 1:2 w p pt 7 lc 'orange' t '64x64  ',\
     'init3' u 1:2 w p pt 7 lc 'red' t '128x128',\
     '' u 1:3 w l lt -1 t 'Reference initial'

Initial and final temperature plots (script)

set key left
set logscale y
set ylabel 'Error profile'
unset xrange
plot 'out0' u 1:4 w p pt 7 lc 'blue' t '16x16',\
     'out1' u 1:4 w p pt 7 lc 'green' t '32x32',\
     'out2' u 1:4 w p pt 7 lc 'orange' t '64x64',\
     'out3' u 1:4 w p pt 7 lc 'red' t '128x128'

Final error plots (script)

#define DOUBLE_EMBED 1
#define LevelSet     1
#define Gibbs_Thomson 0
#define GHIGO 1
#define QUADRATIC 1

#if GHIGO
#include "../../ghigo/src/myembed.h"
#else
#include "embed.h"
#endif
#include "../embed_extrapolate_3.h"
#include "../double_embed-tree.h"

#include "../advection_A.h" // to remove any calculation of timesteps elsewhere
#if GHIGO
#include "../../ghigo/src/mydiffusion.h"
#else
#include "diffusion.h"
#endif

#include "fractions.h"
#include "curvature.h"

#include "../level_set.h"
#include "../LS_curvature.h"
#include "../LS_advection.h"

#include "view.h"


#define T_eq          0.
double TL_inf(double x,double V, double t){
    return -1+exp(-V*(x-V*t));
}
#define TS_inf        0.


int MINLEVEL, MAXLEVEL; 
double latent_heat;

Setup of the physical parameters + level_set variables

redefinition of ‘TL’
redefinition of ‘TS’
redefinition of ‘dist’

scalar TL[], TS[], dist[];

redefinition of ‘vpc’
redefinition of ‘vpcf’

vector vpc[],vpcf[];

redefinition of ‘curve’

scalar curve[];

redefinition of ‘tracers’

scalar * tracers   = {TL};

redefinition of ‘tracers2’

scalar * tracers2  = {TS};

redefinition of ‘level_set’

scalar * level_set = {dist};

redefinition of ‘muv’

face vector muv[];
scalar grad1[], grad2[];
double DT2;

double  epsK = 0., epsV = 0.;
double eps4 = 0.;



double lambda[2];

#define GT_aniso 0
int aniso = 1;

int     nb_cell_NB =  1 << 3 ;  // number of cells for the NB
double  NB_width ;              // length of the NB

redefinition of ‘s_clean’

double s_clean = 1.e-10; // used for fraction cleaning

TL[embed] = dirichlet(T_eq);
TS[embed] = dirichlet(T_eq);

The inteface moves at constant speed 1.

#define V  1.

‘t0’ undeclared (first use in this function); did you mean ‘t’?

TL[top]    = dirichlet(TL_inf(y,V,t+t0)); 
TS[bottom] = dirichlet(TS_inf); 


int j;
int k_loop = 0;
double t0;

#define geo(x, V,t) (x-V*t)


char filename [100];
FILE * fp1, * fp2, * fp3, * fp4;

void myprop(face vector muv, face vector fs, 
  double lambda){
  foreach_face()
  muv.x[] = lambda*fs.x[];
  boundary((scalar *) {muv});
}


int main() {
  L0 = 1.;
  CFL = 0.5;
  origin (-0.5*L0, -0.5*L0);
  periodic(left);  
  j = 1;
  snprintf(filename, 100,  "init_speed");
  fp4 = fopen (filename,"w");
  for (j=0;j<=3;j++){
    MAXLEVEL = 5+j, MINLEVEL = 5 ;
    DT2 = 1.e-6;
    DT  = DT2;
    fprintf(stderr, "##DT %g\n", DT2);
    snprintf(filename, 100,  "log%d", j);
    fp1 = fopen (filename,"w");
    snprintf(filename, 100,  "out%d", j);
    fp2 = fopen (filename,"w");
    snprintf(filename, 100,  "init%d", j);
    fp3 = fopen (filename,"w");

Here we set up the parameters of our simulation. The latent heat L_H, the initial position of the interface h_0 and the resolution of the grid.

    latent_heat  = 1;
    TL.third = true;
    TS.third = true; 
    N = 1 << MAXLEVEL;
    init_grid (N);
    run(); 
    fclose(fp1);
    fclose(fp2);
    fclose(fp3);
  }
  fclose(fp4);
}

event init(t=0){
  t0 = 1.e-6;
  fprintf(stderr, "##%g\n",t0 );

  NB_width = nb_cell_NB * L0 / (1 << MAXLEVEL);

  lambda[0] = 1.;
  lambda[1] = 1.;
  foreach(){
    dist[] = geo(y,V,t0);
  }
  boundary ({dist});
  restriction({dist});

  vertex scalar distn[];
  foreach_vertex(){
    distn[] = geo(y,V,t0);
  }
  boundary({distn});
  restriction({distn});

  fractions (distn, cs, fs);
  boundary({cs,fs});
  restriction({cs,fs});

  foreach(){
    if(cs[]){
      TL[] = TL_inf(y,V,t0);
    }
    else{
      TL[]=  nodata;
    }
    TS[] = 0.;
    boundary({TS,TL});
  }

  boundary({TL,TS});
  restriction({TL,TS});

  foreach(){
    if(x==Delta/2. && TL[] != nodata){
      fprintf(fp3, "%.8g %.8g %.8g %.8g\n", y, TL[], TL_inf(y,V,t0),
       fabs(TL[]-TL_inf(y,V,t0)));}
  }
  int n =0;
  foreach(reduction(+:n)){
    n++;
  }
  myprop(muv,fs,lambda[0]);

#if GHIGO
  u.n[embed] = dirichlet(0);
  u.t[embed] = dirichlet(0);

  // Added to account for moving boundary algorithm
  
  uf.n[embed] = dirichlet (0.);
  uf.t[embed] = dirichlet (0.);
#endif
  // event("final_error");
}

This event is the core the of the hybrid level-set/embedded boundary.

event velocity(i++){
  double lambda1 = lambda[0], lambda2 = lambda[1];

implicit declaration of function ‘LS_speed’ [-Werror=implicit-function-declaration]

  LS_speed(
    dist,latent_heat,cs,fs,TS,TL,T_eq,
    vpc,vpcf,lambda1,lambda2,

‘deltat’ undeclared (first use in this function); did you mean ‘pfloat’?

    epsK,epsV,eps4,deltat=0.45*L0/(1<<MAXLEVEL),
    itrecons = 60,tolrecons = 1.e-10,NB_width);
  tnext = t+DT2;
  dt = DT2;
  // scalar mystat[];
  // foreach(){
  //   if(cs[])mystat[] = vpc.y[];
  //   else mystat[] = nodata;
  // }
  // stats s2 = statsf(mystat);
  // fprintf(stderr, "%g %g %g\n",t, s2.min, s2.max);
}



event movies ( i++,last)
{
  double y_max=0,y_LS = 0.;
  vector h[];
  heights (cs, h);
  boundary((scalar *){h});
  foreach(reduction(max:y_max) reduction(max:y_LS)){
    if(interfacial(point, cs) && y >0){
      double yy = y+Delta*height(h.y[]);
      if(yy < 1.e10)y_max = max(y_max,yy);
      y_LS = max(y_LS,y-dist[]);
    }     
  }

We’ve already done the advection of the interface, so the time is t+dt.

  fprintf(fp1, "%g %g %g\n", t+dt, y_max,y_LS);
}

event tracer_diffusion(i++,last){
  double lambda1 = lambda[0], lambda2 = lambda[1];
  advection_LS(

missing braces around initializer [-Wmissing-braces]

  dist,
  latent_heat,

incompatible types when initializing type ‘int’ using type ‘vector’ {aka ‘struct ’}

  cs,fs,
  TS,TL,
  T_eq,

incompatible types when initializing type ‘int’ using type ‘vector’ {aka ‘struct ’}

  vpc,vpcf,
  lambda1,lambda2,

excess elements in struct initializer

  epsK,epsV,eps4,

excess elements in struct initializer

  curve,

excess elements in struct initializer

  &k_loop,

‘struct LSadv’ has no member named ‘deltat’

  deltat = 0.45*L0 / (1 << grid->maxdepth),
  itredist = 3,

‘struct LSadv’ has no member named ‘tolredist’; did you mean ‘itredist’?

  tolredist = 3.e-3,

‘struct LSadv’ has no member named ‘itrecons’; did you mean ‘itredist’?

  itrecons = 16,

‘struct LSadv’ has no member named ‘tolrecons’

  tolrecons = 1.e-2,
  s_clean = 1.e-10,
  NB_width);

  boundary({TL});
  myprop(muv,fs,lambda[0]);
  diffusion (TL,dt,muv);
  invertcs(cs,fs);
  myprop(muv,fs,lambda[1]);
  boundary({TS});
  diffusion (TS,dt,muv);
  invertcs(cs,fs);
  myprop(muv,fs,lambda[0]);


  // vector vpcf2[];
  // LS_speed(
  //   dist,latent_heat,cs,fs,TS,TL,T_eq,
  //   vpc,vpcf2,lambda1,lambda2,
  //   epsK,epsV,eps4,deltat=0.45*L0/(1<<MAXLEVEL),
  //   itrecons = 30,tolrecons = 1.e-3,NB_width); 

  // foreach()
  //   foreach_dimension()
  //     vpcf.x[] = 0.5*(vpcf.x[] + vpcf2.x[]);
  // boundary((scalar *){vpcf});



  foreach_face(){
    uf.x[] = 0.;
  }
  boundary((scalar *){uf});
  restriction((scalar *){uf});
}

event final_error(i=400){
  scalar e[], ep[], ef[];
  foreach() {
    if (cs[] == 0.){
      e[] = nodata;
      ep[] = nodata;
      ef[] = nodata;
    }
    else {
      e[] = TL[] - TL_inf( y,V,t+t0);
      ep[] = cs[] < 1. ? e[] : nodata;
      ef[] = cs[] >= 1. ? e[] : nodata;
    }
  }
  foreach(){
    if(x==Delta/2. && cs[] > 0){
      fprintf(fp2, "%.8g %.8g %.8g %.8g %.3g %.3g TEST2\n", y, TL[], TL_inf(y,V,t+t0),
        fabs(e[]),log(fabs(e[])),cs[]);
    }
  }
  norm n = normf (e);
  norm np = normf (ep);
  norm nf = normf (ef);
  fprintf (stderr, "%d %.3g %.3g %.3g %.3g %.3g %.3g %g\n",
   N, n.avg, n.max, np.avg, np.max, nf.avg, nf.max, dt);
  fflush (ferr);


  char filename [100];
  scalar loge[];
  foreach(){
    if(cs[] ==0.)
      loge[] = nodata;
    else
      loge[] = log(fabs(e[]));
  }
  sprintf(filename,"error_grid_%d.png",MAXLEVEL);
  cells();
  squares("loge", min = -11 , max = -5);
  save (filename);
}



#if 1
event adapt (i++, last) {

  foreach_cell(){
    cs2[] = 1.-cs[];
  }
  foreach_face(){
    fs2.x[]      = 1.-fs.x[];
  }

  boundary({cs,cs2,fs,fs2});
  fractions_cleanup(cs,fs,smin = 1.e-14);
  fractions_cleanup(cs2,fs2,smin = 1.e-14);
  restriction({cs,cs2,fs,fs2});
  int n=0;
  boundary({TL,TS});
  scalar visu[];
  foreach(){
    visu[] = (cs[])*TL[]+(1.-cs[])*TS[] ;
  }
  boundary({visu});
  restriction({visu});

  adapt_wavelet ({cs,visu,vpcf.x,vpcf.y},
    (double[]){1.e-2,1.e-6,1.e-6,1.e-6},MAXLEVEL, MINLEVEL);
  foreach(reduction(+:n)){
    n++;
  }
}
#endif

Here we first check that the interface is correctly advected, error on the position of the interface should follow f(x) = x

[chen_simple_1997]

S Chen, B Merriman, S Osher, and P Smereka. A simple level set method for solving Stefan problems. Journal of Computational Physics, 135(1):8–29, 1997.