∇-Nabla: Numerical Analysis BAsed LAnguage

with ℝ²,cartesian;

const ℝ π = acos(-1.0); 
const ℝ D2R = π/180.0; // degrees ⇒ radians
const ℝ EarthMeridionalRadius = 6367449.0;
const ℝ D2M = EarthMeridionalRadius*D2R; // degrees ⇒ meters
const ℝ G = 9.80665; // AccelerationDueToGravity m/s²

ℝ deg_to_r(ℝ θ){ return θ*D2R; }
ℝ lat_to_m(ℝ δy){ return δy*D2M; }
ℝ lon_to_m(ℝ δx, ℝ δy){ return δx*D2M*cos(δy*D2R); }

options{
  ℕ NX                    = 32;      // Number of inner X cells
  ℕ NY                    = 32;      // Number of inner Y cells
  ℝ LENGTH                = 1.0;
  ℕ X_EDGE_ELEMS          = NX+2;    // Inner + Fictitious X cells
  ℕ Y_EDGE_ELEMS          = NY+2;    // Inner + Fictitious Y cells
  ℝ option_λ              = 1.e-4;
  ℝ option_stoptime       = 60.0;
  ℕ option_max_iterations = 32;
  ℝ option_epsd           = 1.0e-9;
  ℝ option_δt             = 0.0025;
};

cells{
  ℝ h,hn,hnp; // height
  ℝ x,dx,inv_dx,y,dy,inv_dy;
  ℝ d,d_hn; // depth
  ℝ un,vn,unp,vnp; // velocities
  ℝ deqh,deqh_dx,deqh_dy;
  ℝ coef_gradx_h,coef_grady_h;
  ℝ fc;       // Coriolis force
};

global{
  ℝ xmin,xmax, ymin,ymax;
  ℝ dmax,dxmax,dymax;
  ℝ dx_lon, dy_lat;
  ℕ chk,chkH,chkU,chkV;
};

// Geometry: X, Y and Z Beach + Island
∀ cells @ -19.99 { x=0.0; ∀ nodes x+=coord.x; x*=¼*option_λ;}
∀ cells @ -19.99 { y=0.0; ∀ nodes y+=coord.y; y*=¼*option_λ;}
∀ cells @ -19.99 {
  ℝ i = x/option_λ;
  ℝ j = y/option_λ;
  ℝ islnd = ½*expf(-32.0*π²*((i-½)²+(j-½)²));
  ℝ beach =(2.0/3.0)*expf(-4.0*π*((i-1.0)²+0.1*(j-1.0)²));
  d -= islnd + beach - ¼;
  }

// Analytical inlet functions
ℝ inlet(const ℝ τ){ return -exp(-((-10.0+τ)²/(4.0*π)))*sin(½*τ)/(12.0*π); }

// Geometry Min Max & Coriolis Force
∀ cells xmin <?= x @ -17;
∀ cells xmax >?= x @ -17;
∀ cells ymin <?= y @ -17;
∀ cells ymax >?= y @ -17;

dxLon @ -15 { dx_lon = (xmax-xmin)/(NX-1);}
dyLat @ -15 { dy_lat = (ymax-ymin)/(NY-1);}

∀ cells @ -13 { dx = lon_to_m(dx_lon,y);}
∀ cells @ -13 { dy = lat_to_m(dy_lat);}

∀ cells dxmax >?= dx @ -11;
∀ cells dymax >?= dy @ -11;
∀ cells @ -11 { inv_dx = 1.0/dx; }
∀ cells @ -11 { inv_dy = 1.0/dy; }
∀ cells @ -11 { coef_gradx_h = option_δt*G/dx; }
∀ cells @ -11 { coef_grady_h = option_δt*G/dy; }

∀ cells space_scheme_init_coriolis @ -11 {
  ℝ T_SIDEREAL = 86164.1;
  ℝ ΩT = 2.0*π/T_SIDEREAL;
  ℝ dΩt = 2.0*ΩT;
  fc = dΩt*sin(deg_to_r(y));
}

// Depth Init
∀ cells dmax >?= d @ -17;
∀ outer west cells @ -13 { d=d[→];}
∀ cells @ -11,+1 { d_hn = d + hn; }

// CFL & δt
ini_δt @ -8 { δt=½*option_δt; }

time_scheme_assert_δt_vs_cfl @ -7 {
  ℝ cgmax = √(G*dmax);
  ℝ idxmx = 1.0/dxmax;
  ℝ idymx = 1.0/dymax;
  ℝ cfl = fmax(δt*cgmax*idxmx, δt*cgmax*idymx);
  assert(δt<cfl);
}

set_δt @ -1 { δt=2.0*δt; }

// Velocity U Init
∀ own inner cells nlsw1_eqU @ -4,+2 {
  ℝ tv = ¼*(vn[↓]+vn[↘]+vn+vn[→]);
  ℝ ul = (un>0.0)?un:un[→];
  ℝ ur = (un>0.0)?un[←]:un;
  ℝ uu = (tv>0.0)?un:un[↑];
  ℝ ud = (tv>0.0)?un[↓]:un;
  ℝ tu1 = un*(ul-ur);
  ℝ tu2 = tv*(uu-ud);
  ℝ thu = G*(hn[→]-hn);
  ℝ tfc = -tv*fc;
  ℝ dequ = (tu1+thu)*inv_dx + tu2*inv_dy + tfc;
  unp = un - δt*dequ;
  }

∀ inner cells ini_update_un @ -3,+3 { un = unp; }

// Velocitu V Init
∀ own inner cells nlsw1_eqV @ -4,+2 {
  ℝ tu = ¼*(un[←]+un+un[↖]+un[↑]);
  ℝ vl = (tu>0)?vn:vn[→];
  ℝ vr = (tu>0)?vn[←]:vn;
  ℝ vu = (vn>0)?vn:vn[↑];
  ℝ vd = (vn>0)?vn[↓]:vn;
  ℝ tv1 = tu*(vl-vr);
  ℝ tv2 = vn*(vu-vd);
  ℝ thv = G*(hn[↑]-hn);
  ℝ tfc = tu*fc;
  ℝ deqv = tv1*inv_dx + (tv2+thv)*inv_dy + tfc;
  vnp = vn - δt*deqv;
  }

∀ inner cells ini_update_vn @ -3,+3 { vn = vnp; }

time_init @ 0 { if (time==0) time=δt; }

// Compute loop ]0.0,+∞[
// Iteration
model_iterate @ 1 if (0==(iteration%32)) {
  printf("\n[ %d ] t=%.5fs",iteration,time);
 }

// Quit test
test4quit @ 1 if (iteration==option_max_iterations) { exit; }


// Height Compute loop
∀ own inner cells @ 1.2 {
  const ℝ dhr=(un>0.0)?d_hn:d_hn[→];
  const ℝ dhl=(un[←]>0.0)?d_hn[←]:d_hn;
  deqh_dx = (un*dhr-un[←]*dhl)*inv_dx;
}

∀ own inner cells @ 1.2 {
  ℝ dhu=(vn>0.0)?d_hn:d_hn[↑];
  ℝ dhd=(vn[↓]>0.0)?d_hn[↓]:d_hn;
  deqh_dy = (vn*dhu-vn[↓]*dhd)*inv_dy;
}

∀ own inner cells @ 1.22 { deqh = deqh_dx + deqh_dy; }
∀ inner cells @ 1.23 { hnp = hn - δt*deqh; }
∀ inner cells @ 1.3 { hn = hnp; }

// HN read INLETs
∀ outer west cells @ 7 { hn=inlet(time); }

// HN Boundaries
∀ outer ~west cells @ 7 { ∀ outer south face hn=hn[↑]; }
∀ outer ~west cells @ 7 { ∀ outer north face hn=hn[↓]; }

// H for output
∀ cells @ 9 { h=hn; }

// Velocity U
∀ own inner cells uRunup @ 2.2 {
  ℝ dr = d[→];
  ℝ hnr = hn[→];
  ℝ dhnr = d_hn[→];
  ℝ ε = option_epsd;
  ℝ coef_grad_h = coef_gradx_h;
  ℾ dhnIe = d_hn < ε;
  ℾ dorh = dhnr<ε or hnr<-d;
  ℾ uadh = unp<0.0 and hnr>hn;
  ℾ dIe = dhnr < ε;
  ℾ hid = hn < -dr;
  ℾ hsr = hn>hnr;
  unp=(dhnIe and dorh)?0.0:(dhnIe and uadh)?unp-coef_grad_h*(hn+dr):unp;
  ℾ ft = !dIe and unp>0.0;
  unp=(dIe and hid)?0.0;
  unp+=(dIe and unp>0.0 and hsr)?coef_grad_h*(hnr+d):0.0;
  unp+=(ft and -d>hnr)?coef_grad_h*(hnr+d):0.0;
  unp-=(ft and -dr>hn)?coef_grad_h*(hn+dr):0.0;
}

// UN Boundaries
∀ outer west cells @ -8.7,7.2 { un = hn * √(G/d); }
∀ outer cells @ 7.2 { ∀ outer south faces un=un[↑]; }
∀ outer cells @ 7.2 { ∀ outer north faces un=un[↓]; }

// Velocity V
∀ own inner cells vRunup @ 2.2 {
  ℝ du = d[↑];
  ℝ hnu = hn[↑];
  ℝ dhnu = d_hn[↑];
  ℝ ε = option_epsd;
  ℝ coef_grad_h = coef_grady_h;
  ℾ dorh = dhnu<ε or hnu<-d;
  ℾ vorh = vnp<0.0 and hnu>hn;
  ℾ dhnuIε = dhnu<ε;
  ℾ hnIdu = hn<-du;
  vnp=(d_hn<ε and dorh)?0.0:(d_hn<ε and vorh)?vnp-coef_grad_h*(hn+du):vnp;
  ℾ ft = (!dhnuIε) and vnp>0.0;
  vnp=(dhnuIε and hnIdu)?0.0;
  vnp+=(dhnuIε and vnp>0.0 and hn>hnu)?coef_grad_h*(hnu+d):0.0;
  vnp+=(ft and -d>hnu)?coef_grad_h*(hnu+d):0.0;
  vnp-=(ft and -du>hn)?coef_grad_h*(hn+du):0.0;
}

// VN Boundaries
∀ outer cells @ 7.2 { ∀ outer west faces vn=vn[→]; }
∀ inner north cells @ 7.3 { vn=0.0; }