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

∇ Monotone Non-Linear DDFV

These ∇ jobs are taken from the MNLDDFV example: it is a monotone non-linear finite volume method for approximating diffusion operators on general meshes

// *****************************************************************************
// Si on a des quads et que l'on souhaite un Randomly distorted quadrilateral mesh
// *****************************************************************************
 nodes void randomDistortedQuads(void) @ -20.0 if (option_rdq and option_quads) {
  const  α=option_rdq_α;
  const  ρ1=drand48()+drand48()-1.0;
  const  ρ2=drand48()+drand48()-1.0;
  if (coord.x == 0.0 || coord.x == 1.0) continue;
  if (coord.y == 0.0 || coord.y == 1.0) continue;
  coord.x+=α*ρ1*Δl;
  coord.y+=α*ρ2*Δl;
}

// *****************************************************************************
// *
// *****************************************************************************
 nodes void randomDistortedTriangles(void) @ -20.0 if (option_rdq and option_triangles) {
  const  α=option_rdq_α;
  const  ρ1=drand48()+drand48()-1.0;
  const  ρ2=drand48()+drand48()-1.0;
  if (coord.x == 0.0 || coord.x == 1.0) continue;
  if (coord.y == 0.0 || coord.y == 1.0) continue;
  coord.x+=α*ρ1*Δl;
  coord.y+=α*ρ2*Δl;
 }


// *****************************************************************************
// *
// *****************************************************************************
 nodes void stronglyNonConvexQuads(void) @ -20.0 if (option_sncq and option_quads){
  const Integer nid=uid+1;
  const Integer nbNodesPerLine=sqrtl(globalNbNodes);
  const  θ=M_PI*option_sncq_θ;
  const  Δborders=Δl*(1.0+cos(θ))/2.0;
   Δ=1.5*Δl;
  // On saute un noeud sur deux
  if (!(nid%2)) continue;
  // On saute la dernière colonne
  if ((nid%nbNodesPerLine)==0) continue;
  // On saute la première colonne
  if (((nid-1)%nbNodesPerLine)==0) continue;
  // On saute la première et dernière ligne
  if ((((nid-1)/nbNodesPerLine)==0)) continue;
  if ((((nid-1)/nbNodesPerLine)==(nbNodesPerLine-1))) continue;
  // A l'avant dernière colonne, on change les valeures
  if (((nid+1)%nbNodesPerLine)==0) Δ=Δborders;
  // A l'avant dernière ligne, on change les valeures
  if ((((nid-1)/nbNodesPerLine)==(nbNodesPerLine-2))) Δ=Δborders;
  coord.x+=Δ;
  coord.y+=Δ;
}

///////////////////////////////////////////////
// ANIsotropic Diffusion Square withOUT Hole //
// But with option_θ that must be tied to 0, //
// for A and B to be positives.              //
///////////////////////////////////////////////
 exact_solution_anisotropic_without_hole_and_null_θ(ℝ³ p){
  assert(option_𝜕Ω_temperature==0.0);
  return option_𝜕Ω_temperature+sin(M_PI*p.x)*sin(M_PI*p.y);
}
 f_anisotropic_without_hole_and_null_θ(ℝ³ p){
  assert(option_θ==0.0);
  return M_PI*M_PI*((option_k+1.0)*sin(M_PI*p.x)*sin(M_PI*p.y)
    -(option_k-1.0)*cos(M_PI*p.x)*cos(M_PI*p.y)*sin(2.0*option_θ));}

///////////////////////////////////////////////
// ANIsotropic Diffusion Square withOUT Hole //
// Found a function so that p and f are >=0  //
///////////////////////////////////////////////
 exact_solution_anisotropic_without_hole(ℝ³ p){
  const  γ=1.0/atan(½);
  return γ*atan-((p.x-½)²+(p.y-½)²));
}
 f_anisotropic_without_hole(ℝ³ p){
  const  x=p.x;
  const  y=p.y;
  const  γ=1.0/atan(½);
  const  θ=option_θ;
  const  k=option_k;
  const  dnm = (1.0+(x²-x+y²-y)²)²;
  const  num =
    -16.0*γ*(k-1.0)*(x-½)*(y-½)*(-x+x²+(y-1.0)*y)*cos(θ)*sin(θ)
    -8.0*γ*(y-½)²*(x²-x+y*(y-1.0))*(k*cos(θ)²+sin(θ)²)
    -8.0*γ*(x-½)²*(x²-x+y*(y-1.0))*(cos(θ)²+k*sin(θ)²)
    +2.0*γ*(1.0+(x²-x+(y-1.0)*y)²)*(k*cos(θ)²+sin(θ)²)
    +2.0*γ*(1.0+(x²-x+(y-1.0)*y)²)*(cos(θ)²+k*sin(θ)²);
 return num/dnm;
}

////////////////////////////////////////////
// ANIsotropic Diffusion Square WITH Hole //
////////////////////////////////////////////
 g_with_hole(ℝ³ p){
  const  θ_hole = option_𝜕Ω_temperature+2.0;
  const  θ_bord = option_𝜕Ω_temperature;
  if (p.x==0.0 || p.x== 1.0) return θ_bord;
  if (p.y==0.0 || p.y== 1.0) return θ_bord;
  return θ_hole;
}

// ****************************************************************************
// * Fonction 'qui tourne' pour trouver le bon secteur dans le cas convex
// ****************************************************************************
 convexLoopOnThisNodeToFindPositiveDualFace( dbg,
                                                Node nd, ℝ³x3 kappa, ℝ³ νs, ℝ³ d,
                                                ℝ³ *j, ℝ³ *k,  *pj,  *pk,
                                                 *α,  *β,
                                                int *face_uid,
                                                 *face_swap){
  foreach nd face{
    if (fnd->nbCell()==2){
      const Node n0=fnd->node(0);
      const Node n1=fnd->node(1);
      const Cell bC = fnd->backCell();
      const Cell fC = fnd->frontCell();
      const ℝ³ s=½*(coord[n0]+coord[n1]);
      const ℝ³ p=cell_mass_center[bC];
      const ℝ³ q=cell_mass_center[fC];
      const  swap=Sin(q-d,s-d)<0.0;
      *face_swap=swap;
      *j=swap?p:q;
      *k=swap?q:p;
      *pj=swap?cell_θ[bC]:cell_θ[fC];
      *pk=swap?cell_θ[fC]:cell_θ[bC];
    }else{
      const Cell bC = fnd->cell(0);
      const ℝ³ p=cell_mass_center[bC];
      const ℝ³ q=½*(coord[fnd->node(1)]+coord[fnd->node(0)]);
      const ℝ³ s=½*(p+q);
      const  swap=Sin(q-d,s-d)<0.0;
      *face_swap=swap;
      *j=swap?p:q;
      *k=swap?q:p;
      *pj=swap?cell_θ[bC]:g(q);
      *pk=swap?g(q):cell_θ[bC];
    }   
    *α=n(d,*k)⋅(kappa⨂νs);
    *β=n(*j,d)⋅(kappa⨂νs);
    if (!(*α>=0.0 && *β>=0.0)) continue;
    *face_uid=fnd->uniqueId().asInteger();
    return true;
  }
  return false;
}


// Pour les faces internes, tout est bien déterminé
// Pas de swap à prévoir pour le coté dual non plus
 own inner faces void dDual(void) @ 1.0 if (!option_indirect){
  const Integer nid0 = 1+node(0)->uniqueId().asInteger();
  const Integer nid1 = 1+node(1)->uniqueId().asInteger();
  const  dbg=option_debug_dual;
  ℝ³ j,k,l,m;
   pj,pk,pl,pm;
   ad,bd,ae,be;
  int tail_face_uid, head_face_uid;
   tail_face_swap, head_face_swap;
  const ℝ³ d=coord[0];
  const ℝ³ e=coord[1];
  const ℝ³ p=cell_mass_center[backCell];
  const ℝ³ q=cell_mass_center[frontCell];
  const ℝ³ νs=n(q,p);
  const  convexD=geomComputeTriangleAlgebraicArea(d,q,p)>0.0;
  const  convexE=geomComputeTriangleAlgebraicArea(e,p,q)>0.0;
  const  convex=convexD&&convexE;
  if (!convex)
    fatal("Direct and !convex here");
  if (dbg) 
    info()<<"\33[32m[interiorDualFluxesApproximation] Face #"<<uid<<": "<<nid0<<"-"<<nid1<<"\33[m";
  {
    const  okTail =
      convexLoopOnThisNodeToFindPositiveDualFace(dbg,node(0),κ,νs,d,
                                                 &j,&k,&pj,&pk,&ad,&bd,
                                                 &tail_face_uid,
                                                 &tail_face_swap);
    const  okHead = 
      convexLoopOnThisNodeToFindPositiveDualFace(dbg,node(1),κ,-νs,e,
                                                 &l,&m,&pl,&pm,&ae,&be,
                                                 &head_face_uid,
                                                 &head_face_swap);
    const  Ad=geomComputeTriangleArea(d,j,k);
    //geomComputeTriangleArea(d,j,e)+geomComputeTriangleArea(e,k,d);
    const  Ae=geomComputeTriangleArea(e,l,m);
    //geomComputeTriangleArea(e,l,d)+geomComputeTriangleArea(d,m,e);
    const  μsd_num=ae*pl+be*pm;
    const  μse_num=ad*pj+bd*pk;
    const  μs_denum=Ae*(ad*pj+bd*pk)+Ad*(ae*pl+be*pm);
    const  null=(μs_denum==0.0);
    const  μsd=null?½:μsd_num;
    const  μse=null?½:μse_num;
    const  μsd_denum=null?Ad:μs_denum;
    const  μse_denum=null?Ae:μs_denum;

    if(okHead && okTail){
      //info()<<"\t\33[33m ok\33[m";
      interior_dual_c_sd = ½*(ad+bd)*μsd/μsd_denum;
      interior_dual_c_se = ½*(ae+be)*μse/μse_denum;
      interior_dual_c_sl = interior_dual_c_sm = -∞;
      interior_dual_c_sj = interior_dual_c_sk = -∞;
      continue;
    }

    if (okHead){ // Ici, c'est node[0] qui est on_𝜕Ω, μ=0.0
      //info()<<"\t\33[33m !okTail  \33[m";
      interior_dual_c_lm = true;
      interior_dual_face_uid=head_face_uid;
      interior_dual_face_swap=head_face_swap;
      interior_dual_c_se = ½*(ae+be)/Ae;
      interior_dual_c_sd = -∞;
      interior_dual_c_sl = ½*ae/Ae;
      interior_dual_c_sm = ½*be/Ae;
      continue;
    }

    if (okTail){ // Ici, c'est node[1] qui est on_𝜕Ω, μ=1.0
      //info()<<"\t\33[33m !okHead  \33[m";
      interior_dual_c_jk = true;
      interior_dual_face_uid=tail_face_uid;
      interior_dual_face_swap=tail_face_swap;
      interior_dual_c_sd = ½*(ad+bd)/Ad;
      interior_dual_c_se = -∞; 
      interior_dual_c_sj = ½*ad/Ad;
      interior_dual_c_sk = ½*bd/Ad;
      continue;
    }
    fatal("Should not be there!");
  }
}

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