∇-Nabla: Numerical Analysis BAsed LAnguage

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=fmax((γ-1.0)*πρ*Ξ, εP); }

∀ cells speed_of_sound @ 3.0 { πc=√(γ*π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=√(left_c);
  const ℝ right_w=√(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=√(left_ww2);
      const ℝ right_ww2 =right_c*(1.0+γ6*(p_star-right_p)/right_p);
      const ℝ right_ww=√(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=√(left_w2);
  const ℝ right_w2=right_c*(1.0+γ6*(p_star-right_p)/right_p);
  const ℝ right_w=√(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;
}