sandbox/alimare/LS_speed.h

#include "simple_discretization.h"

#include "LS_recons.h"
#include "phase_change_velocity.h"
#if CURVE_LS
 #include "LS_curvature.h"
#else
#include "curvature.h"
#endif

struct LSspeed{
  scalar dist;
  double L_H;
  scalar cs;
  face vector fs;
  scalar TS;
  scalar TL;
  double T_eq;
  vector vpc;
  vector vpcf;
  double lambda1;
  double lambda2;
  double epsK;
  double epsV;
  double eps4;
  double deltat;
  int itrecons;
  double tolrecons;
  double NB_width;
};

//void LS_speed(struct LSspeed p){
  void LS_speed(scalar dist, double L_H, scalar cs, face vector fs, scalar TS, scalar TL,
                double T_eq, vector vpc, vector vpcf, double lambda1, double lambda2, double epsK,
                double epsV, double eps4, double deltat, int itrecons, double tolrecons, double NB_width){

  // scalar dist      = p.dist;
  // double L_H       = p.L_H;
  // scalar cs        = p.cs;
  // face vector fs   = p.fs;
  // scalar TS        = p.TS;
  // scalar TL        = p.TL;
  // double T_eq      = p.T_eq;
  // vector vpc       = p.vpc;
  // vector vpcf      = p.vpcf;
  // double lambda[2] = {p.lambda1, p.lambda2};
  // double epsK      = p.epsK;
  // double epsV      = p.epsV;
  // double eps4      = p.eps4;
  // double deltat    = p.deltat;
  // int    itrecons  = p.itrecons;
  // double tolrecons = p.tolrecons;
  // double NB_width  = p.NB_width;

  scalar curve[];
#if LS_perf
double start; 
double end; 
#endif
  if(NB_width==0){fprintf(stderr, "Erreur NB_width LS_speed\n" );exit(1);}

#if CURVE_LS
  curvature_LS(dist, curve);
#else
  curvature (cs,curve);
  boundary({curve});
#endif

#if LS_perf
start = omp_get_wtime(); 
#endif 
  phase_change_velocity_LS_embed (cs, fs ,TL, TS, T_eq, vpc, L_H, 
    lambda,epsK, epsV, eps4, curve);
#if LS_perf
end = omp_get_wtime(); 
fprintf(stderr,"phase change velo %f seconds\n", end - start);
#endif

  vector vpcr[];
  foreach(){
    foreach_dimension(){
      if(interfacial(point,cs))vpcr.x[] = vpc.x[];
      else vpcr.x[] = 0.;
    }
  }
  boundary((scalar * ){vpcr});
  restriction((scalar * ){vpcr});
#if dimension == 1
  scalar * speed_recons  = {vpcr.x};
#elif dimension == 2
  scalar * speed_recons  = {vpcr.x,vpcr.y};
#else
  scalar * speed_recons  = {vpcr.x,vpcr.y,vpcr.z};
#endif

  double err = 0;



#if LS_perf
start = omp_get_wtime(); 
#endif

  recons_speed(dist, deltat, speed_recons,
   tolrecons, &err, 
   itrecons, 
   cs, fs,
   NB_width);
#if LS_perf
end = omp_get_wtime(); 
fprintf(stderr,"recons_speed %f seconds\n", end - start);
#endif 
 
  foreach()
    foreach_dimension(){
      if(fabs(dist[])< 0.99*NB_width)
        vpcf.x[] = vpcr.x[];
      else
        vpcf.x[] = 0.;
    }

  boundary((scalar *){vpcf});
  restriction((scalar *){vpcf});
}

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

  LS_speed(
    dist, latent_heat, cs, fs, TS, TL, T_eq,
    vpc, vpcf, lambda1, lambda2,
    epsK, epsV, eps4,
    DT_LS,
    itrecons,
    tolrecons,
    NB_width
  );
  

#if dtLS // only diffusion
  double DT3 = timestep_LS(vpcf,DT,dist,NB_width);
  tnext = t+DT3;
  dt = DT3;
#else // navier stokes
  //dtmax = timestep_LS(vpcf,DT,dist,NB_width);
 DT = timestep_LS(vpcf,DT,dist,NB_width);
#endif
}