∇-Nabla: Numerical Analysis BAsed LAnguage

with ℵ;

// ****************************************************************************
// * Déclaration des options
// ****************************************************************************
options{
  ℾ option_quads             = false;
  ℾ option_triangles         = false;  
  ℝ option_deltat            = 0.001;
  ℝ option_deltat_factor     = 10.0;
  ℝ option_theta_ini         = 300.0;
  ℝ option_theta_hot         = 700.0;
  ℕ option_max_iterations    = 8;
  ℝ alephEpsilon             = 1.e-10;
  ℕ alephUnderlyingSolver    = 0;
  ℕ alephMaxIterations       = 1024;
  ℕ alephPreconditionerMethod = 2;
  ℕ alephSolverMethod        = 0;
  ℕ alephNumberOfCores       = 0;
  ℾ option_aleph_dump_matrix = false;
};

// ****************************************************************************
// * Déclaration des variables aux mailles
// ****************************************************************************
cells{
  ℝ cell_θ;
  ℝ cell_area;
  ℕ cell_nb_nodes;
  ℝ cell_mass_center;
};

// ****************************************************************************
// * Déclaration des variables aux noeuds
// ****************************************************************************
nodes{
  ℝ node_θ;
  ℝ node_area;
  ℾ node_is_an_edge;
};

// ****************************************************************************
// * Déclaration des variables aux faces
// ****************************************************************************
faces{
  ℝ face_θ;
  ℝ sdivs;
  ℝ α,β,δ,γ,σ;
  ℝ Cosθ,Sinθ;
  ℕ cellLowId;
  ℕ nodeLowId;
};

// ****************************************************************************
// * Déclaration des variables globales
// ****************************************************************************
global{
  ℾ full_quads;
  ℾ full_triangles;
};

// ****************************************************************************
// * Partie d'initialisation ]-∞,-0[
// ****************************************************************************
∀ cells iniCells @ -10.0{
  cell_nb_nodes=nbNode; 
  cell_θ=option_theta_ini;
}

∀ nodes iniNodes @ -10.0{
  node_θ=option_theta_ini;
}

// ****************************************************************************
// * Calcul des Aires
// ****************************************************************************
#define computeTriangleArea2D(a, b, c) (½*((b-a)⤫(c-a)))

∀ cells geomComputeTriangleSurface  @ -10.0 if (option_triangles){
  cell_area = computeTriangleArea2D(coord[0],coord[1],coord[2]);
}

∀ cells geomComputeQuadSurface @ -10.0 if (option_quads){
  const ℝ³ fst_edge = coord[2]-coord[0];
  const ℝ³ snd_edge = coord[3]-coord[1];
  cell_area=½*(fst_edge⤫snd_edge);
}

∀ nodes nodeArea @ -9.0{
  node_area=0.0;
  ∀ cell node_area+=cell_area/cell_nb_nodes;
}

// ****************************************************************************
// * Calcul des Centres de Gravité
// ****************************************************************************
∀ cells geomComputeCellTriangleMassCenter @ -10.0 if (option_triangles){
  cell_mass_center=0.0;
  ∀ node cell_mass_center+=coord;
  cell_mass_center/=nbNode;
}

∀ cells geomComputeCellQuadMassCenter @ -10.0 if (option_quads){
  const ℝ³ s0 = coord[0];
  const ℝ³ s1 = coord[1];
  const ℝ³ s2 = coord[2];
  const ℝ³ s3 = coord[3];
  const ℝ³ c = ¼*(s0+s1+s2+s3);
  const ℝ surface[4]= {
    computeTriangleArea2D(s0,s1,c), computeTriangleArea2D(s1,s2,c),
    computeTriangleArea2D(s2,s3,c), computeTriangleArea2D(s3,s0,c)};
  const ℝ³ barycentre_triangle[4] =
    {⅓*(s0+s1+c), ⅓*(s1+s2+c), ⅓*(s2+s3+c), ⅓*(s3+s0+c)};
  ℝ total_surface = 0.0;
  ℝ³ center = ℝ³(0.0,0.0,0.0);
  ∀ node{
    center += barycentre_triangle[#]*surface[#];
    total_surface += surface[#]; 
  }
  cell_mass_center = center / total_surface;
}

// ****************************************************************************
// * Calculs géométriques des σ et sdivs
// ****************************************************************************
∀ faces faceSbySigma @ -9.0{
  ℝ³ fn0 = coord[0];
  ℝ³ fn1 = coord[1];
  ℝ s=norm(fn0-fn1);
  ℝ³ ccc0;
  ℝ³ ccc1;
  if (nbCell==2){
    ccc0=cell_mass_center[frontCell];
    ccc1=cell_mass_center[backCell];
  }
  if (nbCell==1){
    ccc0=cell_mass_center[0];
    ccc1=(fn0+fn1)/2.;
  }
  σ=norm(ccc0-ccc1);
  if (σ==0.0) fatal("faceSbySigma", "σ==0.0");
  sdivs=s/σ;
}

// ****************************************************************************
// * Initialisation du δt en fonction des options
// ****************************************************************************
∀ faces ℝ computeMin_σ(ℝ min_σ) in (face σ){
  min_σ = min(min_σ,σ);
  return min_σ;
}
ini_δt @ -8.5{
  if (option_deltat==0.0){
    ℝ min_σ=computeMin_σ(+∞);
    ℝ reduced_min_σ = mpi_reduce(ReduceMin,min_σ);
    δt=reduced_min_σ/option_deltat_factor;
  }else δt=option_deltat;
  info()<< "\33[7m[ini_δt] δt="<<δt<<"\33[m";
}

// ****************************************************************************
// * Calculs géométriques des Cosθ,Sinθ
// ****************************************************************************
∀ faces faceTheta @ -8.0{
  ℝ³ n0=coord[0];
  ℝ³ n1=coord[1];  
  ℝ³ c1,c0,τ,tan,nrm;
  ℝ sinTauTan,nrm_xτ,nΤ,nTan;
  if (nbCell==1){
    c0=cell_mass_center[0];
    c1=(n0+n1)/2.;
    cellLowId=cell(0).uniqueId();
  }else{
    if (backCellUid>frontCellUid){
      c0=cell_mass_center[frontCell];
      c1=cell_mass_center[backCell];
      cellLowId=frontCell.uniqueId();
   }else{
      c0=cell_mass_center[backCell];
      c1=cell_mass_center[frontCell];
      cellLowId=backCell.uniqueId();
    }
  }
  τ=(c1-c0);      // Vecteur tangentiel des mailles
  nΤ=norm(τ);
  tan=(n1-n0);    // Vecteur tangentiel des noeuds
  nTan=norm(tan);
  // Le sinus de ces deux vecteurs afin d'orienter la normale
  sinTauTan= (τ⤫tan)/(nΤ*nTan);
  nodeLowId=((sinTauTan>0.)?node(0).uniqueId():node(1).uniqueId());
  nrm.x=-tan.y; nrm.y=+tan.x; nrm.z=0.;
  nrm*=(sinTauTan>0)?-1.:+1.;
  nrm_xτ=norm(nrm)*nΤ;
  Cosθ=dot(nrm,τ)/nrm_xτ;
  Sinθ=(nrm⤫τ)/nrm_xτ;
}

// ****************************************************************************
// * updateAlphaBetaGammaDelta
// ****************************************************************************
∀ faces updateAlphaBetaGammaDelta @ -7.0{
  α = δt*sdivs/Cosθ;
  β = γ = δt*(-Sinθ/Cosθ);
  δ = δt*(1./(sdivs*Cosθ));
}

// ****************************************************************************
// * iniNodeIsAnEdge, iniOuterNodeIsAnEdge
// ****************************************************************************
∀ nodes @ -6.0 { node_is_an_edge=false; }
∀ own outer nodes @ -5.0 { node_is_an_edge=true; }

// ****************************************************************************
// * Partie de calcul ]+0,+∞[ 
// ****************************************************************************

rhsInit @ 1.0 {
  alephInitialize();
  ℵ rhs reset;
  ℵ lhs reset;
}

∀ own cells @ 1.1 { ℵ rhs setValue(cell_θ,this, cell_θ); }

∀ own nodes @ 1.2 { ℵ rhs setValue(node_θ,this, node_θ); }

∀ own outer nodes @ 1.3{ ℵ rhs setValue(node_θ,this, option_theta_hot); }

∀ own outer faces @ 1.3{ ℵ rhs setValue(face_θ,this, option_theta_hot); }

// ****************************************************************************
// * alphaCells
// ****************************************************************************
∀ own cells alphaCells @ 3.1{
  ℝ αc, Σα=0.0;
  ∀ face {
    Σα+=sdivs/Cosθ;
    αc=α/cell_area;
    if (nbCell==1){
      ℵ matrix addValue(cell_θ,this, face_θ,f, -αc);     
      continue;
    }
    if (frontCell==*this)
      ℵ matrix addValue(cell_θ,this, cell_θ,backCell, -αc);
    else
      ℵ matrix addValue(cell_θ,this, cell_θ,frontCell,-αc);
  }
  Σα*=(δt/cell_area);
  ℵ matrix addValue(cell_θ,this, cell_θ,this, 1.0+Σα);
}

// ****************************************************************************
// * betaCells
// ****************************************************************************
∀ own cells betaCells @ 3.2{
  Cell cP;
  Node nD,nE;
  ∀ face{
    ℝ βc=β/cell_area;
    nD=(node(0).uniqueId()==nodeLowId)?node(0):node(1);
    nE=(node(0).uniqueId()==nodeLowId)?node(1):node(0);
    cP=(cell(0).uniqueId()==cellLowId)?cell(0):(nbCell==2)?cell(1):cell(0);
    if (*this!=cP) βc*=-1.0;
    ℵ matrix addValue(cell_θ,this, node_θ,nD, βc);
    ℵ matrix addValue(cell_θ,this, node_θ,nE, -βc);
  }
}

// ****************************************************************************
// * Couplage Node-Cell: γ
// ****************************************************************************
∀ own nodes gammaNodes @ 3.3{
  Node nD;
  Cell cP,cQ;
  if (node_is_an_edge) continue;
  ∀ face {
    ℝ γn=γ/node_area;
    nD=(node(0).uniqueId()==nodeLowId)?node(0):node(1);
    cP=(cell(0).uniqueId()==cellLowId)?cell(0):(nbCell==2)?cell(1):cell(0);
    cQ=(cell(0).uniqueId()==cellLowId)?(nbCell==2)?cell(1):cell(0):cell(0);
    if (cP==cQ) fatal("Gamma", "Should have been filtered");
    if (*this!=nD) γn*=-1.0;
    ℵ matrix addValue(node_θ,this, cell_θ,cP, γn);
    ℵ matrix addValue(node_θ,this, cell_θ,cQ, -γn);
  }
}

// ****************************************************************************
// * Node-Node: δ
// ****************************************************************************
∀ own nodes deltaNodes @ 3.4{
  ℝ δn,Σδ=0.0;
  if (node_is_an_edge) continue;
  ∀ face {
    Node other_node = (node(0)==*this)?node(1):node(0);
    δn=δ/node_area;
    Σδ+=1.0/(Cosθ*sdivs);
    ℵ matrix addValue(node_θ,this, node_θ,other_node, -δn);
  }
  Σδ*=δt/node_area;
  ℵ matrix addValue(node_θ,this, node_θ,this, 1.0+Σδ);
}

// ****************************************************************************
// * Conditions de Dirichlet
// ****************************************************************************
∀ own outer faces @ 3.5{ ℵ matrix addValue(face_θ,this, face_θ,this, +1.0); }
∀ own outer nodes @ 3.6{ ℵ matrix addValue(node_θ,this, node_θ,this, +1.0); }

// ****************************************************************************
// * Aleph Solve
// ****************************************************************************
assembleAndSolve @ 4.0 { ℵ solve; }

// ****************************************************************************
// * Récupération des résultats
// ****************************************************************************
∀ own cells @ 4.1 { cell_θ=ℵ lhs getValue(cell_θ,this); }
∀ own nodes @ 4.2 { node_θ=ℵ lhs getValue(node_θ,this); }

// ****************************************************************************
// * Sync & test for Quit
// ****************************************************************************
testForQuit @ 5.0{
  synchronize(cell_θ);
  synchronize(node_θ);
  if (GlobalIteration >= option_max_iterations) exit;
}