∇-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

∇ DDFV

The DDFV method has proved to be efficient for approximating on general meshes a broad class of physical problems involving diffusion operators as well in 2D as in 3D.

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=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;
}

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