/******************************************************************************* * * TVD solver for 2D shock-tube problem * * Alexey Kudriavtsev, Sven-Bodo Scholz, 2008 * ******************************************************************************/ import Array:all; import fluid_2d:all; use Numerical: all except { e}; use StdIO: all; #if defined (DISLIN) use Dislin : all; #define MINS [0d,0d] /* part of the area to be plotted */ #define MAXS [XL,YL] #define STEPS [XL/10d, YL/10d] /* markings on the axes */ objdef DislinCanvas display = createDisplay( ); #if defined (SAVE) #define NSAVE 10 #else /* not defined (SAVE) */ #define NSAVE 10 #endif /* SAVE */ #else /* not defined (DISLIN) */ #define NSAVE 1 #endif /* DISLIN */ #ifndef OUTFILE_TECPLOT #define OUTFILE_TECPLOT "outputs/Tecplot2d.dat" #endif #ifndef OUTFILE_GRID #define OUTFILE_GRID "outputs/grid2d.dat" #endif #ifndef OUTFILE_FLOW #define OUTFILE_FLOW "outputs/flow2d.dat" #endif /******************************************************************************* * * problem-specific constants: */ #ifndef NX /* Number of cells along X */ #define NX 400 #define NX1 401 #define NX4 404 #endif #ifndef NY /* Number of cells along Y */ #define NY 400 #define NY4 404 #endif #ifndef XL /* Size of domain along X */ #define XL 2d #endif #ifndef YL /* Size of domain along Y */ #define YL 2d #endif #ifndef GAM /* Ratio of specific heats */ #define GAM 1.4d #endif #ifndef NJET /* Number of points across nozzle exit */ #define NJET 200 #endif /******************************************************************************* * * algorithm configuration: */ #ifndef IADV /* time integration method */ #define IADV 3 #endif #ifndef IMUSCL /* MUSCL reconstruction method */ #define IMUSCL 1 #endif #ifndef iaxis /* Switch of plane/axisymmetric flow */ #define iaxis 0 #endif /******************************************************************************* * * derived constants: */ #define DX (XL/tod(NX)) /* Spatial increment along X */ #define DY (YL/tod(NY)) /* Spatial increment along Y */ /******************************************************************************* * * fixed constants: */ #define CFL 0.95d /* Courant-Friedrichs-Levy number */ #define MS 2.2d /* Shock wave Mach number */ /******************************************************************************* * * TVD-specific code: */ /** * * @fn double[NX], double[NY] initGrid () * * @brief generate a cell-centered grid */ double[NX], double[NY] initGrid () { x = DX * tod( iota( NX)) + 0.5; y = DY * tod( iota( NY)) + 0.5; return (x,y); } /** * * @fn fluid_pv[NX,NY] initFlow() * * @brief initialize the flow field */ fluid_pv[NX,NY] initFlow() { return( genarray( [NX,NY], genFluid_pv( [ 0d, 0d] , 1d, GAM))); } #if defined (SAVE) /** * * @fn void saveStep(double[+] x, double [+] y, fluid_pv[+] q) * * @brief saves intermediate state either to a file or to a display */ void saveStep(double[+] x, double [+] y, fluid_pv[+] q) { #if defined(DISLIN) clearDisplay( display); createAxes2D( display, tof( MINS), tof( MAXS), tof( STEPS)); printContour( display, "red", tof( rho(q)), 100); destroyAxes2D( display); #else save_grid( x,y); save_flow( q); #endif } #endif /* * Saves grid to file */ void save_grid (double[NX] x, double[NY] y) { File ff; iv,ff = fopen ( OUTFILE_GRID,"w"); for (ix=0; ix <= NX-1; ix++) fprintf(ff, "%lf \n", x[ix]); for (iy=0; iy <= NY-1; iy++) fprintf(ff, "%lf \n", y[iy]); fclose (ff); } /* * Saves flowfield to file */ void save_flow (fluid_pv[+] q) { File ff; iv,ff = fopen ( OUTFILE_FLOW,"w"); save_fluid( ff, q); fclose (ff); } /** * * @fn fluid_pv[+], double stepFlow (fluid_pv[+] qp, double t, double tk) * * @brief Advances solution in time */ fluid_pv[NX,NY], double stepFlow (fluid_pv[NX,NY] qp, double t, double tk) { while (fabs(t-tk) > 0.00000001d){ dt = getDt( qp); dt = min(dt,tk-t); #if IADV==1 qp = rktvd1 (qp, dt); #elif IADV==2 qp = rktvd2 (qp, dt); #elif IADV==3 qp = rktvd3 (qp, dt); #else printf (" Wrong value of IADV! \n"); #endif t = t + dt; printf("t = %1.16e, dt = %1.16e \n",t,dt); } return (qp,t); } /** * * @fn double getDt(fluid_pv[+] qp) * * @brief Evaluates available time step. */ double getDt(fluid_pv[NX,NY] qp) { c = sqrt(GAM * p(qp) / rho(qp)); return( CFL / maxval( ( fabs( u(qp, 0)) + c) / DX + ( fabs( u(qp, 1)) + c) / DY) ); } /** * * @fn fluid_pv[+] rktvd1 ( fluid_pv[+] qp, double dt) * * @brief Time integration with forward Euler method */ fluid_pv[NX,NY] rktvd1 ( fluid_pv[NX,NY] qp, double dt) { qc = prim2cons( qp, GAM); q1c = qc + dt * rhs( qp); return( cons2prim( q1c, GAM)); } /** * * @fn fluid_pv[+] rktvd2 ( fluid_pv[+] qp, double dt) * * @brief Time integration with 2nd order Runge-Kutta method */ fluid_pv[NX,NY] rktvd2 ( fluid_pv[NX,NY] qp, double dt) { qc = prim2cons( qp, GAM); q1c = qc + dt* rhs( qp); q2c = 0.5 * ( qc + q1c + dt* rhs( cons2prim(q1c, GAM))); return( cons2prim( q2c, GAM) ); } /** * * @fn fluid_pv[+] rktvd3 ( fluid_pv[+] qp, double dt) * * @brief Time integration with 3rd order Runge-Kutta TVD method */ fluid_pv[NX,NY] rktvd3 ( fluid_pv[NX,NY] qp, double dt) { qc = prim2cons( qp, GAM); q1c = qc + dt* rhs( qp); q2c = 0.25 * ( 3d * qc + q1c + dt* rhs( cons2prim( q1c, GAM))); q3c = ( qc + 2d*q2c + 2d*dt* rhs( cons2prim( q2c, GAM))) / 3d; return( cons2prim( q3c, GAM)); } /** * * @fn fluid_pv[.], * fluid_pv[.], * fluid_pv[.] computeBoundaries( fluid_pv[.,.] qp, double[2] norm) * * @brief Computes the boundaries l1,l2,u wrt the outermost(!) axis. * Conceptually, they are located as [ l2, l1, qp, u, u]. */ inline fluid_pv[.], fluid_pv[.], fluid_pv[.] computeBoundaries( fluid_pv[.,.] qp, double[2] norm) { shp_x = take( [1], shape(qp)); p_in = genFluid_pv( (2d*(MS-1d/MS) / (GAM+1d)) * norm, (2d*GAM*MS*MS-(GAM-1d)) / (GAM+1d), GAM*(GAM+1d)*MS*MS / ((GAM-1d)*MS*MS+2d)); l1 = genarray( [NJET], p_in) ++ reflect( drop( [NJET], qp[[2]]), norm); l2 = genarray( [NJET], p_in) ++ reflect( drop( [NJET], qp[[3]]), norm); u = qp[ shp_x - 1]; return( l1, l2, u); } /** * * @fn fluid_pv[*] reflect( fluid_pv[*] p, double[.] norm) * * @brief Computes the reflection of all fluid cells according to the norm * given. */ inline fluid_pv reflect( fluid_pv p, double[.] norm) { return( genFluid_pv( u(p) - 2d*norm*u(p), p(p), rho(p))); } inline fluid_pv[*] reflect( fluid_pv[*] p, double[.] norm) { res = with { (. <= iv <= .) : reflect( p[iv], norm); } : genarray( shape(p), zero(p)); return( res); } /** * * @fn fluid_cv[.,.] rhs( fluid_pv[.,.] qp) * * @brief Evaluates right hand side */ fluid_cv[NX,NY] rhs( fluid_pv[NX,NY] qp) { shp_x = take( [1], shape(qp)); shp_y = take( [-1], shape(qp)); l1_x, l2_x, u_x = computeBoundaries( qp, [1d, 0d]); l1_y, l2_y, u_y = computeBoundaries( transpose(qp), [0d, 1d]); dqc_x = transpose( -( dfDxNoBoundary( fluxBoundary( transpose( qp), l1_x, l2_x, u_x, u_x, [1d, 0d]), DX))); dqc_y = -( dfDxNoBoundary( fluxBoundary( qp, l1_y, l2_y, u_y, u_y, [0d, 1d]), DY)); return( dqc_x + dqc_y); } /** * * @fn fluid_cv[+] dfDxNoBoundary( fluid_cv[+] dqc, double delta) * * @brief Differentiate along the innermost axis. * As ther are no boundaries, the result becomes smaller by one element * wrt that axis! */ inline fluid_cv[.] dfDxNoBoundary( fluid_cv[.] dqc, double delta) { return( ( drop([1], dqc) - drop( [-1], dqc) ) / delta); } inline fluid_cv[*] dfDxNoBoundary( fluid_cv[*] dqc, double delta) { res = with { ( . <= iv <= .) : dfDxNoBoundary( dqc[iv], delta); } : genarray( drop( [-1], shape( dqc)), genarray( take( [-1], shape( dqc)) -1, zero( dqc))); return( res); } /** * * @fn fluid_cv[*] fluxBoundary( fluid_pv[*] qp, fluid_pv[*] qpl1, fluid_pv[*] qpl2, * fluid_pv[*] qpr1, fluid_pv[*] qpr2, double[.] norm) * * @brief computes the flux along the innermost axis * */ inline fluid_cv[*] fluxBoundary( fluid_pv[*] qp, fluid_pv[*] qpl1, fluid_pv[*] qpl2, fluid_pv[*] qpr1, fluid_pv[*] qpr2, double[.] norm) { pl, pr = musclBoundary( qp, qpl1, qpl2, qpr1, qpr2, norm); return( riemann( pl, pr, norm)); } /* * * @fn fluid_pv[.,.], * fluid_pv[.,.] musclBoundary( fluid_pv[.,.] qp, fluid_pv[.] qpl1, * fluid_pv[.] qpl2, fluid_pv[.] qpr1, * fluid_pv[.] qpr2, double[2] s) * * @brief Calls different subroutines for reconstructing cell-face values * from cell-centered ones on the innermost axis. It requires two * boundaries being supplied for each axis. * It returns arrays that contain the values left and right of the * cell-faces, respectively. Therefore, the arrays returned are bigger * by one element on the innermost axis. * * NB: this contrasts to the original implementation in several * ways! In particular, the notion of "left" and "right" has * swapped as we no longer look wrt cells but wrt cell-faces!! */ inline fluid_pv[.,.], fluid_pv[.,.] musclBoundary( fluid_pv[.,.] qp, fluid_pv[.] qpl1, fluid_pv[.] qpl2, fluid_pv[.] qpr1, fluid_pv[.] qpr2, double[2] s) { /* * Add debugging to the size of vectors going into and out of * musclboundary and muscl1 add as comments! */ nx = shape( qp)[[0]]; ny = shape( qp)[[1]]; left,right = with { (.<= iv <= .) { pl, pr = musclBoundary( qp[iv], qpl1[iv], qpl2[iv], qpr1[iv], qpr2[iv], s); } : ( pl, pr); } : ( genarray( [nx], genarray( [1 + ny], zero(qp))), genarray( [nx], genarray( [1 + ny], zero(qp))) ); return( left, right); } /* On entry to musclBoundary the length of qp is NX/NY */ inline fluid_pv[.], fluid_pv[.] musclBoundary( fluid_pv[.] qp, fluid_pv qpl1, fluid_pv qpl2, fluid_pv qpr1, fluid_pv qpr2, double[2] s) { #if IMUSCL==1 qpl, qpr= muscl1( qp, qpl1, qpr1); #elif IMUSCL==2 Fpl, Fpr= pmuscl2( Fp); #elif IMUSCL==-2 Fpl, Fpr= xmuscl2( Fp, s); #elif IMUSCL==3 qpl, qpr = weno3( qp, s[[0]], s[[1]], qpl1, qpr1); #else printf (" Wrong value of IMUSCL! \n"); qpl, qpr= muscl1( qp, qpl1, qpr1); #endif return( qpl, qpr); } /** * @fn fluid_pv[.], fluid_pv[.] muscl1 (fluid_pv[.] qp, * fluid_pv qpl1, fluid_pv qpr1) * * @brief Calculates cell-face values using 1st order piecewise constant * reconstruction */ inline fluid_pv[.], fluid_pv[.] muscl1 (fluid_pv[.] qp, fluid_pv qpl1, fluid_pv qpr1) { return( [qpl1] ++ qp , qp ++ [qpr1]); } #if 0 /** * * @fn fluid_pv[NX4], fluid_pv[NX4] pmuscl2( fluid_pv[.] Fp) * * @brief Calculates cell-face values using 2nd order MUSCL reconstruction of * primitive variables */ fluid_pv[NX4], fluid_pv[NX4] pmuscl2( fluid_pv[.] Fp) { dFp = with { ( . <= iv < .) : Fp[iv+1] - Fp[iv]; } : genarray( shape( Fp), zero( Fp)); Fpl, Fpr = with { ( . < iv < .) { gfp = minMod( dFp[iv-1], dFp[iv]); } : ( Fp[iv] - 0.5 * gfp, Fp[iv] + 0.5 * gfp); } : ( genarray( shape( Fp), zero( Fp)), genarray( shape( Fp), zero( Fp))); return( Fpl, Fpl); } #endif /** * * @fn double[+] minMod( double[+] a, double[+] b) * * @brief MIN_MOD limiter */ inline double minMod( double a, double b) { return( ( a*b < 0d ? 0d : (fabs(a) < fabs(b) ? a : b))); } inline double[+] minMod( double[+] a, double[+] b) { return( { iv -> minMod( a[iv], b[iv]) }); } inline fluid_pv minMod(fluid_pv a, fluid_pv b) { return( genFluid_pv( minMod( u(a), u(b)), minMod( p(a), p(b)), minMod( rho(a), rho(b)))); } /** * * @fn fluid_cv[*] riemann( fluid_pv[*] qpl, fluid_pv[*] qpr, double[2] s) * * @brief Evaluates numerical flux along the innermost axis using HLLE * approximate Riemann solver. */ /* NB. This is not the same as flux in tvd2d. Here qpl and qpr are * supplied already shifted so that they can be compared elemenwise. * In flux in tvd2d it is necessary to operate on qpr[i] and qpr[i+1] * for all i that are operated on. */ inline fluid_cv riemann( fluid_pv qpl, fluid_pv qpr, double[2] s) { qcl = prim2cons( qpl, GAM); fl = consPrim2flux( qcl, qpl, s); qcr = prim2cons( qpr, GAM); fr = consPrim2flux( qcr, qpr, s); unl = sum ( s * u( qpl)); cl = sqrt( GAM * p(qpl) / rho(qpl)); unr = sum ( s * u( qpr)); cr = sqrt( GAM * p(qpr) / rho(qpr)); bl = min(unl-cl,unr-cr); br = max(unl+cl,unr+cr); bm = min(bl,0d); bp = max(br,0d); f = ( bp*fl - bm*fr + bp*bm* (qcr - qcl)) / (bp-bm); return( f); } inline fluid_cv[.] riemann( fluid_pv[.] qpl, fluid_pv[.] qpr, double[2] s) { f = with { (. <= iv <= .) : riemann( qpl[iv], qpr[iv], s); } : genarray( shape( qpl), zero( [:fluid_cv])); return(f); } inline fluid_cv[.,.] riemann( fluid_pv[.,.] qpl, fluid_pv[.,.] qpr, double[2] s) { nx = shape( qpl)[[0]]; ny = shape( qpl)[[1]]; f = with { (. <= iv <= .) : riemann( qpl[iv], qpr[iv], s); } : genarray( [nx], genarray( [ny], zero( [:fluid_cv]))); return(f); } /** * * @fn fluid_cv consPrim2flux( fluid_cv q, fluid_pv qp, double[2] s) * * @brief compute flux from conservative and primitive vars. */ inline fluid_cv consPrim2flux( fluid_cv q, fluid_pv qp, double[2] s) { un = sum ( s * u( qp)); return( genFluid_cv( rho_u( q) * un + s * p( qp), un * ( e( q) + p( qp)), rho( q) * un)); } /** * * @fn int main() * * @brief Main program */ int main() { tf = 0.05d; /*tp = 0.05d;*/ tp = tf / tod( NSAVE); x,y = initGrid(); qp = initFlow(); #if defined (SAVE) saveStep( x, y, qp); #endif t = 0d; tk = t; while (t < tf){ tk = tk+tp; qp,t = stepFlow(qp,t,tk); #if defined (SAVE) printf ("\n record at t = %lf \n \n",t); saveStep( x,y,qp); #endif } return(0); } inline fluid_pv[.], fluid_pv[.] weno3 (fluid_pv[.] qp_in, double sx, double sy, fluid_pv qpl1, fluid_pv qpr1) { /* Until we can do this: qpl, qpr= muscl1( qp_in, qpl1, qpr1); we'l do this: */ /* printf( "length of qp_in is %d\n", shape( qp_in)[[0]]);-> 400 */ qp_in = [(:fluid_pv)[0d,0d,0d,0d]] ++ [qpl1] ++ qp_in ++ [qpr1]; qc = (:double[*]) prim2cons( qp_in, GAM); qp = (:double[*]) qp_in; /*qp = with { ( [0] <= [x] < [shape( qp)[[0]]] ) : [qp[[x,0]]] ++ [qp[[x,1]]] ++ [qp[[x,3]]] ++ [qp[[x,2]]]; } : modarray(qp);*/ newqp = [[:double]]; for( i=0 ; i qc[iv,.] ++ drop( [1], reverse( qp[iv,.])) }; printf("qc\n"); print( (:double[*]) qc); n1 = 1; n2 = NX+2; /* ABove is mine */ eps = [0.00000001d, 0.00000001d, 0.00000001d, 0.00000001d]; s13 = 1d/3d; s23 = 2d/3d; qpl = genarray([NX+3,4], 0d); qpr = genarray([NX+3,4], 0d); dq = with{ ([n1-1,0] <= iv <= [n2,3]) : qc[iv+[1,3]]-qc[iv+[0,3]];} : genarray([NX+3,4], 0d); for (i=n1; i <= n2; i++){ r = qc[i,3]; c2 = GAM*qc[i,4]/r; c = sqrt(c2); dunl = sx*dq[i-1,2]+sy*dq[i-1,3]; dutl = -sy*dq[i-1,2]+sx*dq[i-1,3]; dunr = sx*dq[i,2]+sy*dq[i,3]; dutr = -sy*dq[i,2]+sx*dq[i,3]; wql = [dq[i-1,1]-r*c*dunl, dq[i-1,0]-dq[i-1,1]/c2, dutl, dq[i-1,1]+r*c*dunl]; wqr = [dq[i,1]-r*c*dunr, dq[i,0]-dq[i,1]/c2, dutr, dq[i,1]+r*c*dunr]; sl = wql*wql+eps; sr = wqr*wqr+eps; sl = sl*sl; sr = sr*sr; al = s13/sl; ar = s23/sr; bl = s23/sl; br = s13/sr; gwl = (bl*wql+br*wqr)/(bl+br); gwr = (al*wql+ar*wqr)/(al+ar); gunl = 0.5d*(gwl[3]-gwl[0])/(r*c); gutl = gwl[2]; gpl = 0.5d*(gwl[0]+gwl[3]); grl = gpl/c2+gwl[1]; guxl = sx*gunl-sy*gutl; guyl = sy*gunl+sx*gutl; gql = [grl,gpl,guxl,guyl]; gunr = 0.5d*(gwr[3]-gwr[0])/(r*c); gutr = gwr[2]; gpr = 0.5d*(gwr[0]+gwr[3]); grr = gpr/c2+gwr[1]; guxr = sx*gunr-sy*gutr; guyr = sy*gunr+sx*gutr; gqr = [grr,gpr,guxr,guyr]; for (L=0; L <= 3; L++){ qpl[i,L] = qc[i,L+3]-0.5d*gql[L]; qpr[i,L] = qc[i,L+3]+0.5d*gqr[L];} } /***/ /* qpl = (:fluid_cv[.])qpl; qpr = (:fluid_cv[.])qpr; qpl = cons2prim( qpl, GAM); qpr = cons2prim( qpr, GAM); */ qpl = (:fluid_pv[.])qpl; qpr = (:fluid_pv[.])qpr; printf("qpl\n"); print( (:double[*]) qpl); printf("qpr\n"); print( (:double[*]) qpr); qpl = drop( [1], drop( [-1], qpl)); qpr = drop( [2], qpr); /***/ return (qpl,qpr); }