00001 #ifdef CH_LANG_CC 00002 /* 00003 * _______ __ 00004 * / ___/ / ___ __ _ / / ___ 00005 * / /__/ _ \/ _ \/ V \/ _ \/ _ \ 00006 * \___/_//_/\___/_/_/_/_.__/\___/ 00007 * Please refer to Copyright.txt, in Chombo's root directory. 00008 */ 00009 #endif 00010 00011 #ifndef _LEVELRK4_H__ 00012 #define _LEVELRK4_H__ 00013 00014 /// Basic templated implementation of RK4 advance for a single AMRLevel 00015 /** 00016 This implements fourth-order Runge-Kutta time advance of an ODE. 00017 ODE is d(Soln)/dt = RHS 00018 00019 Template types: 00020 \li TSoln is the datatype for the solution 00021 \li TFR is a flux-register datatype 00022 \li TOp is an object encapsulating the actual ODE operations 00023 00024 TOp requirements: TOp must have the following functions: 00025 \code 00026 // evaluate d(soln)/dt at current time based on soln 00027 void evalRHS(TSoln& rhs, // d(soln)/dt based on soln 00028 TSoln& soln, // soln at current time 00029 TFR& fineFR, // flux register w/ finer level 00030 TFR& crseFR, // flux register w/ crse level 00031 const TSoln& oldCrseSoln, // old-time crse solution 00032 Real oldCrseTime, // old crse time 00033 const TSoln& newCrseSoln, // new-time crse solution 00034 Real newCrseTime, // new crse time 00035 Real time, // current time centering of soln 00036 Real fluxWeight // weight to apply to fluxRegister updates 00037 ) 00038 00039 // implements soln += dt*rhs 00040 void updateODE(TSoln& soln, 00041 const TSoln& rhs, 00042 Real dt) 00043 00044 // define data holder newSoln based on existingSoln, 00045 // including ghost cell specification 00046 void defineSolnData(TSoln& newSoln, 00047 const TSoln& existingSoln) 00048 00049 // define data holder for RHS based on existingSoln 00050 // including ghost cell specification 00051 // (which in most cases is no ghost cells) 00052 void defineRHSData(TSoln& newRHS, const TSoln& existingSoln) 00053 00054 /// copy data in src into dest 00055 void copySolnData(TSoln& dest, const TSoln& src) 00056 \endcode 00057 */ 00058 00059 template <class TSoln, class TFR, class TOp> 00060 void RK4LevelAdvance(/// the cell-centered solution at the new time (a_time + a_dt) 00061 TSoln& a_newSoln, 00062 /// the cell-centered solution at the old time (a_time) 00063 TSoln& a_oldSoln, 00064 /// old-time coarse-level solution (if it exists) 00065 const TSoln& a_oldCrseSoln, 00066 /// time-centering of a_oldCrseSoln 00067 Real a_oldCrseTime, 00068 /// new-time coarse-level solution (if it exists) 00069 const TSoln& a_newCrseSoln, 00070 /// time-centering of a_newCrseSoln 00071 Real a_newCrseTime, 00072 /// the pointer to the flux register between this level and the next coarser level 00073 TFR& a_crseFRPtr, 00074 /// the pointer to the flux register between this level and the next finer level 00075 TFR& a_fineFRPtr, 00076 /// time centering of a_oldSoln 00077 Real a_time, 00078 /// time step 00079 Real a_dt, 00080 /// object which encapsulates the ODE functionality (must meet the requirements detailed above) 00081 TOp& a_op, 00082 /// whether a_newSoln should be initialized with a_oldSoln (default true) 00083 const bool a_initializeNewSoln = true) 00084 { 00085 // Advance by one time step using 4th-order Runge-Kutta. 00086 00087 // allocate and define temp storage 00088 TSoln RHSTmp; 00089 TSoln solnTmp; 00090 a_op.defineRHSData(RHSTmp,a_oldSoln); 00091 a_op.defineSolnData(solnTmp,a_oldSoln); 00092 a_op.copySolnData(solnTmp,a_oldSoln); 00093 00094 // start by copying old solution into new solution 00095 if (a_initializeNewSoln) 00096 { 00097 a_op.copySolnData(a_newSoln, a_oldSoln); 00098 } 00099 00100 // RK Stage 1: compute RHS.at old time 00101 // RHSTmp gets RHS(tOld) 00102 // flux registers are also incremented if appropriate 00103 a_op.evalRHS(RHSTmp,a_oldSoln, 00104 a_fineFRPtr,a_crseFRPtr, 00105 a_oldCrseSoln, a_oldCrseTime, 00106 a_newCrseSoln, a_newCrseTime, 00107 a_time,a_dt/6); 00108 00109 // RK Stage 1: compute update {phi}_1, partial update to new vals. 00110 00111 // first part of update... 00112 a_op.updateODE(a_newSoln,RHSTmp,a_dt/6); 00113 00114 00115 // predict half-time value using LofPhiTmp 00116 a_op.updateODE(solnTmp, RHSTmp,a_dt/2); 00117 00118 // RK Stage 2: compute RHS based on first predicted value 00119 // also increment flux registers 00120 a_op.evalRHS(RHSTmp,solnTmp, 00121 a_fineFRPtr,a_crseFRPtr, 00122 a_oldCrseSoln, a_oldCrseTime, 00123 a_newCrseSoln, a_newCrseTime, 00124 a_time+a_dt/2,a_dt/3); 00125 00126 00127 // RK Stage 2: compute update {phi}_2, partial update to new vals. 00128 00129 a_op.updateODE(a_newSoln,RHSTmp,a_dt/3); 00130 00131 // now predict half-time value using latest estimate of RHS 00132 a_op.copySolnData(solnTmp,a_oldSoln); 00133 a_op.updateODE(solnTmp,RHSTmp,a_dt/2); 00134 00135 // RK Stage 3: compute new RHS. 00136 00137 a_op.evalRHS(RHSTmp,solnTmp, 00138 a_fineFRPtr,a_crseFRPtr, 00139 a_oldCrseSoln, a_oldCrseTime, 00140 a_newCrseSoln, a_newCrseTime, 00141 a_time+a_dt/2,a_dt/3); 00142 00143 00144 // RK Stage 3: compute update {phi}_3, partial update to new vals. 00145 00146 a_op.updateODE(a_newSoln,RHSTmp,a_dt/3); 00147 00148 // predict value at new time 00149 a_op.copySolnData(solnTmp,a_oldSoln); 00150 a_op.updateODE(solnTmp,RHSTmp,a_dt); 00151 00152 // RK Stage 4: compute RHS. 00153 00154 a_op.evalRHS(RHSTmp,solnTmp, 00155 a_fineFRPtr,a_crseFRPtr, 00156 a_oldCrseSoln, a_oldCrseTime, 00157 a_newCrseSoln, a_newCrseTime, 00158 a_time+a_dt,a_dt/6); 00159 00160 // RK Stage 4: compute final update of solution. 00161 00162 a_op.updateODE(a_newSoln,RHSTmp,a_dt/6); 00163 } 00164 00165 00166 /// Templated implementation of RK4 advance for a single AMRLevel, allowing interpolation in time 00167 /** 00168 This implements fourth-order Runge-Kutta time advance of an ODE, 00169 and saves intermediate results for use in interpolation in time. 00170 ODE is d(Soln)/dt = RHS 00171 00172 Template types: 00173 \li TSoln is the datatype for the solution 00174 \li TInterp is a storage class for intermediate values 00175 \li TFR is a flux-register datatype 00176 \li TOp is an object encapsulating the actual ODE operations 00177 00178 TInterp requirements: TInterp must have the following functions: 00179 \code 00180 // set time step 00181 void setDt(Real dt) 00182 00183 // save initial solution 00184 void saveInitialSoln(TSoln& soln) 00185 00186 // save this value of d(soln)/dt 00187 void saveRHS(TSoln& rhs) 00188 \endcode 00189 TOp requirements: TOp must have the following functions: 00190 \code 00191 // evaluate d(soln)/dt at current time based on soln 00192 void evalRHS(TSoln& rhs, // d(soln)/dt based on soln 00193 TSoln& soln, // soln at current time 00194 TFR& fineFR, // flux register w/ finer level 00195 TFR& crseFR, // flux register w/ crse level 00196 const TSoln& oldCrseSoln, // old-time crse solution 00197 Real oldCrseTime, // old crse time 00198 const TSoln& newCrseSoln, // new-time crse solution 00199 Real newCrseTime, // new crse time 00200 Real time, // current time centering of soln 00201 Real fluxWeight // weight to apply to fluxRegister updates 00202 ) 00203 00204 // implements soln += dt*rhs 00205 void updateODE(TSoln& soln, 00206 const TSoln& rhs, 00207 Real dt) 00208 00209 // define data holder newSoln based on existingSoln, 00210 // including ghost cell specification 00211 void defineSolnData(TSoln& newSoln, 00212 const TSoln& existingSoln) 00213 00214 // define data holder for RHS based on existingSoln 00215 // including ghost cell specification 00216 // (which in most cases is no ghost cells) 00217 void defineRHSData(TSoln& newRHS, const TSoln& existingSoln) 00218 00219 /// copy data in src into dest 00220 void copySolnData(TSoln& dest, const TSoln& src) 00221 \endcode 00222 */ 00223 00224 template <class TSoln, class TInterp, class TFR, class TOp> 00225 void RK4LevelAdvance(/// the cell-centered solution at the new time (a_time + a_dt) 00226 TSoln& a_newSoln, 00227 /// the cell-centered solution at the old time (a_time) 00228 TSoln& a_oldSoln, 00229 /// object that encapsulates time interpolation 00230 TInterp& a_interpPtr, 00231 /// old-time coarse-level solution (if it exists) 00232 const TSoln& a_oldCrseSoln, 00233 /// time-centering of a_oldCrseSoln 00234 Real a_oldCrseTime, 00235 /// new-time coarse-level solution (if it exists) 00236 const TSoln& a_newCrseSoln, 00237 /// time-centering of a_newCrseSoln 00238 Real a_newCrseTime, 00239 /// the pointer to the flux register between this level and the next coarser level 00240 TFR& a_crseFRPtr, 00241 /// the pointer to the flux register between this level and the next finer level 00242 TFR& a_fineFRPtr, 00243 /// time centering of a_oldSoln 00244 Real a_time, 00245 /// time step 00246 Real a_dt, 00247 /// object which encapsulates the ODE functionality (must meet the requirements detailed above) 00248 TOp& a_op, 00249 /// whether a_newSoln is initialized with a_oldSoln (default true) 00250 const bool a_initializeNewSoln = true) 00251 00252 00253 { 00254 00255 // Advance by one time step using 4th-order Runge-Kutta. 00256 00257 // allocate and define temp storage 00258 TSoln RHSTmp; 00259 TSoln solnTmp; 00260 a_op.defineRHSData(RHSTmp,a_oldSoln); 00261 a_op.defineSolnData(solnTmp,a_oldSoln); 00262 a_op.copySolnData(solnTmp,a_oldSoln); 00263 00264 // start by copying old solution into new solution 00265 if (a_initializeNewSoln) 00266 { 00267 a_op.copySolnData(a_newSoln, a_oldSoln); 00268 } 00269 a_interpPtr.setDt(a_dt); 00270 a_interpPtr.saveInitialSoln(a_oldSoln); 00271 00272 // RK Stage 1: compute RHS.at old time 00273 // RHSTmp gets RHS(tOld) 00274 // flux registers are also incremented if appropriate 00275 a_op.evalRHS(RHSTmp,a_oldSoln, 00276 a_fineFRPtr,a_crseFRPtr, 00277 a_oldCrseSoln, a_oldCrseTime, 00278 a_newCrseSoln, a_newCrseTime, 00279 a_time,a_dt/6); 00280 a_interpPtr.saveRHS(RHSTmp); 00281 00282 // RK Stage 1: compute update {phi}_1, partial update to new vals. 00283 00284 // first part of update... 00285 a_op.updateODE(a_newSoln,RHSTmp,a_dt/6); 00286 00287 00288 // predict half-time value using LofPhiTmp 00289 a_op.updateODE(solnTmp, RHSTmp,a_dt/2); 00290 00291 // RK Stage 2: compute RHS based on first predicted value 00292 // also increment flux registers 00293 a_op.evalRHS(RHSTmp,solnTmp, 00294 a_fineFRPtr,a_crseFRPtr, 00295 a_oldCrseSoln, a_oldCrseTime, 00296 a_newCrseSoln, a_newCrseTime, 00297 a_time+a_dt/2,a_dt/3); 00298 a_interpPtr.saveRHS(RHSTmp); 00299 00300 00301 // RK Stage 2: compute update {phi}_2, partial update to new vals. 00302 a_op.updateODE(a_newSoln,RHSTmp,a_dt/3); 00303 00304 // now predict half-time value using latest estimate of RHS 00305 a_op.copySolnData(solnTmp,a_oldSoln); 00306 a_op.updateODE(solnTmp,RHSTmp,a_dt/2); 00307 00308 // RK Stage 3: compute new RHS. 00309 00310 a_op.evalRHS(RHSTmp,solnTmp, 00311 a_fineFRPtr,a_crseFRPtr, 00312 a_oldCrseSoln, a_oldCrseTime, 00313 a_newCrseSoln, a_newCrseTime, 00314 a_time+a_dt/2,a_dt/3); 00315 a_interpPtr.saveRHS(RHSTmp); 00316 00317 00318 // RK Stage 3: compute update {phi}_3, partial update to new vals. 00319 00320 a_op.updateODE(a_newSoln,RHSTmp,a_dt/3); 00321 00322 // predict value at new time 00323 a_op.copySolnData(solnTmp,a_oldSoln); 00324 a_op.updateODE(solnTmp,RHSTmp,a_dt); 00325 00326 // RK Stage 4: compute RHS. 00327 00328 a_op.evalRHS(RHSTmp,solnTmp, 00329 a_fineFRPtr,a_crseFRPtr, 00330 a_oldCrseSoln, a_oldCrseTime, 00331 a_newCrseSoln, a_newCrseTime, 00332 a_time+a_dt,a_dt/6); 00333 a_interpPtr.saveRHS(RHSTmp); 00334 00335 // RK Stage 4: compute final update of solution. 00336 00337 a_op.updateODE(a_newSoln,RHSTmp,a_dt/6); 00338 } 00339 00340 #endif // end multiple-include preventer
1.5.5