1 /* 2 Code for time stepping with the multi-rate Runge-Kutta method 3 4 Notes: 5 1) The general system is written as 6 Udot = F(t,U) for the nonsplit version of multi-rate RK, 7 user should give the indexes for both slow and fast components; 8 2) The general system is written as 9 Usdot = Fs(t,Us,Uf) 10 Ufdot = Ff(t,Us,Uf) for multi-rate RK with RHS splits, 11 user should partioned RHS by themselves and also provide the indexes for both slow and fast components. 12 */ 13 14 #include <petsc/private/tsimpl.h> 15 #include <petscdm.h> 16 #include <../src/ts/impls/explicit/rk/rk.h> 17 #include <../src/ts/impls/explicit/rk/mrk.h> 18 19 static PetscErrorCode TSReset_RK_MultirateNonsplit(TS ts) 20 { 21 TS_RK *rk = (TS_RK*)ts->data; 22 RKTableau tab = rk->tableau; 23 24 PetscFunctionBegin; 25 PetscCall(VecDestroy(&rk->X0)); 26 PetscCall(VecDestroyVecs(tab->s,&rk->YdotRHS_slow)); 27 PetscFunctionReturn(0); 28 } 29 30 static PetscErrorCode TSInterpolate_RK_MultirateNonsplit(TS ts,PetscReal itime,Vec X) 31 { 32 TS_RK *rk = (TS_RK*)ts->data; 33 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 34 PetscReal h = ts->time_step; 35 PetscReal tt,t; 36 PetscScalar *b; 37 const PetscReal *B = rk->tableau->binterp; 38 39 PetscFunctionBegin; 40 PetscCheck(B,PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 41 t = (itime - rk->ptime)/h; 42 PetscCall(PetscMalloc1(s,&b)); 43 for (i=0; i<s; i++) b[i] = 0; 44 for (j=0,tt=t; j<p; j++,tt*=t) { 45 for (i=0; i<s; i++) { 46 b[i] += h * B[i*p+j] * tt; 47 } 48 } 49 PetscCall(VecCopy(rk->X0,X)); 50 PetscCall(VecMAXPY(X,s,b,rk->YdotRHS_slow)); 51 PetscCall(PetscFree(b)); 52 PetscFunctionReturn(0); 53 } 54 55 static PetscErrorCode TSStepRefine_RK_MultirateNonsplit(TS ts) 56 { 57 TS previousts,subts; 58 TS_RK *rk = (TS_RK*)ts->data; 59 RKTableau tab = rk->tableau; 60 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 61 Vec vec_fast,subvec_fast; 62 const PetscInt s = tab->s; 63 const PetscReal *A = tab->A,*c = tab->c; 64 PetscScalar *w = rk->work; 65 PetscInt i,j,k; 66 PetscReal t = ts->ptime,h = ts->time_step; 67 68 PetscFunctionBegin; 69 PetscCall(VecDuplicate(ts->vec_sol,&vec_fast)); 70 previousts = rk->subts_current; 71 PetscCall(TSRHSSplitGetSubTS(rk->subts_current,"fast",&subts)); 72 PetscCall(TSRHSSplitGetSubTS(subts,"fast",&subts)); 73 for (k=0; k<rk->dtratio; k++) { 74 for (i=0; i<s; i++) { 75 PetscCall(TSInterpolate_RK_MultirateNonsplit(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i])); 76 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 77 /* update the fast components in the stage value, the slow components will be ignored, so we do not care the slow part in vec_fast */ 78 PetscCall(VecCopy(ts->vec_sol,vec_fast)); 79 PetscCall(VecMAXPY(vec_fast,i,w,YdotRHS)); 80 /* update the fast components in the stage value */ 81 PetscCall(VecGetSubVector(vec_fast,rk->is_fast,&subvec_fast)); 82 PetscCall(VecISCopy(Y[i],rk->is_fast,SCATTER_FORWARD,subvec_fast)); 83 PetscCall(VecRestoreSubVector(vec_fast,rk->is_fast,&subvec_fast)); 84 /* compute the stage RHS */ 85 PetscCall(TSComputeRHSFunction(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS[i])); 86 } 87 PetscCall(VecCopy(ts->vec_sol,vec_fast)); 88 PetscCall(TSEvaluateStep(ts,tab->order,vec_fast,NULL)); 89 /* update the fast components in the solution */ 90 PetscCall(VecGetSubVector(vec_fast,rk->is_fast,&subvec_fast)); 91 PetscCall(VecISCopy(ts->vec_sol,rk->is_fast,SCATTER_FORWARD,subvec_fast)); 92 PetscCall(VecRestoreSubVector(vec_fast,rk->is_fast,&subvec_fast)); 93 94 if (subts) { 95 Vec *YdotRHS_copy; 96 PetscCall(VecDuplicateVecs(ts->vec_sol,s,&YdotRHS_copy)); 97 rk->subts_current = rk->subts_fast; 98 ts->ptime = t+k*h/rk->dtratio; 99 ts->time_step = h/rk->dtratio; 100 PetscCall(TSRHSSplitGetIS(rk->subts_current,"fast",&rk->is_fast)); 101 for (i=0; i<s; i++) { 102 PetscCall(VecCopy(rk->YdotRHS_slow[i],YdotRHS_copy[i])); 103 PetscCall(VecCopy(YdotRHS[i],rk->YdotRHS_slow[i])); 104 } 105 106 PetscCall(TSStepRefine_RK_MultirateNonsplit(ts)); 107 108 rk->subts_current = previousts; 109 ts->ptime = t; 110 ts->time_step = h; 111 PetscCall(TSRHSSplitGetIS(previousts,"fast",&rk->is_fast)); 112 for (i=0; i<s; i++) { 113 PetscCall(VecCopy(YdotRHS_copy[i],rk->YdotRHS_slow[i])); 114 } 115 PetscCall(VecDestroyVecs(s,&YdotRHS_copy)); 116 } 117 } 118 PetscCall(VecDestroy(&vec_fast)); 119 PetscFunctionReturn(0); 120 } 121 122 static PetscErrorCode TSStep_RK_MultirateNonsplit(TS ts) 123 { 124 TS_RK *rk = (TS_RK*)ts->data; 125 RKTableau tab = rk->tableau; 126 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS,*YdotRHS_slow = rk->YdotRHS_slow; 127 Vec stage_slow,sol_slow; /* vectors store the slow components */ 128 Vec subvec_slow; /* sub vector to store the slow components */ 129 IS is_slow = rk->is_slow; 130 const PetscInt s = tab->s; 131 const PetscReal *A = tab->A,*c = tab->c; 132 PetscScalar *w = rk->work; 133 PetscInt i,j,dtratio = rk->dtratio; 134 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 135 136 PetscFunctionBegin; 137 rk->status = TS_STEP_INCOMPLETE; 138 PetscCall(VecDuplicate(ts->vec_sol,&stage_slow)); 139 PetscCall(VecDuplicate(ts->vec_sol,&sol_slow)); 140 PetscCall(VecCopy(ts->vec_sol,rk->X0)); 141 for (i=0; i<s; i++) { 142 rk->stage_time = t + h*c[i]; 143 PetscCall(TSPreStage(ts,rk->stage_time)); 144 PetscCall(VecCopy(ts->vec_sol,Y[i])); 145 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 146 PetscCall(VecMAXPY(Y[i],i,w,YdotRHS_slow)); 147 PetscCall(TSPostStage(ts,rk->stage_time,i,Y)); 148 /* compute the stage RHS */ 149 PetscCall(TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS_slow[i])); 150 } 151 /* update the slow components in the solution */ 152 rk->YdotRHS = YdotRHS_slow; 153 rk->dtratio = 1; 154 PetscCall(TSEvaluateStep(ts,tab->order,sol_slow,NULL)); 155 rk->dtratio = dtratio; 156 rk->YdotRHS = YdotRHS; 157 /* update the slow components in the solution */ 158 PetscCall(VecGetSubVector(sol_slow,is_slow,&subvec_slow)); 159 PetscCall(VecISCopy(ts->vec_sol,is_slow,SCATTER_FORWARD,subvec_slow)); 160 PetscCall(VecRestoreSubVector(sol_slow,is_slow,&subvec_slow)); 161 162 rk->subts_current = ts; 163 rk->ptime = t; 164 rk->time_step = h; 165 PetscCall(TSStepRefine_RK_MultirateNonsplit(ts)); 166 167 ts->ptime = t + ts->time_step; 168 ts->time_step = next_time_step; 169 rk->status = TS_STEP_COMPLETE; 170 171 /* free memory of work vectors */ 172 PetscCall(VecDestroy(&stage_slow)); 173 PetscCall(VecDestroy(&sol_slow)); 174 PetscFunctionReturn(0); 175 } 176 177 static PetscErrorCode TSSetUp_RK_MultirateNonsplit(TS ts) 178 { 179 TS_RK *rk = (TS_RK*)ts->data; 180 RKTableau tab = rk->tableau; 181 182 PetscFunctionBegin; 183 PetscCall(TSRHSSplitGetIS(ts,"slow",&rk->is_slow)); 184 PetscCall(TSRHSSplitGetIS(ts,"fast",&rk->is_fast)); 185 PetscCheck(rk->is_slow && rk->is_fast,PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up RHSSplits with TSRHSSplitSetIS() using split names 'slow' and 'fast' respectively in order to use multirate RK"); 186 PetscCall(TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow)); 187 PetscCall(TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast)); 188 PetscCheck(rk->subts_slow && rk->subts_fast,PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up the RHSFunctions for 'slow' and 'fast' components using TSRHSSplitSetRHSFunction() or calling TSSetRHSFunction() for each sub-TS"); 189 PetscCall(VecDuplicate(ts->vec_sol,&rk->X0)); 190 PetscCall(VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS_slow)); 191 rk->subts_current = rk->subts_fast; 192 193 ts->ops->step = TSStep_RK_MultirateNonsplit; 194 ts->ops->interpolate = TSInterpolate_RK_MultirateNonsplit; 195 PetscFunctionReturn(0); 196 } 197 198 /* 199 Copy DM from tssrc to tsdest, while keeping the original DMTS and DMSNES in tsdest. 200 */ 201 static PetscErrorCode TSCopyDM(TS tssrc,TS tsdest) 202 { 203 DM newdm,dmsrc,dmdest; 204 205 PetscFunctionBegin; 206 PetscCall(TSGetDM(tssrc,&dmsrc)); 207 PetscCall(DMClone(dmsrc,&newdm)); 208 PetscCall(TSGetDM(tsdest,&dmdest)); 209 PetscCall(DMCopyDMTS(dmdest,newdm)); 210 PetscCall(DMCopyDMSNES(dmdest,newdm)); 211 PetscCall(TSSetDM(tsdest,newdm)); 212 PetscCall(DMDestroy(&newdm)); 213 PetscFunctionReturn(0); 214 } 215 216 static PetscErrorCode TSReset_RK_MultirateSplit(TS ts) 217 { 218 TS_RK *rk = (TS_RK*)ts->data; 219 220 PetscFunctionBegin; 221 if (rk->subts_slow) { 222 PetscCall(PetscFree(rk->subts_slow->data)); 223 rk->subts_slow = NULL; 224 } 225 if (rk->subts_fast) { 226 PetscCall(PetscFree(rk->YdotRHS_fast)); 227 PetscCall(PetscFree(rk->YdotRHS_slow)); 228 PetscCall(VecDestroy(&rk->X0)); 229 PetscCall(TSReset_RK_MultirateSplit(rk->subts_fast)); 230 PetscCall(PetscFree(rk->subts_fast->data)); 231 rk->subts_fast = NULL; 232 } 233 PetscFunctionReturn(0); 234 } 235 236 static PetscErrorCode TSInterpolate_RK_MultirateSplit(TS ts,PetscReal itime,Vec X) 237 { 238 TS_RK *rk = (TS_RK*)ts->data; 239 Vec Xslow; 240 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 241 PetscReal h; 242 PetscReal tt,t; 243 PetscScalar *b; 244 const PetscReal *B = rk->tableau->binterp; 245 246 PetscFunctionBegin; 247 PetscCheck(B,PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 248 249 switch (rk->status) { 250 case TS_STEP_INCOMPLETE: 251 case TS_STEP_PENDING: 252 h = ts->time_step; 253 t = (itime - ts->ptime)/h; 254 break; 255 case TS_STEP_COMPLETE: 256 h = ts->ptime - ts->ptime_prev; 257 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 258 break; 259 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 260 } 261 PetscCall(PetscMalloc1(s,&b)); 262 for (i=0; i<s; i++) b[i] = 0; 263 for (j=0,tt=t; j<p; j++,tt*=t) { 264 for (i=0; i<s; i++) { 265 b[i] += h * B[i*p+j] * tt; 266 } 267 } 268 for (i=0; i<s; i++) { 269 PetscCall(VecGetSubVector(rk->YdotRHS[i],rk->is_slow,&rk->YdotRHS_slow[i])); 270 } 271 PetscCall(VecGetSubVector(X,rk->is_slow,&Xslow)); 272 PetscCall(VecISCopy(rk->X0,rk->is_slow,SCATTER_REVERSE,Xslow)); 273 PetscCall(VecMAXPY(Xslow,s,b,rk->YdotRHS_slow)); 274 PetscCall(VecRestoreSubVector(X,rk->is_slow,&Xslow)); 275 for (i=0; i<s; i++) { 276 PetscCall(VecRestoreSubVector(rk->YdotRHS[i],rk->is_slow,&rk->YdotRHS_slow[i])); 277 } 278 PetscCall(PetscFree(b)); 279 PetscFunctionReturn(0); 280 } 281 282 /* 283 This is for partitioned RHS multirate RK method 284 The step completion formula is 285 286 x1 = x0 + h b^T YdotRHS 287 288 */ 289 static PetscErrorCode TSEvaluateStep_RK_MultirateSplit(TS ts,PetscInt order,Vec X,PetscBool *done) 290 { 291 TS_RK *rk = (TS_RK*)ts->data; 292 RKTableau tab = rk->tableau; 293 Vec Xslow,Xfast; /* subvectors of X which store slow components and fast components respectively */ 294 PetscScalar *w = rk->work; 295 PetscReal h = ts->time_step; 296 PetscInt s = tab->s,j; 297 298 PetscFunctionBegin; 299 PetscCall(VecCopy(ts->vec_sol,X)); 300 if (rk->slow) { 301 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 302 PetscCall(VecGetSubVector(ts->vec_sol,rk->is_slow,&Xslow)); 303 PetscCall(VecMAXPY(Xslow,s,w,rk->YdotRHS_slow)); 304 PetscCall(VecRestoreSubVector(ts->vec_sol,rk->is_slow,&Xslow)); 305 } else { 306 for (j=0; j<s; j++) w[j] = h/rk->dtratio*tab->b[j]; 307 PetscCall(VecGetSubVector(X,rk->is_fast,&Xfast)); 308 PetscCall(VecMAXPY(Xfast,s,w,rk->YdotRHS_fast)); 309 PetscCall(VecRestoreSubVector(X,rk->is_fast,&Xfast)); 310 } 311 PetscFunctionReturn(0); 312 } 313 314 static PetscErrorCode TSStepRefine_RK_MultirateSplit(TS ts) 315 { 316 TS_RK *rk = (TS_RK*)ts->data; 317 TS subts_fast = rk->subts_fast,currentlevelts; 318 TS_RK *subrk_fast = (TS_RK*)subts_fast->data; 319 RKTableau tab = rk->tableau; 320 Vec *Y = rk->Y; 321 Vec *YdotRHS = rk->YdotRHS,*YdotRHS_fast = rk->YdotRHS_fast; 322 Vec Yfast,Xfast; 323 const PetscInt s = tab->s; 324 const PetscReal *A = tab->A,*c = tab->c; 325 PetscScalar *w = rk->work; 326 PetscInt i,j,k; 327 PetscReal t = ts->ptime,h = ts->time_step; 328 329 PetscFunctionBegin; 330 for (k=0; k<rk->dtratio; k++) { 331 PetscCall(VecGetSubVector(ts->vec_sol,rk->is_fast,&Xfast)); 332 for (i=0; i<s; i++) { 333 PetscCall(VecGetSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i])); 334 } 335 /* propagate fast component using small time steps */ 336 for (i=0; i<s; i++) { 337 /* stage value for slow components */ 338 PetscCall(TSInterpolate_RK_MultirateSplit(rk->ts_root,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i])); 339 currentlevelts = rk->ts_root; 340 while (currentlevelts != ts) { /* all the slow parts need to be interpolated separated */ 341 currentlevelts = ((TS_RK*)currentlevelts->data)->subts_fast; 342 PetscCall(TSInterpolate_RK_MultirateSplit(currentlevelts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i])); 343 } 344 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 345 subrk_fast->stage_time = t + h/rk->dtratio*c[i]; 346 PetscCall(TSPreStage(subts_fast,subrk_fast->stage_time)); 347 /* stage value for fast components */ 348 PetscCall(VecGetSubVector(Y[i],rk->is_fast,&Yfast)); 349 PetscCall(VecCopy(Xfast,Yfast)); 350 PetscCall(VecMAXPY(Yfast,i,w,YdotRHS_fast)); 351 PetscCall(VecRestoreSubVector(Y[i],rk->is_fast,&Yfast)); 352 PetscCall(TSPostStage(subts_fast,subrk_fast->stage_time,i,Y)); 353 /* compute the stage RHS for fast components */ 354 PetscCall(TSComputeRHSFunction(subts_fast,t+k*h*rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS_fast[i])); 355 } 356 PetscCall(VecRestoreSubVector(ts->vec_sol,rk->is_fast,&Xfast)); 357 /* update the value of fast components using fast time step */ 358 rk->slow = PETSC_FALSE; 359 PetscCall(TSEvaluateStep_RK_MultirateSplit(ts,tab->order,ts->vec_sol,NULL)); 360 for (i=0; i<s; i++) { 361 PetscCall(VecRestoreSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i])); 362 } 363 364 if (subrk_fast->subts_fast) { 365 subts_fast->ptime = t+k*h/rk->dtratio; 366 subts_fast->time_step = h/rk->dtratio; 367 PetscCall(TSStepRefine_RK_MultirateSplit(subts_fast)); 368 } 369 /* update the fast components of the solution */ 370 PetscCall(VecGetSubVector(ts->vec_sol,rk->is_fast,&Xfast)); 371 PetscCall(VecISCopy(rk->X0,rk->is_fast,SCATTER_FORWARD,Xfast)); 372 PetscCall(VecRestoreSubVector(ts->vec_sol,rk->is_fast,&Xfast)); 373 } 374 PetscFunctionReturn(0); 375 } 376 377 static PetscErrorCode TSStep_RK_MultirateSplit(TS ts) 378 { 379 TS_RK *rk = (TS_RK*)ts->data; 380 RKTableau tab = rk->tableau; 381 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 382 Vec *YdotRHS_fast = rk->YdotRHS_fast,*YdotRHS_slow = rk->YdotRHS_slow; 383 Vec Yslow,Yfast; /* subvectors store the stges of slow components and fast components respectively */ 384 const PetscInt s = tab->s; 385 const PetscReal *A = tab->A,*c = tab->c; 386 PetscScalar *w = rk->work; 387 PetscInt i,j; 388 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 389 390 PetscFunctionBegin; 391 rk->status = TS_STEP_INCOMPLETE; 392 for (i=0; i<s; i++) { 393 PetscCall(VecGetSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i])); 394 PetscCall(VecGetSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i])); 395 } 396 PetscCall(VecCopy(ts->vec_sol,rk->X0)); 397 /* propagate both slow and fast components using large time steps */ 398 for (i=0; i<s; i++) { 399 rk->stage_time = t + h*c[i]; 400 PetscCall(TSPreStage(ts,rk->stage_time)); 401 PetscCall(VecCopy(ts->vec_sol,Y[i])); 402 PetscCall(VecGetSubVector(Y[i],rk->is_fast,&Yfast)); 403 PetscCall(VecGetSubVector(Y[i],rk->is_slow,&Yslow)); 404 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 405 PetscCall(VecMAXPY(Yfast,i,w,YdotRHS_fast)); 406 PetscCall(VecMAXPY(Yslow,i,w,YdotRHS_slow)); 407 PetscCall(VecRestoreSubVector(Y[i],rk->is_fast,&Yfast)); 408 PetscCall(VecRestoreSubVector(Y[i],rk->is_slow,&Yslow)); 409 PetscCall(TSPostStage(ts,rk->stage_time,i,Y)); 410 PetscCall(TSComputeRHSFunction(rk->subts_slow,t+h*c[i],Y[i],YdotRHS_slow[i])); 411 PetscCall(TSComputeRHSFunction(rk->subts_fast,t+h*c[i],Y[i],YdotRHS_fast[i])); 412 } 413 rk->slow = PETSC_TRUE; 414 /* update the slow components of the solution using slow time step */ 415 PetscCall(TSEvaluateStep_RK_MultirateSplit(ts,tab->order,ts->vec_sol,NULL)); 416 /* YdotRHS will be used for interpolation during refinement */ 417 for (i=0; i<s; i++) { 418 PetscCall(VecRestoreSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i])); 419 PetscCall(VecRestoreSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i])); 420 } 421 422 PetscCall(TSStepRefine_RK_MultirateSplit(ts)); 423 424 ts->ptime = t + ts->time_step; 425 ts->time_step = next_time_step; 426 rk->status = TS_STEP_COMPLETE; 427 PetscFunctionReturn(0); 428 } 429 430 static PetscErrorCode TSSetUp_RK_MultirateSplit(TS ts) 431 { 432 TS_RK *rk = (TS_RK*)ts->data,*nextlevelrk,*currentlevelrk; 433 TS nextlevelts; 434 Vec X0; 435 436 PetscFunctionBegin; 437 PetscCall(TSRHSSplitGetIS(ts,"slow",&rk->is_slow)); 438 PetscCall(TSRHSSplitGetIS(ts,"fast",&rk->is_fast)); 439 PetscCheck(rk->is_slow && rk->is_fast,PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up RHSSplits with TSRHSSplitSetIS() using split names 'slow' and 'fast' respectively in order to use -ts_type bsi"); 440 PetscCall(TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow)); 441 PetscCall(TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast)); 442 PetscCheck(rk->subts_slow && rk->subts_fast,PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must set up the RHSFunctions for 'slow' and 'fast' components using TSRHSSplitSetRHSFunction() or calling TSSetRHSFunction() for each sub-TS"); 443 444 PetscCall(VecDuplicate(ts->vec_sol,&X0)); 445 /* The TS at each level share the same tableau, work array, solution vector, stage values and stage derivatives */ 446 currentlevelrk = rk; 447 while (currentlevelrk->subts_fast) { 448 PetscCall(PetscMalloc1(rk->tableau->s,¤tlevelrk->YdotRHS_fast)); 449 PetscCall(PetscMalloc1(rk->tableau->s,¤tlevelrk->YdotRHS_slow)); 450 PetscCall(PetscObjectReference((PetscObject)X0)); 451 currentlevelrk->X0 = X0; 452 currentlevelrk->ts_root = ts; 453 454 /* set up the ts for the slow part */ 455 nextlevelts = currentlevelrk->subts_slow; 456 PetscCall(PetscNewLog(nextlevelts,&nextlevelrk)); 457 nextlevelrk->tableau = rk->tableau; 458 nextlevelrk->work = rk->work; 459 nextlevelrk->Y = rk->Y; 460 nextlevelrk->YdotRHS = rk->YdotRHS; 461 nextlevelts->data = (void*)nextlevelrk; 462 PetscCall(TSCopyDM(ts,nextlevelts)); 463 PetscCall(TSSetSolution(nextlevelts,ts->vec_sol)); 464 465 /* set up the ts for the fast part */ 466 nextlevelts = currentlevelrk->subts_fast; 467 PetscCall(PetscNewLog(nextlevelts,&nextlevelrk)); 468 nextlevelrk->tableau = rk->tableau; 469 nextlevelrk->work = rk->work; 470 nextlevelrk->Y = rk->Y; 471 nextlevelrk->YdotRHS = rk->YdotRHS; 472 nextlevelrk->dtratio = rk->dtratio; 473 PetscCall(TSRHSSplitGetIS(nextlevelts,"slow",&nextlevelrk->is_slow)); 474 PetscCall(TSRHSSplitGetSubTS(nextlevelts,"slow",&nextlevelrk->subts_slow)); 475 PetscCall(TSRHSSplitGetIS(nextlevelts,"fast",&nextlevelrk->is_fast)); 476 PetscCall(TSRHSSplitGetSubTS(nextlevelts,"fast",&nextlevelrk->subts_fast)); 477 nextlevelts->data = (void*)nextlevelrk; 478 PetscCall(TSCopyDM(ts,nextlevelts)); 479 PetscCall(TSSetSolution(nextlevelts,ts->vec_sol)); 480 481 currentlevelrk = nextlevelrk; 482 } 483 PetscCall(VecDestroy(&X0)); 484 485 ts->ops->step = TSStep_RK_MultirateSplit; 486 ts->ops->evaluatestep = TSEvaluateStep_RK_MultirateSplit; 487 ts->ops->interpolate = TSInterpolate_RK_MultirateSplit; 488 PetscFunctionReturn(0); 489 } 490 491 PetscErrorCode TSRKSetMultirate_RK(TS ts,PetscBool use_multirate) 492 { 493 TS_RK *rk = (TS_RK*)ts->data; 494 495 PetscFunctionBegin; 496 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 497 rk->use_multirate = use_multirate; 498 if (use_multirate) { 499 rk->dtratio = 2; 500 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateSplit_C",TSSetUp_RK_MultirateSplit)); 501 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateSplit_C",TSReset_RK_MultirateSplit)); 502 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateNonsplit_C",TSSetUp_RK_MultirateNonsplit)); 503 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateNonsplit_C",TSReset_RK_MultirateNonsplit)); 504 } else { 505 rk->dtratio = 0; 506 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateSplit_C",NULL)); 507 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateSplit_C",NULL)); 508 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateNonsplit_C",NULL)); 509 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateNonsplit_C",NULL)); 510 } 511 PetscFunctionReturn(0); 512 } 513 514 PetscErrorCode TSRKGetMultirate_RK(TS ts,PetscBool *use_multirate) 515 { 516 TS_RK *rk = (TS_RK*)ts->data; 517 518 PetscFunctionBegin; 519 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 520 *use_multirate = rk->use_multirate; 521 PetscFunctionReturn(0); 522 } 523