∇-Nabla: Numerical Analysis BAsed LAnguage

∇-Nabla: Numerical Analysis BAsed LAnguage

nablaLulesh-512.png nablaCoMD-512.png nablaDDFV-512.png nablaMNLDDFV-512.png nablaSethi-512.png nablaSPH-512.png

∇ Sethi

This ∇ port of Hydro which has been extracted from a real astrophysics code RAMSES, developed in CEA Saclay (France) by Romain Teyssier to study large scale structure and galaxy formation.

Sethi solves the compressible Euler equations of hydrodynamics, relying on a finite volume numerical method using a second order Godunov scheme.

with cartesian;

// ****************************************************************************
// * OPTIONS
// ****************************************************************************
options{
   δx             = 0.001;
   time_end       = 50.0;
   courant_number = 0.8; 
   σ_type         = 1.0;
   nstepmax       = 2;
   testcase       = 0;
   uid_bubble_one = 0;
   uid_bubble_two = 0;
   γ              = 1.4;
   εC             = 1e-10;
   εR             = 1e-10;
   ε              = 1e-8;
   nRiemann       = 10;
};

// ****************************************************************************
// * Déclaration des variables aux mailles
// ****************************************************************************
cells{
   old_ρ,old_u,old_v;
   old_E,old_p,old_c;
   ρ,u,v,E,p,c;
   πρ,inv_πρ,πu,πv,Ξ;
   πp,πc,σρ,σu,σv,σp;
   flux_ρ_left,flux_ρ_right;
   flux_u_left,flux_u_right;
   flux_v_left,flux_v_right;
   flux_p_left,flux_p_right;
   q_ρ_left,q_ρ_right;
   q_u_left,q_u_right;
   q_v_left,q_v_right;
   q_p_left,q_p_right;
   φρ,φu,φv,φE;
   γρ,γu,γv,γp;
};

// ****************************************************************************
// * Déclaration des variables globales
// ****************************************************************************
global{
   dt_δx;
   γ6,εP,εPP;
   ZEROL,ZEROR;
};

// ****************************************************************************
// * Partie d'initialisation ]-∞,-0] 
// ****************************************************************************
print_params @ -10.0{
  info()<< "TESTCASE            =" << testcase;
  info()<< "NSTEPMAX            =" << nstepmax;
  info()<< "COURANT_NUMBER      =" << courant_number;
  info()<< "DX                  =" << δx;
  info()<< "TEND                =" << time_end;
}

ini @ -10.0{
  γ6=+1.0)/(2.0*γ);
  εP=εC*εC/γ;
  εPP=εR*εP;
}

 cells init_hydro @ -9.1{
  old_ρ=ρ=1.0;
  old_u=u=old_v=v=0.0;
  old_E=E=1.0e-5;
  if (testcase==0 && uid==0) old_E=E=1.0/δx/δx;
  if (testcase==1 && uid==uid_bubble_one) old_E=E=½/δx/δx;
  if (testcase==1 && uid==uid_bubble_two) old_E=E=0.25/δx/δx;
  old_p=p=old_c=c=0.0;
}

init_δt @ -4.0{
  δt=ini_δt;
  dt_δx=δt/δx;
  ZEROL  =-100.0/dt_δx;
  ZEROR  = 100.0/dt_δx;
  info() << "[init_δt] δt="<<δt<<", dt/δx="<<dt_δx<<"";
}

// ****************************************************************************
// * Partie de calcul ]+0,+∞[ 
// ****************************************************************************
 cells gatherConservativeVariables @ 0.5{
  ρ=old_ρ;
  if (isX()){
    u=old_u; v=old_v;
  }else{
    v=old_u; u=old_v;
  }
  E=old_E;
}

 cells conservative_to_primitive @ 1.0{
  πρ=max(ρ, εR);
  inv_πρ=1.0/πρ;
  πu=u*inv_πρ;
  πv=v*inv_πρ;
  const  E_kinetic=½*(πu*πu+πv*πv);
  Ξ=E*inv_πρ-E_kinetic;
}

 cells equation_of_state @ 2.0 { πp=max((γ-1.0)*πρ*Ξ, εP); }

 cells speed_of_sound @ 3.0 { πc=sqrt*πp*inv_πρ); }


// ****************************************************************************
// * Compute Slopes @ 4.0
// ****************************************************************************
 slope(const  nbv_m, const  nbv_0, const  nbv_p){
   rtn;
  const  left=σ_type*(nbv_0-nbv_m);
  const  right=σ_type*(nbv_p-nbv_0);
  const  center=½*(left+right)/σ_type;
  const  sign=(center>0.0)?1.0:-1.0;
  const  llftrgt=(left*right)<=0.0;
  const  t1=fmin(fabs(left), fabs(right));
  rtn=sign*fmin((1.0-llftrgt)*t1, fabs(center));
  return rtn;
}
 inner compute_inner_slopes_xyz @ 4.0{
  σρ=slope(πρ[prevCell],πρ,πρ[nextCell]);
  σu=slope(πu[prevCell],πu,πu[nextCell]);
  σv=slope(πv[prevCell],πv,πv[nextCell]);
  σp=slope(πp[prevCell],πp,πp[nextCell]);
}
 outer compute_outer_slopes_xyz @ 4.0{
  if (prevCell.null()){ // Frontière gauche
    σρ=slope(πρ,πρ,πρ[nextCell]);
    σu=slope(-πu[nextCell],-πu,πu);
    σv=slope(πv,πv,πv[nextCell]);
    σp=slope(πp,πp,πp[nextCell]);
    continue;
  }
  if (nextCell.null()){ // Frontière droite
    σρ=slope(πρ[prevCell],πρ,πρ);
    σu=slope(πu,-πu,-πu[prevCell]);
    σv=slope(πv[prevCell],πv,πv);
    σp=slope(πp[prevCell],πp,πp);
    continue;
  }
}

// ****************************************************************************
// * Compute fluxes
// ****************************************************************************
 cells compute_right_flux @ 5.0{
  const  inv_c=1.0/πc;
  const  αm=½*(σp*( inv_πρ*inv_c)-σu)*πρ*inv_c;
  const  αp=½*(σp*( inv_πρ*inv_c)+σu)*πρ*inv_c;
  const  αr=σρ-σp*(inv_c*inv_c);
  const  αv=σv;
  const  sp_m=((πu-πc)>=ZEROR)?1.0:(πu-πc)*dt_δx+1.0;
  const  sp_p=((πu+πc)>=ZEROR)?1.0:(πu+πc)*dt_δx+1.0;
  const  sp_0=(πu >= ZEROR)?1.0:πu*dt_δx+1.0;
  const  a_p=-½*sp_p*αp;
  const  a_m=-½*sp_m*αm;
  const  a_0r=-½*sp_0*αr;
  const  a_0v=-½*sp_0*αv;
  q_ρ_right=φρ_right=πρ+(a_p+a_m+a_0r);
  q_u_right=φu_right=πu+(a_p-a_m)*πc*inv_πρ;
  q_v_right=φv_right=πv+a_0v;
  q_p_right=flux_p_right=πp+(a_p+a_m)*πc*πc;
}
 cells compute_left_flux @ 5.0{
  const  inv_c=1.0/πc;
  const  αm=½*(σp*( inv_πρ*inv_c)-σu)*πρ*inv_c;
  const  αp=½*(σp*( inv_πρ*inv_c)+σu)*πρ*inv_c;
  const  αr=σρ-σp*(inv_c*inv_c);
  const  αv=σv;
  const  sp_m=((πu-πc) <= ZEROL)?-1.0:(πu-πc)*dt_δx-1.0;
  const  sp_p=((πu+πc) <= ZEROL)?-1.0:(πu+πc)*dt_δx-1.0;
  const  sp_0=(πu <= ZEROL)?-1.0:πu*dt_δx-1.0;
  const  a_p=-½*sp_p*αp;
  const  a_m=-½*sp_m*αm;
  const  a_0r=-½*sp_0*αr;
  const  a_0v=-½*sp_0*αv;
  q_ρ_left=φρ_left=πρ+(a_p+a_m+a_0r);
  q_u_left=φu_left=πu+(a_p-a_m)*πc*inv_πρ;
  q_v_left=φv_left=πv+a_0v;
  q_p_left=flux_p_left=πp+(a_p+a_m)*πc*πc;
}

// ****************************************************************************
// * Compute qleftright
// ****************************************************************************
 inner cells inner_qleftright_std(xyz direction){
  q_ρ_left=φρ_left[prevCell];
  q_u_left=φu_left[prevCell];
  q_v_left=φv_left[prevCell];
  q_p_left=flux_p_left[prevCell];
}
 outer cells outer_qleftright(xyz direction){
  if (prevCell.null()){ // Frontière gauche
    q_u_left=-φu_right;
    continue;
  }
  if (nextCell.null()){ // Frontière droite
    q_ρ_left=φρ_left[prevCell];
    q_u_left=φu_left[prevCell];
    q_v_left=φv_left[prevCell];
    q_p_left=flux_p_left[prevCell];
    continue;
  }
}
do_qleftright @ 5.5{
  if (isX()){
    inner_qleftright_std(MD_DirX);
    outer_qleftright(MD_DirX);
  }else{
    inner_qleftright_std(MD_DirY);
    outer_qleftright(MD_DirY);
  }
}

// ****************************************************************************
// * Compute Riemann @ 6.0
// ****************************************************************************
 cells riemann @ 6.0{
  const  in_left_ρ=q_ρ_left;
  const  in_left_u=q_u_left;
  const  in_left_v=q_v_left;
  const  in_left_p=q_p_left;
  const  in_right_ρ=q_ρ_right;
  const  in_right_u=q_u_right;
  const  in_right_v=q_v_right;
  const  in_right_p=q_p_right;
  const  left_ρ=fmax(in_left_ρ, εR);
  const  left_u=in_left_u;
  const  left_v=in_left_v;
  const  left_p=fmax(in_left_p, left_ρ*εP);
  const  left_c=γ*left_p*left_ρ;
  const  right_ρ=fmax(in_right_ρ, εR);
  const  right_u=in_right_u;
  const  right_v=in_right_v;
  const  right_p=fmax(in_right_p, right_ρ*εP);
  const  right_c=γ*right_p*right_ρ;
  const  left_w=sqrt(left_c);
  const  right_w=sqrt(right_c);
   p_star;
  Bool goon=true;
  p_star=fmax((right_w*left_p+left_w*right_p+left_w*right_w*(left_u-right_u))
              /(left_w+right_w), 0.0);
  for(int i=0; i<nRiemann; ++i){
    if(goon){
      const  left_ww2=left_c*(1.0+γ6*(p_star-left_p)/left_p);
      const  left_ww=sqrt(left_ww2);
      const  right_ww2 =right_c*(1.0+γ6*(p_star-right_p)/right_p);
      const  right_ww=sqrt(right_ww2);
      const  tmp_num=2.0*left_ww2*right_ww2*(left_ww*right_ww*(left_u-right_u)
                                              -left_ww*(p_star-right_p)
                                              -right_ww*(p_star-left_p));
      const  tmp_den =right_ww2*right_ww*(left_ww2+left_c)+left_ww2*left_ww*(right_ww2+right_c);
      const  tmp=tmp_num /(tmp_den);
      const  deleft_p=fmax(tmp, -p_star);
      p_star+= deleft_p;
      {
        const  uo=fabs(deleft_p/(p_star+εPP));
        goon=uo>ε;
      }
    }
  }
  const  left_w2=left_c*(1.0+γ6*(p_star-left_p)/left_p);
  const  left_w=sqrt(left_w2);
  const  right_w2=right_c*(1.0+γ6*(p_star-right_p)/right_p);
  const  right_w=sqrt(right_w2);
  const  u_star=½*(left_u+(left_p-p_star)/left_w+right_u-(right_p-p_star)/right_w);
  const  sgnm=(u_star>0.0)?1.0:-1.0;
  const  ρ_0=(u_star>0.0)?left_ρ:right_ρ;
  const  u_0=(u_star>0.0)?left_u:right_u;
  const  p_0=(u_star>0.0)?left_p:right_p;
  const  w_0=(u_star>0.0)?left_w:right_w;
  const  inv_ρ_0=1.0/(ρ_0);
  const  c_0=fmax(√(fabs*p_0*inv_ρ_0)), εC);
  const  ρ_star=fmax(ρ_0/(1.0+ρ_0*(p_0-p_star)/(w_0*w_0)), εR);
  const  c_star=fmax(√(fabs*p_star/ρ_star)), εC);
  const  ushock=w_0*inv_ρ_0-sgnm*u_0;
  const  spout=(p_star>=p_0)?ushock:c_0-sgnm*u_0;
  const  spin=(p_star>=p_0)?ushock:c_star-sgnm*u_star;
   frac;
  if(spout < 0.0) frac=0.0;
  else if(spin>0.0) frac=1.0;
  else {
    const  scr=fmax(spout-spin, εC+fabs(spout+spin));
    frac=fmax(fmin((1.0+(spout+spin)/(scr))*½, 1.0), 0.0);
  }
  γρ=(frac*ρ_star+(1.0-frac)*ρ_0);
  γu=(frac*u_star+(1.0-frac)*u_0);
  γv=(u_star>0.0)?left_v:right_v;
  γp=(frac*p_star+(1.0-frac)*p_0);
}

// ****************************************************************************
// * Compute fluxes from solution @ 7.0
// ****************************************************************************
 cells cmpflx @ 7.0{
  const  mass_density=γρ*γu;
  φρ=mass_density;
  φu=mass_density*γu+γp;
  φv=mass_density*γv;
  const  E_kinetic=½*γρ*(γu*γu+γv*γv);
  const  E_total=γp*1.0/-1.0)+E_kinetic;
  φE=γu*(E_total+γp);
}

// ****************************************************************************
// * Update fluxes from solution @ 8.0
// ****************************************************************************
 cells updateXYZ(xyz direction) {
   dtSδx=δt/δx;
  if (!this->isOwn()) continue;
  if (prevCell.null()){ // Frontière prev
    old_ρ=ρ+(φρ-φρ[nextCell])*dtSδx;
    if (isX()){
      old_u=u+(φu-φu[nextCell])*dtSδx;
      old_v=v+(φv-φv[nextCell])*dtSδx;
    }else{
      old_v=u+(φu-φu[nextCell])*dtSδx;
      old_u=v+(φv-φv[nextCell])*dtSδx;
    }
    old_E=E+(φE-φE[nextCell])*dtSδx;
    continue;
  }
  if (nextCell.null()){ // Frontière next
    old_ρ=ρ+φρ*dtSδx;
    if (isX()){
      old_u=u+φu*dtSδx;
      old_v=v+φv*dtSδx;
    }else{
      old_v=u+φu*dtSδx;
      old_u=v+φv*dtSδx;
    }
    old_E=E+(φE)*dtSδx;
    continue;
  }
  old_ρ=ρ+(φρ-φρ[nextCell])*dtSδx;
  if (isX()){
    old_u=u+(φu-φu[nextCell])*dtSδx;
    old_v=v+(φv-φv[nextCell])*dtSδx;
  }else{
    old_v=u+(φu-φu[nextCell])*dtSδx;
    old_u=v+(φv-φv[nextCell])*dtSδx;
  }
  old_E=E+(φE-φE[nextCell])*dtSδx;
}
update @ 8.0{
  if (isX()) updateXYZ(MD_DirX);
  else updateXYZ(MD_DirY);
}

// ****************************************************************************
// * recompute_δt & test for quit @ [10.0,+∞[ 
// ****************************************************************************
recompute_δt @ 10.0{
  if ((GlobalIteration%4) || (GlobalIteration==1)) return;
  δt=compute_δt();
  dt_δx=δt/δx;
}

isItTimeToQuit @ 11.0{
  if (GlobalIteration >= (nstepmax<<1)) exit;
  if (time < time_end) return;
  exit;
}

camierjs@nabla-lang.org, Jean-Sylvain.Camier@cea.fr