∇-Nabla: Numerical Analysis BAsed LAnguage
 |
 |
 |
 |
 |
 |
 |
 |
∇ Pennant
This ∇ port of the PENNANT Mini-App, an unstructured mesh mini-app for advanced architecture research, from Charles R. Ferenbaugh at Los Alamos National Laboratory.
PENNANT is an unstructured mesh physics mini-app designed for advanced architecture research. It contains mesh data structures and a few physics algorithms adapted from the LANL rad-hydro code FLAG, and gives a sample of the typical memory access patterns of FLAG.
Further documentation can be found in the 'doc' directory of the PENNANT distribution.
with ℝ²; options{ // Mesh options ℝ LENGTH = 1.125; ℕ X_EDGE_ELEMS = 3; ℕ Y_EDGE_ELEMS = 3; ℕ Z_EDGE_ELEMS = 1; ℕ option_max_iterations = 0; ℝ option_δt_initial = 1e-7; ℝ option_stoptime = 1.0; // Pennant options ℾ trace = false; ℕ cstop = 64; // simulation stop cycle ℝ tstop = 1.0; // simulation stop time ℝ dtmax = 1.e99; // maximum timestep size ℝ dtinit = 0.0025; // initial timestep size ℝ dtfac = 1.2; // factor limiting timestep growth ℕ dtreport = 10; // frequency for timestep reports ℝ cfl = 0.6; // Courant number, limits timestep ℝ cflv = 0.1; // volume change limit for timestep ℝ rinit = 1.0; // initial density for main mesh ℝ einit = 0.0; // initial energy for main mesh ℾ impacts = false; ℾ subregion = true; ℝ sbxmin = 0.0; ℝ sbxmax = 0.3; ℝ sbymin = 0.0; ℝ sbymax = 0.3; ℝ rinitsub = 1.0; // initial density in subregion ℝ einitsub = 40.222; // initial energy in subregion ℝ uinitradial = 0.0; // initial velocity in radial direction ℝ γ = 5./3.; // coeff. for ideal gas equation ℝ ssmin = 0.1; // minimum sound speed for gas ℝ alpha = 0.5; // alpha coefficient for TTS model ℝ q1 = 0.1; // linear coefficients for Q model ℝ q2 = 1.0; // quadratic coefficients for Q model ℝ ε = 1.e-12; ℝ ι = 1.e-99; }; // **************************************************************************** // * Zone/Sides/Corners Variables // **************************************************************************** cells{ ℝ² zx; // zone center coordinates ℝ² zxp; // zone ctr coords, middle of cycle ℝ zarea; // zone area ℝ zvol; // zone volume ℝ zareap; // zone area, middle of cycle ℝ zvolp; // zone volume, middle of cycle ℝ zvol0; // zone volume, start of cycle ℝ zdl; // zone characteristic length ℝ zm; // zone mass ℝ zr; // zone density ℝ zrp; // zone density, middle of cycle ℝ ze; // zone specific internal energy (energy per unit mass) ℝ zetot; // zone total internal energy ℝ zw; // zone work done in cycle ℝ zwrate; // zone work rate ℝ zp; // zone pressure ℝ zss; // zone sound speed ℝ zdu; // zone velocity difference ℝ z0per; // zone tmp variable in PolyGas calcStateAtHalf ℝ² z0uc; // zone tmp centered velocity // Sides ℝ sarea[faces]; // area ℝ svol[faces]; // volume ℝ sareap[faces]; // area, middle of cycle ℝ svolp[faces]; // volume, middle of cycle ℝ² ssurfp[faces];// surface vector ℝ smf[faces]; // mass fraction ℝ² sfp[faces]; // force from pressure ℝ² sfq[faces]; // force from artificial visc. ℝ² sft[faces]; // force from tts // Corners ℝ² cftot[nodes]; // force, total from all sources ℝ cmaswt[nodes]; // contribution to pmaswt ℝ carea[nodes]; // tmp ℝ cevol[nodes]; // tmp ℝ cdu[nodes]; // tmp ℝ cdiv[nodes]; // tmp ℝ ccos[nodes]; // tmp ℝ² cqe0[nodes]; // tmp ℝ² cqe1[nodes]; // tmp ℝ cw[nodes]; // tmp ℝ crmu[nodes]; // tmp // Utilisé par legion ℕ znump; }; // **************************************************************************** // * Point Variables // **************************************************************************** nodes{ ℝ² px,px0,pxp; // point coordinates, start & middle of cycle ℝ² pu,pu0; // point velocity ℝ² pap,pf; // point acceleration, force ℝ pmaswt; // point mass, weighted by 1/r ℾ has_bcx; // Utilisé par legion ℾ has_bcy; // Utilisé par legion }; // **************************************************************************** // * Edge Variables // **************************************************************************** faces{ ℝ² ex,exp; // edge center coordinates & middle of cycle ℝ elen; // edge length }; // **************************************************************************** // * Global Variables // **************************************************************************** global{ ℕ cycle; // simulation cycle number ℝ dt; // current timestep ℝ dtrec; // maximum timestep for hydro ℝ dvovmax;// }; // **************************************************************************** // * Initialization ]-∞,-0.0[ // **************************************************************************** ∀ nodes iniMeshPx @ -33 { px = coord; } ∀ cells iniMeshZx @ -32 { zx=0; ∀ nodes zx+=px; zx/=#nodes; } ∀ faces iniMeshEx @ -32 { ex=0; ∀ nodes ex+=px; ex/=#nodes; } ∀ cells iniMeshCalcVols @ -31.8 { zvol=0; zarea=0; ∀ nodes{ sarea = ½ * (px[#+1]-px)⤫(zx-px); svol = ⅓* sarea * (px.x+px[#+1].x+zx.x); zarea += sarea; zvol += svol; } } ∀ cells iniMeshCalcSideFracs @ -31.7 { ∀ faces smf = sarea/zarea; } iniGlobals @ -20.0 { time=0; cycle=0; } iniTestCase @ -20 { if (!impacts) return; subregion=false; cfl=0.15; } ∀ cells iniZr @ -14.9 { zr=rinit; } ∀ cells iniZe @ -14.9 { ze=einit; } ∀ cells iniZwrate @ -14.9 { zwrate=0; } ∀ cells iniImpacts @ -14.8 if (impacts) { int n=NABLA_NB_CELLS; int l=√(n); int center=n/2+l/2; int k0=n/4+l/3; int k1=k0; if (uid==center || uid==k0 || uid==k1){ zr = rinitsub; ze = einitsub; } } ∀ cells iniSubRegion @ -14.8 if (subregion) { if (!(zx.x > (sbxmin-ε) && zx.x < (sbxmax+ε) && zx.y > (sbymin-ε) && zx.y < (sbymax+ε))) continue; zr = rinitsub; ze = einitsub; } ∀ cells iniZm @ -14.7 { zm=zr*zvol; } ∀ cells iniZetot @ -14.6 { zetot=ze*zm; } ∀ nodes iniPuStd @ -14.5 if (uinitradial==0.0) { pu=0.0; } iniRstDtHydro @ -14.4 { dtrec = +∞; } // **************************************************************************** // * Compute loop ⊂ ]+0.0,+∞[ // **************************************************************************** incCycle @ 1.0 { cycle += 1; } calcGlobalDt @ 2.0 { ℝ dtlast = dt; dt = dtmax; if (cycle == 1) { dt=(dtinit<dt)?dtinit; } else { ℝ dtrecover = dtfac * dtlast; if (dtrecover<dt) dt = dtrecover; } dt=((tstop-time)<dt)?tstop-time; dt=(dtrec<dt)?dtrec; } ∀ nodes savePx @ 3.0 { px0=px; } ∀ nodes savePu @ 3.0 { pu0=pu; } ∀ cells saveZvol @ 3.0 { zvol0=zvol; } ∀ nodes advPosHalf0 @ 3.1 { pxp = px0+½*pu0*dt; } ∀ cells calcZxp @ 3.2 { zxp=0; ∀ nodes zxp+=pxp; zxp/=#nodes; } ∀ faces calcExp @ 3.2 { exp=0; ∀ nodes exp+=pxp; exp/=#nodes; } ∀ cells calcVols @ 3.4 { zareap=0; zvolp=0; ∀ nodes { sareap = ½*(pxp[#+1]-pxp)⤫(zxp-pxp); svolp = ⅓*sareap*(pxp.x+pxp[#+1].x+zxp.x); zareap+=sareap; zvolp +=svolp; } } ∀ cells calcSurfVecs @ 3.5 { ∀ faces ssurfp = rotateCCW(exp-zxp); } ∀ faces calcEdgeLen @ 3.5 { elen = length(pxp[1]-pxp[0]); } ∀ cells calcCharLen @ 3.6 { zdl = 1.e99; ∀ faces zdl=fmin(zdl,4.0*sareap/elen); } ∀ cells calcZrp @ 4.0 { zrp = zm/zvolp; } ∀ cells calcCrnrMass @ 4.1 { ∀ nodes cmaswt = ½*zrp*zareap*(smf+smf[#-1]); } ∀ cells calcEOS1 @ 5.0 { ℝ γm1 = γ - 1.0; ℝ ss2 = fmax(ssmin*ssmin,1.e-99); ℝ rx = zr; ℝ local_ex = fmax(ze,0.0); ℝ local_px = γm1*rx*local_ex; ℝ prex = γm1*local_ex; ℝ perx = γm1*rx; ℝ csqd = fmax(ss2, prex+perx*local_px/(rx*rx)); zp =local_px; z0per = perx; zss = √(csqd); } ∀ cells calcEOS2 @ 5.1 { ℝ dth = ½ * dt; ℝ zminv = 1.0 / zm; ℝ dv = (zvolp - zvol0) * zminv; ℝ bulk = zr * zss * zss; ℝ denom = 1.0 + ½ * z0per * dv; ℝ src = zwrate * dth * zminv; zp += (z0per * src - zr * bulk * dv) / denom; } ∀ cells pgasCalcForce @ 6.0 { ∀ faces sfp=-zp * ssurfp; } ∀ cells ttsCalcForce @ 6.1 { ∀ faces { ℝ svfacinv = zareap / sareap; ℝ srho = zrp * smf * svfacinv; ℝ sstmp = fmax(zss, ssmin); sstmp = alpha * sstmp * sstmp; ℝ sdp = sstmp * (srho - zrp); ℝ³ sqq = -sdp * ssurfp; sft = sqq; } } ∀ cells qcsZuc @ 6.2 { z0uc=0; ∀ nodes z0uc += pu; z0uc/=#nodes; } ∀ cells qcsCornerDiv @ 6.3 { ∀ nodes { ℝ³ up0 = pu; ℝ³ xp0 = pxp; ℝ³ up1 = ½ * (pu + pu[#+1]); ℝ³ xp1 = exp; ℝ³ up2 = z0uc; ℝ³ xp2 = zxp; ℝ³ up3 = ½ * (pu[#-1] + pu); ℝ³ xp3 = exp[#-1]; ℝ cvolume = ½ * (xp2-xp0)⤫(xp3-xp1); carea = cvolume; ℝ³ v1 = xp3-xp0; ℝ³ v2 = xp1-xp0; ℝ de1 = elen[#-1]; ℝ de2 = elen[#]; ℝ minelen = fmin(de1, de2); ccos = ((minelen<ε)?0.0:4.0*(v1⋅v2)/(de1*de2)); cdiv = ((up2-up0)⤫(xp3-xp1)-(up3-up1)⤫(xp2-xp0))/(2.0*cvolume); ℝ³ dxx1 = ½ * (xp1+xp2-xp0-xp3); ℝ³ dxx2 = ½ * (xp2+xp3-xp0-xp1); ℝ dx1 = length(dxx1); ℝ dx2 = length(dxx2); ℝ³ duav = ¼*(up0+up1+up2+up3); ℝ test1 = fabs((dxx1⋅duav)*dx2); ℝ test2 = fabs((dxx2⋅duav)*dx1); ℝ num = (test1>test2?dx1:dx2); ℝ den = (test1>test2?dx2:dx1); ℝ r = num / den; ℝ evol = fmin(√(4.0*cvolume*r),2.0*minelen); ℝ dv1 = length2(up1+up2-up0-up3); ℝ dv2 = length2(up2+up3-up0-up1); ℝ du = √(fmax(dv1, dv2)); cevol = (cdiv < 0.0 ? evol : 0.); cdu = (cdiv < 0.0 ? du : 0.); } } ∀ cells qcsCrmu @ 6.4 { // Kurapatenko viscous scalar ℝ γp1 = γ + 1.0; ∀ nodes { ℝ ztmp2 = ¼*q2*γp1*cdu; ℝ ztmp1 = q1 * zss; ℝ zkur = ztmp2 + √(ztmp2² + ztmp1²); ℝ rmu = zkur * zrp * cevol; crmu = (cdiv>0.0)?0.0:rmu; } } ∀ cells qcsC0qe @ 6.5 { ∀ nodes { cqe0 = crmu * (pu-pu[#-1])/elen[#-1]; cqe1 = crmu * (pu[#+1]-pu)/elen; } } ∀ cells qcsExtraVars @ 6.6 { ∀ nodes { ℝ csin2 = 1.0-ccos²; cw = ((csin2 < 1.e-4) ? 0.0:carea/csin2); ccos = ((csin2 < 1.e-4) ? 0.0:ccos); } } ∀ cells qcsForcesCorners @ 6.7 { ∀ faces sfq = (cw*(cqe1+ccos*cqe0)+ cw[#+1]*(cqe0[#+1]+ccos[#+1]*cqe1[#+1]))/elen; } ∀ cells qcsSetVelDiff @ 6.8 { ℝ z0tmp=0; ∀ nodes { ℝ³ dx = pxp[#+1]-pxp; ℝ³ du = pu[#+1]-pu; ℝ lenx = elen; ℝ dux = (lenx > 0.0)?fabs((du⋅dx))/lenx:0.0; z0tmp = fmax(z0tmp, dux); } zdu = q1 * zss + 2.0 * q2 * z0tmp; } ∀ cells sumCrnrForce @ 7.0 { ∀ nodes cftot=(sfp+sfq+sft)-(sfp[#-1]+sfq[#-1]+sft[#-1]); } ∀ nodes sumCrnrMasses @ 8.0 { pmaswt=0; ∀ cells pmaswt+=cmaswt; } ∀ nodes sumCrnrForces @ 8.1 { pf=0; ∀ cells pf+=cftot; } ∀ /*outer*/ nodes applyFixedBC_X @ 9.0 { if (px.x!=0.0 && px.x!=LENGTH) continue; ℝ³ vfixx = ℝ³(1,0,0); pu = project(pu, vfixx); pf = project(pf, vfixx); } ∀ /*outer*/ nodes applyFixedBC_Y @ 9.1 { if (px.y!=0.0 && px.y!=LENGTH) continue; ℝ³ vfixy = ℝ³(0,1,0); pu = project(pu, vfixy); pf = project(pf, vfixy); } ∀ nodes calcAccel @ 10.0 { pap = pf/fmax(pmaswt,ι); } ∀ nodes advPosFullPu @ 11.0 { pu = pu0 + pap*dt; } ∀ nodes advPosFullPx @ 11.1 { px = px0 + ½*(pu+pu0)*dt; } ∀ cells updateZx @ 13.0 { zx=0; ∀ nodes zx+=px; zx/=#nodes; } ∀ faces updateEx @ 13.0 { ex=0; ∀ nodes ex+=px; ex/=#nodes; } ∀ cells updateVols @ 13.2 { zvol=zarea=0; ∀ nodes { sarea = ½ * (px[#+1]-px)⤫(zx-px); svol = ⅓* sarea * (px.x+px[#+1].x+zx.x); zarea += sarea; zvol += svol; } } ∀ cells calcWork @ 14.0 { zw=0; ℝ dth = ½ * dt; ∀ nodes { ℝ³ sftot = sfp + sfq; ℝ sd1 = sftot⋅(pu0+pu); ℝ sd2 = -sftot⋅(pu0[#+1]+pu[#+1]); ℝ dwork = -dth*(sd1*pxp.x+sd2*pxp[#+1].x); zetot += dwork; zw += dwork; } } ∀ cells calcWorkRate @ 15.0 { ℝ dvol=zvol-zvol0; zwrate = (zw+zp*dvol)/dt; } ∀ cells calcEnergy @ 16.0 { ze = zetot/(zm+ι); } ∀ cells calcRho @ 16.0 { zr = zm / zvol; } resetDtHydro @ 12.0 { dtrec=dtmax; dvovmax=ι;} ∀ cells calcDtCourant @ 19.0 { ℝ du = fmax(zdu, fmax(zss,ι)); ℝ zdthyd = zdl * cfl / du; ℝ dtnew = (zdthyd < dtmax)?zdthyd:dtmax; if (dtnew < dtrec) dtrec = dtnew; } ∀ cells calcDtVolume @ 20.0 { ℝ zdvov = fabs((zvol-zvol0)/zvol0); dvovmax = (zdvov>dvovmax)?zdvov; } calcDtVolumeTst @ 20.1 { ℝ dtnew = dt * cflv / dvovmax; if (dtnew < dtrec) dtrec = dtnew; } incTime @ 30.0 { time += dt; cout<<"End cycle "<<cycle <<", time = "<<time<<", dt = "<<dt<<endl; } tstExit @ 40.0 { if (cycle < cstop && time < tstop) return; if (X_EDGE_ELEMS==3 && Y_EDGE_ELEMS==3 && Z_EDGE_ELEMS==1) assert(iteration==48); //exit; } // **************************************************************************** // * Mathematical functions // **************************************************************************** ℝ length(ℝ² v){ return √(v.x*v.x + v.y*v.y); } ℝ length2(ℝ² v){ return(v.x*v.x + v.y*v.y); } ℝ² rotateCCW(ℝ² v){ return ℝ²(-v.y, v.x,0); } ℝ² project(ℝ² v, ℝ² u){ return v-(v⋅u)*u; }
