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 18 static PetscErrorCode TSReset_MRKNONSPLIT(TS ts) 19 { 20 TS_RK *rk = (TS_RK*)ts->data; 21 RKTableau tab = rk->tableau; 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 ierr = VecDestroy(&rk->X0);CHKERRQ(ierr); 26 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS_slow);CHKERRQ(ierr); 27 PetscFunctionReturn(0); 28 } 29 30 static PetscErrorCode TSInterpolate_MRKNONSPLIT(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 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 42 t = (itime - rk->ptime)/h; 43 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 44 for (i=0; i<s; i++) b[i] = 0; 45 for (j=0,tt=t; j<p; j++,tt*=t) { 46 for (i=0; i<s; i++) { 47 b[i] += h * B[i*p+j] * tt; 48 } 49 } 50 ierr = VecCopy(rk->X0,X);CHKERRQ(ierr); 51 ierr = VecMAXPY(X,s,b,rk->YdotRHS_slow);CHKERRQ(ierr); 52 ierr = PetscFree(b);CHKERRQ(ierr); 53 PetscFunctionReturn(0); 54 } 55 56 static PetscErrorCode TSStepRefine_MRKNONSPLIT(TS ts) 57 { 58 TS previousts,subts; 59 TS_RK *rk = (TS_RK*)ts->data; 60 RKTableau tab = rk->tableau; 61 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 62 Vec vec_fast,subvec_fast; 63 const PetscInt s = tab->s; 64 const PetscReal *A = tab->A,*c = tab->c; 65 PetscScalar *w = rk->work; 66 PetscInt i,j,k; 67 PetscReal t = ts->ptime,h = ts->time_step; 68 PetscErrorCode ierr; 69 70 ierr = VecDuplicate(ts->vec_sol,&vec_fast);CHKERRQ(ierr); 71 previousts = rk->subts_current; 72 ierr = TSRHSSplitGetSubTS(rk->subts_current,"fast",&subts);CHKERRQ(ierr); 73 ierr = TSRHSSplitGetSubTS(subts,"fast",&subts);CHKERRQ(ierr); 74 for (k=0; k<rk->dtratio; k++) { 75 for (i=0; i<s; i++) { 76 ierr = TSInterpolate_MRKNONSPLIT(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i]);CHKERRQ(ierr); 77 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 78 /* 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 */ 79 ierr = VecCopy(ts->vec_sol,vec_fast);CHKERRQ(ierr); 80 ierr = VecMAXPY(vec_fast,i,w,YdotRHS);CHKERRQ(ierr); 81 /* update the fast components in the stage value */ 82 ierr = VecGetSubVector(vec_fast,rk->is_fast,&subvec_fast);CHKERRQ(ierr); 83 ierr = VecISCopy(Y[i],rk->is_fast,SCATTER_FORWARD,subvec_fast);CHKERRQ(ierr); 84 ierr = VecRestoreSubVector(vec_fast,rk->is_fast,&subvec_fast);CHKERRQ(ierr); 85 /* compute the stage RHS */ 86 ierr = TSComputeRHSFunction(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 87 } 88 ierr = VecCopy(ts->vec_sol,vec_fast);CHKERRQ(ierr); 89 ierr = TSEvaluateStep(ts,tab->order,vec_fast,NULL);CHKERRQ(ierr); 90 /* update the fast components in the solution */ 91 ierr = VecGetSubVector(vec_fast,rk->is_fast,&subvec_fast);CHKERRQ(ierr); 92 ierr = VecISCopy(ts->vec_sol,rk->is_fast,SCATTER_FORWARD,subvec_fast);CHKERRQ(ierr); 93 ierr = VecRestoreSubVector(vec_fast,rk->is_fast,&subvec_fast);CHKERRQ(ierr); 94 95 if (subts) { 96 Vec *YdotRHS_copy; 97 ierr = VecDuplicateVecs(ts->vec_sol,s,&YdotRHS_copy);CHKERRQ(ierr); 98 rk->subts_current = rk->subts_fast; 99 ts->ptime = t+k*h/rk->dtratio; 100 ts->time_step = h/rk->dtratio; 101 ierr = TSRHSSplitGetIS(rk->subts_current,"fast",&rk->is_fast);CHKERRQ(ierr); 102 for (i=0; i<s; i++) { 103 ierr = VecCopy(rk->YdotRHS_slow[i],YdotRHS_copy[i]);CHKERRQ(ierr); 104 ierr = VecCopy(YdotRHS[i],rk->YdotRHS_slow[i]);CHKERRQ(ierr); 105 } 106 107 ierr = TSStepRefine_MRKNONSPLIT(ts);CHKERRQ(ierr); 108 109 rk->subts_current = previousts; 110 ts->ptime = t; 111 ts->time_step = h; 112 ierr = TSRHSSplitGetIS(previousts,"fast",&rk->is_fast);CHKERRQ(ierr); 113 for (i=0; i<s; i++) { 114 ierr = VecCopy(YdotRHS_copy[i],rk->YdotRHS_slow[i]);CHKERRQ(ierr); 115 } 116 ierr = VecDestroyVecs(s,&YdotRHS_copy);CHKERRQ(ierr); 117 } 118 } 119 ierr = VecDestroy(&vec_fast);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 static PetscErrorCode TSStep_MRKNONSPLIT(TS ts) 124 { 125 TS_RK *rk = (TS_RK*)ts->data; 126 RKTableau tab = rk->tableau; 127 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS,*YdotRHS_slow = rk->YdotRHS_slow; 128 Vec stage_slow,sol_slow; /* vectors store the slow components */ 129 Vec subvec_slow; /* sub vector to store the slow components */ 130 IS is_slow = rk->is_slow; 131 const PetscInt s = tab->s; 132 const PetscReal *A = tab->A,*c = tab->c; 133 PetscScalar *w = rk->work; 134 PetscInt i,j,dtratio = rk->dtratio; 135 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 136 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 rk->status = TS_STEP_INCOMPLETE; 140 ierr = VecDuplicate(ts->vec_sol,&stage_slow);CHKERRQ(ierr); 141 ierr = VecDuplicate(ts->vec_sol,&sol_slow);CHKERRQ(ierr); 142 ierr = VecCopy(ts->vec_sol,rk->X0);CHKERRQ(ierr); 143 for (i=0; i<s; i++) { 144 rk->stage_time = t + h*c[i]; 145 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 146 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 147 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 148 ierr = VecMAXPY(Y[i],i,w,YdotRHS_slow);CHKERRQ(ierr); 149 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 150 /* compute the stage RHS */ 151 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 152 } 153 /* update the slow components in the solution */ 154 rk->YdotRHS = YdotRHS_slow; 155 rk->dtratio = 1; 156 ierr = TSEvaluateStep(ts,tab->order,sol_slow,NULL);CHKERRQ(ierr); 157 rk->dtratio = dtratio; 158 rk->YdotRHS = YdotRHS; 159 /* update the slow components in the solution */ 160 ierr = VecGetSubVector(sol_slow,is_slow,&subvec_slow);CHKERRQ(ierr); 161 ierr = VecISCopy(ts->vec_sol,is_slow,SCATTER_FORWARD,subvec_slow);CHKERRQ(ierr); 162 ierr = VecRestoreSubVector(sol_slow,is_slow,&subvec_slow);CHKERRQ(ierr); 163 164 rk->subts_current = ts; 165 rk->ptime = t; 166 rk->time_step = h; 167 ierr = TSStepRefine_MRKNONSPLIT(ts);CHKERRQ(ierr); 168 169 ts->ptime = t + ts->time_step; 170 ts->time_step = next_time_step; 171 rk->status = TS_STEP_COMPLETE; 172 173 /* free memory of work vectors */ 174 ierr = VecDestroy(&stage_slow);CHKERRQ(ierr); 175 ierr = VecDestroy(&sol_slow);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 179 static PetscErrorCode TSSetUp_MRKNONSPLIT(TS ts) 180 { 181 TS_RK *rk = (TS_RK*)ts->data; 182 RKTableau tab = rk->tableau; 183 PetscErrorCode ierr; 184 185 PetscFunctionBegin; 186 ierr = TSRHSSplitGetIS(ts,"slow",&rk->is_slow);CHKERRQ(ierr); 187 ierr = TSRHSSplitGetIS(ts,"fast",&rk->is_fast);CHKERRQ(ierr); 188 if (!rk->is_slow || !rk->is_fast) SETERRQ(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"); 189 ierr = TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow);CHKERRQ(ierr); 190 ierr = TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast);CHKERRQ(ierr); 191 if (!rk->subts_slow || !rk->subts_fast) SETERRQ(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"); 192 ierr = VecDuplicate(ts->vec_sol,&rk->X0);CHKERRQ(ierr); 193 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS_slow);CHKERRQ(ierr); 194 rk->subts_current = rk->subts_fast; 195 196 ts->ops->step = TSStep_MRKNONSPLIT; 197 ts->ops->interpolate = TSInterpolate_MRKNONSPLIT; 198 PetscFunctionReturn(0); 199 } 200 201 /* 202 Copy DM from tssrc to tsdest, while keeping the original DMTS and DMSNES in tsdest. 203 */ 204 static PetscErrorCode TSCopyDM(TS tssrc,TS tsdest) 205 { 206 DM newdm,dmsrc,dmdest; 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 ierr = TSGetDM(tssrc,&dmsrc);CHKERRQ(ierr); 211 ierr = DMClone(dmsrc,&newdm);CHKERRQ(ierr); 212 ierr = TSGetDM(tsdest,&dmdest);CHKERRQ(ierr); 213 ierr = DMCopyDMTS(dmdest,newdm);CHKERRQ(ierr); 214 ierr = DMCopyDMSNES(dmdest,newdm);CHKERRQ(ierr); 215 ierr = TSSetDM(tsdest,newdm);CHKERRQ(ierr); 216 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 220 static PetscErrorCode TSReset_MRKSPLIT(TS ts) 221 { 222 TS_RK *rk = (TS_RK*)ts->data; 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 if (rk->subts_slow) { 227 ierr = PetscFree(rk->subts_slow->data);CHKERRQ(ierr); 228 rk->subts_slow = NULL; 229 } 230 if (rk->subts_fast) { 231 ierr = PetscFree(rk->YdotRHS_fast);CHKERRQ(ierr); 232 ierr = PetscFree(rk->YdotRHS_slow);CHKERRQ(ierr); 233 ierr = VecDestroy(&rk->X0);CHKERRQ(ierr); 234 ierr = TSReset_MRKSPLIT(rk->subts_fast); 235 ierr = PetscFree(rk->subts_fast->data);CHKERRQ(ierr); 236 rk->subts_fast = NULL; 237 } 238 PetscFunctionReturn(0); 239 } 240 241 static PetscErrorCode TSInterpolate_MRKSPLIT(TS ts,PetscReal itime,Vec X) 242 { 243 TS_RK *rk = (TS_RK*)ts->data; 244 Vec Xslow; 245 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 246 PetscReal h; 247 PetscReal tt,t; 248 PetscScalar *b; 249 const PetscReal *B = rk->tableau->binterp; 250 PetscErrorCode ierr; 251 252 PetscFunctionBegin; 253 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 254 255 switch (rk->status) { 256 case TS_STEP_INCOMPLETE: 257 case TS_STEP_PENDING: 258 h = ts->time_step; 259 t = (itime - ts->ptime)/h; 260 break; 261 case TS_STEP_COMPLETE: 262 h = ts->ptime - ts->ptime_prev; 263 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 264 break; 265 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 266 } 267 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 268 for (i=0; i<s; i++) b[i] = 0; 269 for (j=0,tt=t; j<p; j++,tt*=t) { 270 for (i=0; i<s; i++) { 271 b[i] += h * B[i*p+j] * tt; 272 } 273 } 274 for (i=0; i<s; i++) { 275 ierr = VecGetSubVector(rk->YdotRHS[i],rk->is_slow,&rk->YdotRHS_slow[i]);CHKERRQ(ierr); 276 } 277 ierr = VecGetSubVector(X,rk->is_slow,&Xslow);CHKERRQ(ierr); 278 ierr = VecISCopy(rk->X0,rk->is_slow,SCATTER_REVERSE,Xslow);CHKERRQ(ierr); 279 ierr = VecMAXPY(Xslow,s,b,rk->YdotRHS_slow);CHKERRQ(ierr); 280 ierr = VecRestoreSubVector(X,rk->is_slow,&Xslow);CHKERRQ(ierr); 281 for (i=0; i<s; i++) { 282 ierr = VecRestoreSubVector(rk->YdotRHS[i],rk->is_slow,&rk->YdotRHS_slow[i]);CHKERRQ(ierr); 283 } 284 ierr = PetscFree(b);CHKERRQ(ierr); 285 PetscFunctionReturn(0); 286 } 287 288 /* 289 This is for partitioned RHS multirate RK method 290 The step completion formula is 291 292 x1 = x0 + h b^T YdotRHS 293 294 */ 295 static PetscErrorCode TSEvaluateStep_MRKSPLIT(TS ts,PetscInt order,Vec X,PetscBool *done) 296 { 297 TS_RK *rk = (TS_RK*)ts->data; 298 RKTableau tab = rk->tableau; 299 Vec Xslow,Xfast; /* subvectors of X which store slow components and fast components respectively */ 300 PetscScalar *w = rk->work; 301 PetscReal h = ts->time_step; 302 PetscInt s = tab->s,j; 303 PetscErrorCode ierr; 304 305 PetscFunctionBegin; 306 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 307 if (rk->slow) { 308 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 309 ierr = VecGetSubVector(ts->vec_sol,rk->is_slow,&Xslow);CHKERRQ(ierr); 310 ierr = VecMAXPY(Xslow,s,w,rk->YdotRHS_slow);CHKERRQ(ierr); 311 ierr = VecRestoreSubVector(ts->vec_sol,rk->is_slow,&Xslow);CHKERRQ(ierr);; 312 } else { 313 for (j=0; j<s; j++) w[j] = h/rk->dtratio*tab->b[j]; 314 ierr = VecGetSubVector(X,rk->is_fast,&Xfast);CHKERRQ(ierr); 315 ierr = VecMAXPY(Xfast,s,w,rk->YdotRHS_fast);CHKERRQ(ierr); 316 ierr = VecRestoreSubVector(X,rk->is_fast,&Xfast);CHKERRQ(ierr); 317 } 318 PetscFunctionReturn(0); 319 } 320 321 static PetscErrorCode TSStepRefine_MRKSPLIT(TS ts) 322 { 323 TS_RK *rk = (TS_RK*)ts->data; 324 TS subts_fast = rk->subts_fast,currentlevelts; 325 TS_RK *subrk_fast = (TS_RK*)subts_fast->data; 326 RKTableau tab = rk->tableau; 327 Vec *Y = rk->Y; 328 Vec *YdotRHS = rk->YdotRHS,*YdotRHS_fast = rk->YdotRHS_fast; 329 Vec Yfast,Xfast; 330 const PetscInt s = tab->s; 331 const PetscReal *A = tab->A,*c = tab->c; 332 PetscScalar *w = rk->work; 333 PetscInt i,j,k; 334 PetscReal t = ts->ptime,h = ts->time_step; 335 PetscErrorCode ierr; 336 337 for (k=0; k<rk->dtratio; k++) { 338 ierr = VecGetSubVector(ts->vec_sol,rk->is_fast,&Xfast);CHKERRQ(ierr); 339 for (i=0; i<s; i++) { 340 ierr = VecGetSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);CHKERRQ(ierr); 341 } 342 /* propogate fast component using small time steps */ 343 for (i=0; i<s; i++) { 344 /* stage value for slow components */ 345 ierr = TSInterpolate_MRKSPLIT(rk->ts_root,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i]);CHKERRQ(ierr); 346 currentlevelts = rk->ts_root; 347 while (currentlevelts != ts) { /* all the slow parts need to be interpolated separated */ 348 currentlevelts = ((TS_RK*)currentlevelts->data)->subts_fast; 349 ierr = TSInterpolate_MRKSPLIT(currentlevelts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i]);CHKERRQ(ierr); 350 } 351 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 352 subrk_fast->stage_time = t + h/rk->dtratio*c[i]; 353 ierr = TSPreStage(subts_fast,subrk_fast->stage_time);CHKERRQ(ierr); 354 /* stage value for fast components */ 355 ierr = VecGetSubVector(Y[i],rk->is_fast,&Yfast);CHKERRQ(ierr); 356 ierr = VecCopy(Xfast,Yfast);CHKERRQ(ierr); 357 ierr = VecMAXPY(Yfast,i,w,YdotRHS_fast);CHKERRQ(ierr); 358 ierr = VecRestoreSubVector(Y[i],rk->is_fast,&Yfast);CHKERRQ(ierr); 359 ierr = TSPostStage(subts_fast,subrk_fast->stage_time,i,Y);CHKERRQ(ierr); 360 /* compute the stage RHS for fast components */ 361 ierr = TSComputeRHSFunction(subts_fast,t+k*h*rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 362 } 363 ierr = VecRestoreSubVector(ts->vec_sol,rk->is_fast,&Xfast);CHKERRQ(ierr); 364 /* update the value of fast components using fast time step */ 365 rk->slow = PETSC_FALSE; 366 ierr = TSEvaluateStep_MRKSPLIT(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 367 for (i=0; i<s; i++) { 368 ierr = VecRestoreSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);CHKERRQ(ierr); 369 } 370 371 if (subrk_fast->subts_fast) { 372 subts_fast->ptime = t+k*h/rk->dtratio; 373 subts_fast->time_step = h/rk->dtratio; 374 ierr = TSStepRefine_MRKSPLIT(subts_fast);CHKERRQ(ierr); 375 } 376 /* update the fast components of the solution */ 377 ierr = VecGetSubVector(ts->vec_sol,rk->is_fast,&Xfast);CHKERRQ(ierr); 378 ierr = VecISCopy(rk->X0,rk->is_fast,SCATTER_FORWARD,Xfast);CHKERRQ(ierr); 379 ierr = VecRestoreSubVector(ts->vec_sol,rk->is_fast,&Xfast);CHKERRQ(ierr); 380 } 381 PetscFunctionReturn(0); 382 } 383 384 static PetscErrorCode TSStep_MRKSPLIT(TS ts) 385 { 386 TS_RK *rk = (TS_RK*)ts->data; 387 RKTableau tab = rk->tableau; 388 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 389 Vec *YdotRHS_fast = rk->YdotRHS_fast,*YdotRHS_slow = rk->YdotRHS_slow; 390 Vec Yslow,Yfast; /* subvectors store the stges of slow components and fast components respectively */ 391 const PetscInt s = tab->s; 392 const PetscReal *A = tab->A,*c = tab->c; 393 PetscScalar *w = rk->work; 394 PetscInt i,j; 395 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 396 PetscErrorCode ierr; 397 398 PetscFunctionBegin; 399 rk->status = TS_STEP_INCOMPLETE; 400 for (i=0; i<s; i++) { 401 ierr = VecGetSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i]);CHKERRQ(ierr); 402 ierr = VecGetSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);CHKERRQ(ierr); 403 } 404 ierr = VecCopy(ts->vec_sol,rk->X0);CHKERRQ(ierr); 405 /* propogate both slow and fast components using large time steps */ 406 for (i=0; i<s; i++) { 407 rk->stage_time = t + h*c[i]; 408 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 409 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 410 ierr = VecGetSubVector(Y[i],rk->is_fast,&Yfast);CHKERRQ(ierr); 411 ierr = VecGetSubVector(Y[i],rk->is_slow,&Yslow);CHKERRQ(ierr); 412 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 413 ierr = VecMAXPY(Yfast,i,w,YdotRHS_fast);CHKERRQ(ierr); 414 ierr = VecMAXPY(Yslow,i,w,YdotRHS_slow);CHKERRQ(ierr); 415 ierr = VecRestoreSubVector(Y[i],rk->is_fast,&Yfast);CHKERRQ(ierr); 416 ierr = VecRestoreSubVector(Y[i],rk->is_slow,&Yslow);CHKERRQ(ierr); 417 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 418 ierr = TSComputeRHSFunction(rk->subts_slow,t+h*c[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 419 ierr = TSComputeRHSFunction(rk->subts_fast,t+h*c[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 420 } 421 rk->slow = PETSC_TRUE; 422 /* update the slow components of the solution using slow time step */ 423 ierr = TSEvaluateStep_MRKSPLIT(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 424 /* YdotRHS will be used for interpolation during refinement */ 425 for (i=0; i<s; i++) { 426 ierr = VecRestoreSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i]);CHKERRQ(ierr); 427 ierr = VecRestoreSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);CHKERRQ(ierr); 428 } 429 430 ierr = TSStepRefine_MRKSPLIT(ts);CHKERRQ(ierr); 431 432 ts->ptime = t + ts->time_step; 433 ts->time_step = next_time_step; 434 rk->status = TS_STEP_COMPLETE; 435 PetscFunctionReturn(0); 436 } 437 438 static PetscErrorCode TSSetUp_MRKSPLIT(TS ts) 439 { 440 TS_RK *rk = (TS_RK*)ts->data,*nextlevelrk,*currentlevelrk; 441 TS nextlevelts; 442 Vec X0; 443 PetscErrorCode ierr; 444 445 PetscFunctionBegin; 446 ierr = TSRHSSplitGetIS(ts,"slow",&rk->is_slow);CHKERRQ(ierr); 447 ierr = TSRHSSplitGetIS(ts,"fast",&rk->is_fast);CHKERRQ(ierr); 448 if (!rk->is_slow || !rk->is_fast) SETERRQ(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"); 449 ierr = TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow);CHKERRQ(ierr); 450 ierr = TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast);CHKERRQ(ierr); 451 if (!rk->subts_slow || !rk->subts_fast) SETERRQ(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"); 452 453 ierr = VecDuplicate(ts->vec_sol,&X0);CHKERRQ(ierr); 454 /* The TS at each level share the same tableau, work array, solution vector, stage values and stage derivatives */ 455 currentlevelrk = rk; 456 while (currentlevelrk->subts_fast) { 457 ierr = PetscMalloc1(rk->tableau->s,¤tlevelrk->YdotRHS_fast);CHKERRQ(ierr); 458 ierr = PetscMalloc1(rk->tableau->s,¤tlevelrk->YdotRHS_slow);CHKERRQ(ierr); 459 ierr = PetscObjectReference((PetscObject)X0);CHKERRQ(ierr); 460 currentlevelrk->X0 = X0; 461 currentlevelrk->ts_root = ts; 462 463 /* set up the ts for the slow part */ 464 nextlevelts = currentlevelrk->subts_slow; 465 ierr = PetscNewLog(nextlevelts,&nextlevelrk);CHKERRQ(ierr); 466 nextlevelrk->tableau = rk->tableau; 467 nextlevelrk->work = rk->work; 468 nextlevelrk->Y = rk->Y; 469 nextlevelrk->YdotRHS = rk->YdotRHS; 470 nextlevelts->data = (void*)nextlevelrk; 471 ierr = TSCopyDM(ts,nextlevelts);CHKERRQ(ierr); 472 ierr = TSSetSolution(nextlevelts,ts->vec_sol);CHKERRQ(ierr); 473 474 /* set up the ts for the fast part */ 475 nextlevelts = currentlevelrk->subts_fast; 476 ierr = PetscNewLog(nextlevelts,&nextlevelrk);CHKERRQ(ierr); 477 nextlevelrk->tableau = rk->tableau; 478 nextlevelrk->work = rk->work; 479 nextlevelrk->Y = rk->Y; 480 nextlevelrk->YdotRHS = rk->YdotRHS; 481 nextlevelrk->dtratio = rk->dtratio; 482 ierr = TSRHSSplitGetIS(nextlevelts,"slow",&nextlevelrk->is_slow);CHKERRQ(ierr); 483 ierr = TSRHSSplitGetSubTS(nextlevelts,"slow",&nextlevelrk->subts_slow);CHKERRQ(ierr); 484 ierr = TSRHSSplitGetIS(nextlevelts,"fast",&nextlevelrk->is_fast);CHKERRQ(ierr); 485 ierr = TSRHSSplitGetSubTS(nextlevelts,"fast",&nextlevelrk->subts_fast);CHKERRQ(ierr); 486 nextlevelts->data = (void*)nextlevelrk; 487 ierr = TSCopyDM(ts,nextlevelts);CHKERRQ(ierr); 488 ierr = TSSetSolution(nextlevelts,ts->vec_sol);CHKERRQ(ierr); 489 490 currentlevelrk = nextlevelrk; 491 } 492 ierr = VecDestroy(&X0);CHKERRQ(ierr); 493 494 ts->ops->step = TSStep_MRKSPLIT; 495 ts->ops->evaluatestep = TSEvaluateStep_MRKSPLIT; 496 ts->ops->interpolate = TSInterpolate_MRKSPLIT; 497 PetscFunctionReturn(0); 498 } 499 500 /*@C 501 TSRKSetMultirate - Use the interpolation-based multirate RK method 502 503 Logically collective 504 505 Input Parameter: 506 + ts - timestepping context 507 - use_multirate - PETSC_TRUE enables the multirate RK method, sets the basic method to be RK2A and sets the ratio between slow stepsize and fast stepsize to be 2 508 509 Options Database: 510 . -ts_rk_multirate - <true,false> 511 512 Notes: 513 The multirate method requires interpolation. The default interpolation works for 1st- and 2nd- order RK, but not for high-order RKs except TSRK5DP which comes with the interpolation coeffcients (binterp). 514 515 Level: intermediate 516 517 .seealso: TSRKGetMultirate() 518 @*/ 519 PetscErrorCode TSRKSetMultirate(TS ts,PetscBool use_multirate) 520 { 521 TS_RK *rk = (TS_RK*)ts->data; 522 PetscErrorCode ierr; 523 524 PetscFunctionBegin; 525 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 526 rk->use_multirate = use_multirate; 527 if (use_multirate) { 528 rk->dtratio = 2; 529 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_MRKSPLIT_C",TSSetUp_MRKSPLIT);CHKERRQ(ierr); 530 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSReset_MRKSPLIT_C",TSReset_MRKSPLIT);CHKERRQ(ierr); 531 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_MRKNONSPLIT_C",TSSetUp_MRKNONSPLIT);CHKERRQ(ierr); 532 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSReset_MRKNONSPLIT_C",TSReset_MRKNONSPLIT);CHKERRQ(ierr); 533 } else { 534 rk->dtratio = 0; 535 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_MRKSPLIT_C",NULL);CHKERRQ(ierr); 536 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSReset_MRKSPLIT_C",NULL);CHKERRQ(ierr); 537 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_MRKNONSPLIT_C",NULL);CHKERRQ(ierr); 538 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSReset_MRKNONSPLIT_C",NULL);CHKERRQ(ierr); 539 } 540 PetscFunctionReturn(0); 541 } 542 543 /*@C 544 TSRKGetMultirate - Gets whether to Use the interpolation-based multirate RK method 545 546 Not collective 547 548 Input Parameter: 549 . ts - timestepping context 550 551 Output Parameter: 552 . use_multirate - PETSC_TRUE if the multirate RK method is enabled, PETSC_FALSE otherwise 553 554 Level: intermediate 555 556 .seealso: TSRKSetMultirate() 557 @*/ 558 PetscErrorCode TSRKGetMultirate(TS ts,PetscBool *use_multirate) 559 { 560 TS_RK *rk = (TS_RK*)ts->data; 561 562 PetscFunctionBegin; 563 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 564 *use_multirate = rk->use_multirate; 565 PetscFunctionReturn(0); 566 } 567