Athena: integrators/integrate_1d

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file integrate_1d_ctu.c
00004  *  \brief Integrate MHD equations using 1D version of the CTU integrator.
00005  *
00006  * PURPOSE: Integrate MHD equations using 1D version of the CTU integrator.
00007  *   Updates U.[d,M1,M2,M3,E,B2c,B3c,s] in Grid structure, where U is of type
00008  *   ConsS. Adds gravitational source terms, self-gravity, and optically-thin
00009  *   cooling.
00010  *
00011  * CONTAINS PUBLIC FUNCTIONS: 
00012  * - integrate_1d_ctu()
00013  * - integrate_init_1d()
00014  * - integrate_destruct_1d() */
00015 /*============================================================================*/
00016 
00017 #include <math.h>
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include <string.h>
00021 #include "../defs.h"
00022 #include "../athena.h"
00023 #include "../globals.h"
00024 #include "prototypes.h"
00025 #include "../prototypes.h"
00026 #ifdef PARTICLES
00027 #include "../particles/particle.h"
00028 #endif
00029 
00030 #if defined(CTU_INTEGRATOR)
00031 #ifdef SPECIAL_RELATIVITY
00032 #error : The CTU integrator cannot be used for special relativity.
00033 #endif /* SPECIAL_RELATIVITY */
00034 
00035 /* The L/R states of conserved variables and fluxes at each cell face */
00036 static Cons1DS *Ul_x1Face=NULL, *Ur_x1Face=NULL, *x1Flux=NULL;
00037 
00038 /* 1D scratch vectors used by lr_states and flux functions */
00039 static Real *Bxc=NULL, *Bxi=NULL;
00040 static Prim1DS *W=NULL, *Wl=NULL, *Wr=NULL;
00041 static Cons1DS *U1d=NULL;
00042 
00043 /* Variables at t^{n+1/2} used for source terms */
00044 static Real *dhalf = NULL, *phalf = NULL;
00045 
00046 /* Variables needed for cylindrical coordinates */
00047 #ifdef CYLINDRICAL
00048 static Real *geom_src=NULL;
00049 #endif
00050 
00051 /*=========================== PUBLIC FUNCTIONS ===============================*/
00052 /*----------------------------------------------------------------------------*/
00053 /*! \fn void integrate_1d_ctu(DomainS *pD)
00054  *  \brief 1D version of CTU unsplit integrator for MHD
00055  *
00056  *   The numbering of steps follows the numbering in the 3D version.
00057  *   NOT ALL STEPS ARE NEEDED IN 1D.
00058  */
00059 
00060 void integrate_1d_ctu(DomainS *pD)
00061 {
00062   GridS *pG=(pD->Grid);
00063   Real dtodx1 = pG->dt/pG->dx1, hdtodx1 = 0.5*pG->dt/pG->dx1;
00064   int i,il,iu, is = pG->is, ie = pG->ie;
00065   int js = pG->js;
00066   int ks = pG->ks;
00067   Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic,M1h,M2h,M3h;
00068 #ifndef BAROTROPIC
00069   Real coolfl,coolfr,coolf,Eh=0.0;
00070 #endif
00071 #if defined(MHD) 
00072   Real B1ch,B2ch,B3ch;
00073 #endif
00074 #if (NSCALARS > 0)
00075   int n;
00076 #endif
00077 #ifdef SELF_GRAVITY
00078   Real gxl,gxr,flux_m1l,flux_m1r;
00079 #endif
00080 #ifdef FEEDBACK
00081   Real dt1 = 1.0/pG->dt;
00082 #endif
00083 #ifdef STATIC_MESH_REFINEMENT
00084   int ncg,npg,dim;
00085 #endif
00086 
00087 #ifdef CYLINDRICAL
00088 #ifndef ISOTHERMAL
00089   Real Pavgh;
00090 #endif
00091   Real g,gl,gr,rinv;
00092   Real hdt = 0.5*pG->dt;
00093   Real geom_src_d,geom_src_Vx,geom_src_Vy,geom_src_P,geom_src_By,geom_src_Bz;
00094   const Real *r=pG->r, *ri=pG->ri;
00095 #endif /* CYLINDRICAL */
00096   Real lsf=1.0, rsf=1.0;
00097 
00098 /* With particles, one more ghost cell must be updated in predict step */
00099 #ifdef PARTICLES
00100   Real d1;
00101   il = is - 3;
00102   iu = ie + 3;
00103 #else
00104   il = is - 1;
00105   iu = ie + 1;
00106 #endif
00107 
00108 /* Compute predictor feedback from particle drag */
00109 #ifdef FEEDBACK
00110   feedback_predictor(pG);
00111 #endif
00112 
00113 /*=== STEP 1: Compute L/R x1-interface states and 1D x1-Fluxes ===============*/
00114 
00115 /*--- Step 1a ------------------------------------------------------------------
00116  * Load 1D vector of conserved variables;  
00117  * U1d = (d, M1, M2, M3, E, B2c, B3c, s[n])
00118  */
00119 
00120   for (i=is-nghost; i<=ie+nghost; i++) {
00121     U1d[i].d  = pG->U[ks][js][i].d;
00122     U1d[i].Mx = pG->U[ks][js][i].M1;
00123     U1d[i].My = pG->U[ks][js][i].M2;
00124     U1d[i].Mz = pG->U[ks][js][i].M3;
00125 #ifndef BAROTROPIC
00126     U1d[i].E  = pG->U[ks][js][i].E;
00127 #endif /* BAROTROPIC */
00128 #ifdef MHD
00129     U1d[i].By = pG->U[ks][js][i].B2c;
00130     U1d[i].Bz = pG->U[ks][js][i].B3c;
00131     Bxc[i] = pG->U[ks][js][i].B1c;
00132     Bxi[i] = pG->B1i[ks][js][i];
00133 #endif /* MHD */
00134 #if (NSCALARS > 0)
00135     for (n=0; n<NSCALARS; n++) U1d[i].s[n] = pG->U[ks][js][i].s[n];
00136 #endif
00137   }
00138 
00139 /*--- Step 1b ------------------------------------------------------------------
00140  * Compute L and R states at X1-interfaces.
00141  */
00142 
00143   for (i=is-nghost; i<=ie+nghost; i++) {
00144     W[i] = Cons1D_to_Prim1D(&U1d[i], &Bxc[i]);
00145   }
00146 
00147   lr_states(pG,W,Bxc,pG->dt,pG->dx1,il+1,iu-1,Wl,Wr,1);
00148 
00149 /*--- Step 1c ------------------------------------------------------------------
00150  * Add source terms from static gravitational potential for 0.5*dt to L/R states
00151  */
00152 
00153   if (StaticGravPot != NULL){
00154     for (i=il+1; i<=iu; i++) {
00155       cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00156 #ifdef CYLINDRICAL
00157       gl = (*x1GravAcc)(x1vc(pG,i-1),x2,x3);
00158       gr = (*x1GravAcc)(x1vc(pG,i),x2,x3);
00159       /* APPLY GRAV. SOURCE TERMS TO V1 USING ACCELERATION FOR (dt/2) */
00160       Wl[i].Vx -= hdt*gl;
00161       Wr[i].Vx -= hdt*gr;
00162 #else
00163       phicr = (*StaticGravPot)( x1             ,x2,x3);
00164       phicl = (*StaticGravPot)((x1-    pG->dx1),x2,x3);
00165       phifc = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00166 
00167       Wl[i].Vx -= dtodx1*(phifc - phicl);
00168       Wr[i].Vx -= dtodx1*(phicr - phifc);
00169 #endif /* CYLINDRICAL */
00170     }
00171   }
00172 
00173 /*--- Step 1c (cont) -----------------------------------------------------------
00174  * Add source terms for self-gravity for 0.5*dt to L/R states
00175  */
00176 
00177 #ifdef SELF_GRAVITY
00178   for (i=il+1; i<=iu; i++) {
00179     Wl[i].Vx -= hdtodx1*(pG->Phi[ks][js][i] - pG->Phi[ks][js][i-1]);
00180     Wr[i].Vx -= hdtodx1*(pG->Phi[ks][js][i] - pG->Phi[ks][js][i-1]);
00181   }
00182 #endif
00183 
00184 /*--- Step 1c (cont) -----------------------------------------------------------
00185  * Add source terms from optically-thin cooling for 0.5*dt to L/R states
00186  */
00187 
00188 #ifndef BAROTROPIC
00189   if (CoolingFunc != NULL){
00190     for (i=il+1; i<=iu; i++) {
00191       coolfl = (*CoolingFunc)(Wl[i].d,Wl[i].P,(0.5*pG->dt));
00192       coolfr = (*CoolingFunc)(Wr[i].d,Wr[i].P,(0.5*pG->dt));
00193 
00194       Wl[i].P -= 0.5*pG->dt*Gamma_1*coolfl;
00195       Wr[i].P -= 0.5*pG->dt*Gamma_1*coolfr;
00196     }
00197   }
00198 #endif /* BAROTROPIC */
00199 
00200 /*--- Step 1c (cont) -----------------------------------------------------------
00201  * Add source terms for particle feedback for 0.5*dt to L/R states
00202  */
00203 
00204 #ifdef FEEDBACK
00205     for (i=il+1; i<=iu; i++) {
00206       d1 = 1.0/W[i-1].d;
00207       Wl[i].Vx -= pG->Coup[ks][js][i-1].fb1*d1;
00208       Wl[i].Vy -= pG->Coup[ks][js][i-1].fb2*d1;
00209       Wl[i].Vz -= pG->Coup[ks][js][i-1].fb3*d1;
00210 
00211       d1 = 1.0/W[i].d;
00212       Wr[i].Vx -= pG->Coup[ks][js][i].fb1*d1;
00213       Wr[i].Vy -= pG->Coup[ks][js][i].fb2*d1;
00214       Wr[i].Vz -= pG->Coup[ks][js][i].fb3*d1;
00215 
00216 #ifndef BAROTROPIC
00217       Wl[i].P += pG->Coup[ks][js][i-1].Eloss*Gamma_1;
00218       Wr[i].P += pG->Coup[ks][js][i].Eloss*Gamma_1;
00219 #endif
00220     }
00221 #endif /* FEEDBACK */
00222 
00223 /*--- Step 1c (cont) -----------------------------------------------------------
00224  * ADD THE GEOMETRIC SOURCE-TERMS NOW USING CELL-CENTERED PRIMITIVE
00225  * VARIABLES AT TIME t^n
00226  */
00227 
00228 #ifdef CYLINDRICAL
00229       for (i=il+1; i<=iu; i++) {
00230         // LEFT STATE GEOMETRIC SOURCE TERM (USES W[i-1])
00231         rinv = 1.0/x1vc(pG,i-1);
00232         geom_src_d  = -W[i-1].d*W[i-1].Vx*rinv;
00233         geom_src_Vx =  SQR(W[i-1].Vy);
00234         geom_src_Vy = -W[i-1].Vx*W[i-1].Vy;
00235 #ifdef MHD
00236         geom_src_Vx -= SQR(W[i-1].By)/W[i-1].d;
00237         geom_src_Vy += Bxc[i-1]*W[i-1].By/W[i-1].d;
00238         geom_src_By =  -W[i-1].Vy*Bxc[i-1]*rinv;
00239         geom_src_Bz =  -W[i-1].Vx*W[i-1].Bz*rinv;
00240 #endif /* MHD */
00241         geom_src_Vx *= rinv;
00242         geom_src_Vy *= rinv;
00243 #ifndef ISOTHERMAL
00244         geom_src_P  = -Gamma*W[i-1].P*W[i-1].Vx*rinv;
00245 #endif /* ISOTHERMAL */
00246 
00247         // ADD SOURCE TERM TO LEFT STATE
00248         Wl[i].d  += hdt*geom_src_d;
00249         Wl[i].Vx += hdt*geom_src_Vx;
00250         Wl[i].Vy += hdt*geom_src_Vy;
00251 #ifdef MHD
00252         Wl[i].By += hdt*geom_src_By;
00253         Wl[i].Bz += hdt*geom_src_Bz;
00254 #endif /* MHD */
00255 #ifndef ISOTHERMAL
00256         Wl[i].P  += hdt*geom_src_P;
00257 #endif /* ISOTHERMAL */
00258 
00259         // RIGHT STATE GEOMETRIC SOURCE TERM (USES W[i])
00260         rinv = 1.0/x1vc(pG,i);
00261         geom_src_d  = -W[i].d*W[i].Vx*rinv;
00262         geom_src_Vx =  SQR(W[i].Vy);
00263         geom_src_Vy = -W[i].Vx*W[i].Vy;
00264 #ifdef MHD
00265         geom_src_Vx -= SQR(W[i].By)/W[i].d;
00266         geom_src_Vy += Bxc[i]*W[i].By/W[i].d;
00267         geom_src_By =  -W[i].Vy*Bxc[i]*rinv;
00268         geom_src_Bz =  -W[i].Vx*W[i].Bz*rinv;
00269 #endif /* MHD */
00270         geom_src_Vx *= rinv;
00271         geom_src_Vy *= rinv;
00272 #ifndef ISOTHERMAL
00273         geom_src_P  = -Gamma*W[i].P*W[i].Vx*rinv;
00274 #endif /* ISOTHERMAL */
00275 
00276         // ADD SOURCE TERM TO RIGHT STATE
00277         Wr[i].d  += hdt*geom_src_d;
00278         Wr[i].Vx += hdt*geom_src_Vx;
00279         Wr[i].Vy += hdt*geom_src_Vy;
00280 #ifdef MHD
00281         Wr[i].By += hdt*geom_src_By;
00282         Wr[i].Bz += hdt*geom_src_Bz;
00283 #endif /* MHD */
00284 #ifndef ISOTHERMAL
00285         Wr[i].P  += hdt*geom_src_P;
00286 #endif /* ISOTHERMAL */
00287       }
00288 #endif /* CYLINDRICAL */
00289 
00290 /*--- Step 1d ------------------------------------------------------------------
00291  * Compute 1D fluxes in x1-direction
00292  */
00293 
00294   for (i=il+1; i<=iu; i++) {
00295     Ul_x1Face[i] = Prim1D_to_Cons1D(&Wl[i], &Bxi[i]);
00296     Ur_x1Face[i] = Prim1D_to_Cons1D(&Wr[i], &Bxi[i]);
00297 
00298     fluxes(Ul_x1Face[i],Ur_x1Face[i],Wl[i],Wr[i],Bxi[i],&x1Flux[i]);
00299   }
00300 
00301 /*=== STEPS 2-7: Not needed in 1D ===*/
00302 
00303 /*=== STEP 8: Compute cell-centered values at n+1/2 ==========================*/
00304 
00305 /*--- Step 8a ------------------------------------------------------------------
00306  * Calculate d^{n+1/2} (needed with static potential, or cooling)
00307  */
00308 
00309 #ifndef PARTICLES
00310   if ((StaticGravPot != NULL) || (CoolingFunc != NULL))
00311 #endif
00312   {
00313     for (i=il+1; i<=iu-1; i++) {
00314       dhalf[i] = pG->U[ks][js][i].d - hdtodx1*(x1Flux[i+1].d - x1Flux[i].d );
00315 #ifdef PARTICLES
00316       pG->Coup[ks][js][i].grid_d = dhalf[i];
00317 #endif
00318     }
00319   }
00320 
00321 /*--- Step 8b ------------------------------------------------------------------
00322  * Calculate P^{n+1/2} (needed with cooling)
00323  */
00324 
00325 #ifndef PARTICLES
00326   if (CoolingFunc != NULL)
00327 #endif /* PARTICLES */
00328   {
00329     for (i=il+1; i<=iu-1; i++) {
00330       M1h = pG->U[ks][js][i].M1 - hdtodx1*(x1Flux[i+1].Mx - x1Flux[i].Mx);
00331       M2h = pG->U[ks][js][i].M2 - hdtodx1*(x1Flux[i+1].My - x1Flux[i].My);
00332       M3h = pG->U[ks][js][i].M3 - hdtodx1*(x1Flux[i+1].Mz - x1Flux[i].Mz);
00333 #ifndef BAROTROPIC
00334       Eh  = pG->U[ks][js][i].E  - hdtodx1*(x1Flux[i+1].E  - x1Flux[i].E );
00335 #endif
00336 
00337 /* Add source terms for fixed gravitational potential */
00338       if (StaticGravPot != NULL){
00339         cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00340         phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00341         phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00342         M1h -= hdtodx1*(phir-phil)*pG->U[ks][js][i].d;
00343       }
00344 
00345 /* Add source terms due to self-gravity  */
00346 #ifdef SELF_GRAVITY
00347       phir = 0.5*(pG->Phi[ks][js][i] + pG->Phi[ks][js][i+1]);
00348       phil = 0.5*(pG->Phi[ks][js][i] + pG->Phi[ks][js][i-1]);
00349       M1h -= hdtodx1*(phir-phil)*pG->U[ks][js][i].d;
00350 #endif /* SELF_GRAVITY */
00351 
00352 /* Add the particle feedback terms */
00353 #ifdef FEEDBACK
00354       M1h -= pG->Coup[ks][js][i].fb1;
00355       M2h -= pG->Coup[ks][js][i].fb2;
00356       M3h -= pG->Coup[ks][js][i].fb3;
00357 #endif /* FEEDBACK */
00358 
00359 #ifndef BAROTROPIC
00360       phalf[i] = Eh - 0.5*(M1h*M1h + M2h*M2h + M3h*M3h)/dhalf[i];
00361 
00362 #ifdef MHD
00363       B1ch = pG->U[ks][js][i].B1c;
00364       B2ch = pG->U[ks][js][i].B2c - hdtodx1*(x1Flux[i+1].By - x1Flux[i].By);
00365       B3ch = pG->U[ks][js][i].B3c - hdtodx1*(x1Flux[i+1].Bz - x1Flux[i].Bz);
00366       phalf[i] -= 0.5*(B1ch*B1ch + B2ch*B2ch + B3ch*B3ch);
00367 #endif /* MHD */
00368 
00369       phalf[i] *= Gamma_1;
00370 #endif /* BAROTROPIC */
00371 
00372 #ifdef PARTICLES
00373       d1 = 1.0/dhalf[i];
00374       pG->Coup[ks][js][i].grid_v1 = M1h*d1;
00375       pG->Coup[ks][js][i].grid_v2 = M2h*d1;
00376       pG->Coup[ks][js][i].grid_v3 = M3h*d1;
00377 #ifndef BAROTROPIC
00378       pG->Coup[ks][js][i].grid_cs = sqrt(Gamma*phalf[i]*d1);
00379 #endif  /* BAROTROPIC */
00380 #endif /* PARTICLES */
00381 
00382     }
00383   }
00384 
00385 /*=== STEP 8.5: Integrate the particles, compute the feedback ================*/
00386 
00387 #ifdef PARTICLES
00388   Integrate_Particles(pG,pD);
00389 #ifdef FEEDBACK
00390   exchange_feedback(pG, pD);
00391 #endif
00392 #endif
00393 
00394 /*=== STEPS 9-10: Not needed in 1D ===*/
00395 
00396 /*=== STEP 11: Add source terms for a full timestep using n+1/2 states =======*/
00397 
00398 /*--- Step 11a -----------------------------------------------------------------
00399  * ADD GEOMETRIC SOURCE TERMS.
00400  */
00401 
00402 #ifdef CYLINDRICAL
00403   for (i=il+1; i<=iu-1; i++) {
00404     rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00405 
00406     /* CALCULATE DENSITY AT TIME n+1/2 */
00407     dhalf[i] = pG->U[ks][js][i].d
00408              - hdtodx1*(rsf*x1Flux[i+1].d - lsf*x1Flux[i].d);
00409 
00410     /* CALCULATE X2-MOMENTUM AT TIME n+1/2 */
00411     M2h = pG->U[ks][js][i].M2
00412         - hdtodx1*(SQR(rsf)*x1Flux[i+1].My - SQR(lsf)*x1Flux[i].My);
00413 
00414     /* COMPUTE GEOMETRIC SOURCE TERM AT TIME n+1/2 */
00415     geom_src[i] = SQR(M2h)/dhalf[i];
00416 #ifdef MHD
00417     B2ch = pG->U[ks][js][i].B2c - hdtodx1*(x1Flux[i+1].By - x1Flux[i].By);
00418     geom_src[i] -= SQR(B2ch);
00419 #endif
00420 #ifdef ISOTHERMAL
00421     geom_src[i] += Iso_csound2*dhalf[i];
00422 #ifdef MHD
00423     B1ch = pG->U[ks][js][i].B1c;
00424     B3ch = pG->U[ks][js][i].B3c - hdtodx1*(rsf*x1Flux[i+1].Bz - lsf*x1Flux[i].Bz);
00425     geom_src[i] += 0.5*(SQR(B1ch)+SQR(B2ch)+SQR(B3ch));
00426 #endif
00427 #else /* ISOTHERMAL */
00428     Pavgh = 0.5*(lsf*x1Flux[i].Pflux + rsf*x1Flux[i+1].Pflux);
00429     geom_src[i] += Pavgh;
00430 #endif
00431     geom_src[i] /= x1vc(pG,i);
00432 
00433     /* ADD TIME-CENTERED GEOMETRIC SOURCE TERM FOR FULL dt */
00434     pG->U[ks][js][i].M1 += pG->dt*geom_src[i];
00435   }
00436 #endif /* CYLINDRICAL */
00437 
00438 /*--- Step 11a -----------------------------------------------------------------
00439  * Add gravitational source terms as a Static Potential.
00440  *   The energy source terms computed at cell faces are averaged to improve
00441  * conservation of total energy.
00442  *    S_{M} = -(\rho)^{n+1/2} Grad(Phi);   S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
00443  */
00444 
00445   if (StaticGravPot != NULL){
00446     for (i=is; i<=ie; i++) {
00447       cc_pos(pG,i,js,ks,&x1,&x2,&x3);
00448       phic = (*StaticGravPot)((x1            ),x2,x3);
00449       phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
00450       phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00451 
00452 #ifdef CYLINDRICAL
00453       g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
00454       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00455       pG->U[ks][js][i].M1 -= pG->dt*dhalf[i]*g;
00456 #else
00457       pG->U[ks][js][i].M1 -= dtodx1*dhalf[i]*(phir-phil);
00458 #endif
00459 
00460 #ifndef BAROTROPIC
00461       pG->U[ks][js][i].E -= dtodx1*(lsf*x1Flux[i  ].d*(phic - phil) +
00462                                     rsf*x1Flux[i+1].d*(phir - phic));
00463 #endif
00464     }
00465   }
00466 
00467 /*--- Step 11b -----------------------------------------------------------------
00468  * Add source terms for self-gravity.
00469  * A flux correction using Phi^{n+1} in the main loop is required to make
00470  * the source terms 2nd order: see selfg_flux_correction().
00471  */
00472 
00473 #ifdef SELF_GRAVITY
00474   for (i=is; i<=ie; i++) {
00475       phic = pG->Phi[ks][js][i];
00476       phil = 0.5*(pG->Phi[ks][js][i-1] + pG->Phi[ks][js][i  ]);
00477       phir = 0.5*(pG->Phi[ks][js][i  ] + pG->Phi[ks][js][i+1]);
00478 
00479       gxl = (pG->Phi[ks][js][i-1] - pG->Phi[ks][js][i  ])/(pG->dx1);
00480       gxr = (pG->Phi[ks][js][i  ] - pG->Phi[ks][js][i+1])/(pG->dx1);
00481 
00482 /* 1-momentum flux.  2nd term is needed only if Jean's swindle used */
00483       flux_m1l = 0.5*(gxl*gxl)/four_pi_G + grav_mean_rho*phil;
00484       flux_m1r = 0.5*(gxr*gxr)/four_pi_G + grav_mean_rho*phir;
00485 
00486       pG->U[ks][js][i].M1 -= dtodx1*(flux_m1r - flux_m1l);
00487 #ifndef BAROTROPIC
00488       pG->U[ks][js][i].E -= dtodx1*(x1Flux[i  ].d*(phic - phil) +
00489                                     x1Flux[i+1].d*(phir - phic));
00490 #endif
00491   }
00492 
00493 /* Save mass fluxes in Grid structure for source term correction in main loop */
00494   for (i=is; i<=ie+1; i++) {
00495     pG->x1MassFlux[ks][js][i] = x1Flux[i].d;
00496   }
00497 #endif /* SELF_GRAVITY */
00498 
00499 /*--- Step 11c -----------------------------------------------------------------
00500  * Add source terms for optically thin cooling
00501  */
00502 
00503 #ifndef BAROTROPIC
00504   if (CoolingFunc != NULL){
00505     for (i=is; i<=ie; i++) {
00506       coolf = (*CoolingFunc)(dhalf[i],phalf[i],pG->dt);
00507       pG->U[ks][js][i].E -= pG->dt*coolf;
00508     }
00509   }
00510 #endif /* BAROTROPIC */
00511 
00512 /*--- Step 11d -----------------------------------------------------------------
00513  * Add source terms for particle feedback
00514  */
00515 
00516 #ifdef FEEDBACK
00517   for (i=is; i<=ie; i++) {
00518     pG->U[ks][js][i].M1 -= pG->Coup[ks][js][i].fb1;
00519     pG->U[ks][js][i].M2 -= pG->Coup[ks][js][i].fb2;
00520     pG->U[ks][js][i].M3 -= pG->Coup[ks][js][i].fb3;
00521 #ifndef BAROTROPIC
00522     pG->U[ks][js][i].E += pG->Coup[ks][js][i].Eloss;
00523     pG->Coup[ks][js][i].Eloss *= dt1;
00524 #endif
00525   }
00526 #endif
00527 
00528 /*=== STEP 12: Update cell-centered values for a full timestep ===============*/
00529 
00530 /*--- Step 12a -----------------------------------------------------------------
00531  * Update cell-centered variables in pG using 1D x1-fluxes
00532  */
00533 
00534   for (i=is; i<=ie; i++) {
00535 #ifdef CYLINDRICAL
00536     rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00537 #endif
00538     pG->U[ks][js][i].d  -= dtodx1*(rsf*x1Flux[i+1].d  - lsf*x1Flux[i].d );
00539     pG->U[ks][js][i].M1 -= dtodx1*(rsf*x1Flux[i+1].Mx - lsf*x1Flux[i].Mx);
00540     pG->U[ks][js][i].M2 -= dtodx1*(SQR(rsf)*x1Flux[i+1].My - SQR(lsf)*x1Flux[i].My);
00541     pG->U[ks][js][i].M3 -= dtodx1*(rsf*x1Flux[i+1].Mz - lsf*x1Flux[i].Mz);
00542 #ifndef BAROTROPIC
00543     pG->U[ks][js][i].E  -= dtodx1*(rsf*x1Flux[i+1].E  - lsf*x1Flux[i].E );
00544 #endif /* BAROTROPIC */
00545 #ifdef MHD
00546     pG->U[ks][js][i].B2c -= dtodx1*(x1Flux[i+1].By - x1Flux[i].By);
00547     pG->U[ks][js][i].B3c -= dtodx1*(rsf*x1Flux[i+1].Bz - lsf*x1Flux[i].Bz);
00548 /* For consistency, set B2i and B3i to cell-centered values.  */
00549     pG->B2i[ks][js][i] = pG->U[ks][js][i].B2c;
00550     pG->B3i[ks][js][i] = pG->U[ks][js][i].B3c;
00551 #endif /* MHD */
00552 #if (NSCALARS > 0)
00553     for (n=0; n<NSCALARS; n++)
00554       pG->U[ks][js][i].s[n] -= dtodx1*(rsf*x1Flux[i+1].s[n] - lsf*x1Flux[i].s[n]);
00555 #endif
00556   }
00557 
00558 /*--- Step 12b: Not needed in 1D ---*/
00559 /*--- Step 12c: Not needed in 1D ---*/
00560 /*--- Step 12d: Not needed in 1D ---*/
00561 
00562 #ifdef STATIC_MESH_REFINEMENT
00563 /*--- Step 12e -----------------------------------------------------------------
00564  * With SMR, store fluxes at boundaries of child and parent grids.
00565  */
00566 
00567 /* x1-boundaries of child Grids (interior to THIS Grid) */
00568   for (ncg=0; ncg<pG->NCGrid; ncg++) {
00569     for (dim=0; dim<2; dim++){
00570       if (pG->CGrid[ncg].myFlx[dim] != NULL) {
00571 
00572         if (dim==0) i = pG->CGrid[ncg].ijks[0];
00573         if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
00574 
00575         pG->CGrid[ncg].myFlx[dim][ks][js].d  = x1Flux[i].d; 
00576         pG->CGrid[ncg].myFlx[dim][ks][js].M1 = x1Flux[i].Mx; 
00577         pG->CGrid[ncg].myFlx[dim][ks][js].M2 = x1Flux[i].My;
00578         pG->CGrid[ncg].myFlx[dim][ks][js].M3 = x1Flux[i].Mz; 
00579 #ifndef BAROTROPIC
00580         pG->CGrid[ncg].myFlx[dim][ks][js].E  = x1Flux[i].E; 
00581 #endif /* BAROTROPIC */
00582 #ifdef MHD
00583         pG->CGrid[ncg].myFlx[dim][ks][js].B1c = 0.0;
00584         pG->CGrid[ncg].myFlx[dim][ks][js].B2c = x1Flux[i].By; 
00585         pG->CGrid[ncg].myFlx[dim][ks][js].B3c = x1Flux[i].Bz; 
00586 #endif /* MHD */
00587 #if (NSCALARS > 0)
00588         for (n=0; n<NSCALARS; n++)
00589           pG->CGrid[ncg].myFlx[dim][ks][js].s[n]  = x1Flux[i].s[n]; 
00590 #endif
00591       }
00592     }
00593   }
00594 
00595 /* x1-boundaries of parent Grids (at boundaries of THIS Grid)  */
00596   for (npg=0; npg<pG->NPGrid; npg++) {
00597     for (dim=0; dim<2; dim++){
00598       if (pG->PGrid[npg].myFlx[dim] != NULL) {
00599 
00600         if (dim==0) i = pG->PGrid[npg].ijks[0];
00601         if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
00602 
00603         pG->PGrid[npg].myFlx[dim][ks][js].d  = x1Flux[i].d; 
00604         pG->PGrid[npg].myFlx[dim][ks][js].M1 = x1Flux[i].Mx; 
00605         pG->PGrid[npg].myFlx[dim][ks][js].M2 = x1Flux[i].My;
00606         pG->PGrid[npg].myFlx[dim][ks][js].M3 = x1Flux[i].Mz; 
00607 #ifndef BAROTROPIC
00608         pG->PGrid[npg].myFlx[dim][ks][js].E  = x1Flux[i].E; 
00609 #endif /* BAROTROPIC */
00610 #ifdef MHD
00611         pG->PGrid[npg].myFlx[dim][ks][js].B1c = 0.0;
00612         pG->PGrid[npg].myFlx[dim][ks][js].B2c = x1Flux[i].By; 
00613         pG->PGrid[npg].myFlx[dim][ks][js].B3c = x1Flux[i].Bz; 
00614 #endif /* MHD */
00615 #if (NSCALARS > 0)
00616         for (n=0; n<NSCALARS; n++)
00617           pG->PGrid[npg].myFlx[dim][ks][js].s[n]  = x1Flux[i].s[n]; 
00618 #endif
00619       }
00620     }
00621   }
00622 
00623 #endif /* STATIC_MESH_REFINEMENT */
00624 
00625   return;
00626 }
00627 
00628 /*----------------------------------------------------------------------------*/
00629 /*! \fn void integrate_init_1d(MeshS *pM)
00630  *  \brief Allocate temporary integration arrays */
00631 void integrate_init_1d(MeshS *pM)
00632 {
00633   int size1=0,nl,nd;
00634 
00635 /* Cycle over all Grids on this processor to find maximum Nx1 */
00636   for (nl=0; nl<(pM->NLevels); nl++){
00637     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00638       if (pM->Domain[nl][nd].Grid != NULL) {
00639         if ((pM->Domain[nl][nd].Grid->Nx[0]) > size1){
00640           size1 = pM->Domain[nl][nd].Grid->Nx[0];
00641         }
00642       }
00643     }
00644   }
00645 
00646   size1 = size1 + 2*nghost;
00647 
00648   if ((Bxc = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00649   if ((Bxi = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00650 
00651   if ((U1d       =(Cons1DS*)malloc(size1*sizeof(Cons1DS)))==NULL) goto on_error;
00652   if ((Ul_x1Face =(Cons1DS*)malloc(size1*sizeof(Cons1DS)))==NULL) goto on_error;
00653   if ((Ur_x1Face =(Cons1DS*)malloc(size1*sizeof(Cons1DS)))==NULL) goto on_error;
00654   if ((x1Flux    =(Cons1DS*)malloc(size1*sizeof(Cons1DS)))==NULL) goto on_error;
00655 
00656   if ((W  = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00657   if ((Wl = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00658   if ((Wr = (Prim1DS*)malloc(size1*sizeof(Prim1DS))) == NULL) goto on_error;
00659 
00660 #ifdef CYLINDRICAL
00661   if((StaticGravPot != NULL) || (CoolingFunc != NULL))
00662 #endif
00663   {
00664     if ((dhalf  = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00665   }
00666   if(CoolingFunc != NULL){
00667     if ((phalf  = (Real*)malloc(size1*sizeof(Real))) == NULL) goto on_error;
00668   }
00669 
00670 #ifdef CYLINDRICAL
00671   if ((geom_src = (Real*)calloc_1d_array(size1, sizeof(Real))) == NULL)
00672     goto on_error;
00673 #endif
00674 
00675   return;
00676 
00677   on_error:
00678     integrate_destruct();
00679     ath_error("[integrate_init_1d]: malloc returned a NULL pointer\n");
00680 }
00681 
00682 /*----------------------------------------------------------------------------*/
00683 /*! \fn void integrate_destruct_1d(void)
00684  *  \brief Free temporary integration arrays  */
00685 void integrate_destruct_1d(void)
00686 {
00687   if (Bxc != NULL) free(Bxc);
00688   if (Bxi != NULL) free(Bxi);
00689 
00690   if (U1d != NULL) free(U1d);
00691   if (Ul_x1Face != NULL) free(Ul_x1Face);
00692   if (Ur_x1Face != NULL) free(Ur_x1Face);
00693   if (x1Flux != NULL) free(x1Flux);
00694 
00695   if (W  != NULL) free(W);
00696   if (Wl != NULL) free(Wl);
00697   if (Wr != NULL) free(Wr);
00698 
00699   if (dhalf != NULL) free(dhalf);
00700   if (phalf != NULL) free(phalf);
00701 
00702 #ifdef CYLINDRICAL
00703   if (geom_src != NULL) free_1d_array((void*)geom_src);
00704 #endif
00705 
00706   return;
00707 }
00708 #endif /* CTU_INTEGRATOR */
integrators/integrate_1d_ctu.c
