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 TSRKType TSMRKDefault = TSRK2A; 19 20 static PetscErrorCode TSSetUp_MRKSPLIT(TS ts) 21 { 22 TS_RK *rk = (TS_RK*)ts->data,*nextlevelrk,*currentlevelrk; 23 TS nextlevelts; 24 DM dm,subdm,newdm; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 ierr = TSRHSSplitGetIS(ts,"slow",&rk->is_slow);CHKERRQ(ierr); 29 ierr = TSRHSSplitGetIS(ts,"fast",&rk->is_fast);CHKERRQ(ierr); 30 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"); 31 ierr = TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow);CHKERRQ(ierr); 32 ierr = TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast);CHKERRQ(ierr); 33 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"); 34 35 /* Only copy */ 36 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 37 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 38 ierr = TSGetDM(rk->subts_fast,&subdm);CHKERRQ(ierr); 39 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 40 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 41 ierr = TSSetDM(rk->subts_fast,newdm);CHKERRQ(ierr); 42 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 43 ierr = DMClone(dm,&newdm);CHKERRQ(ierr); 44 ierr = TSGetDM(rk->subts_slow,&subdm);CHKERRQ(ierr); 45 ierr = DMCopyDMTS(subdm,newdm);CHKERRQ(ierr); 46 ierr = DMCopyDMSNES(subdm,newdm);CHKERRQ(ierr); 47 ierr = TSSetDM(rk->subts_slow,newdm);CHKERRQ(ierr); 48 ierr = DMDestroy(&newdm);CHKERRQ(ierr); 49 50 ierr = VecDuplicate(ts->vec_sol,&rk->X0);CHKERRQ(ierr); 51 /* The TS at each level share the same tableau, work array, solution vector, stage values and stage derivatives */ 52 currentlevelrk = rk; 53 while (currentlevelrk->subts_fast) { 54 /* set up the ts for the slow part */ 55 nextlevelts = currentlevelrk->subts_slow; 56 ierr = PetscNewLog(nextlevelts,&nextlevelrk);CHKERRQ(ierr); 57 nextlevelrk->tableau = rk->tableau; 58 nextlevelrk->work = rk->work; 59 nextlevelrk->Y = rk->Y; 60 nextlevelrk->X0 = rk->X0; 61 ierr = PetscMalloc1(rk->tableau->s,&nextlevelrk->YdotRHS);CHKERRQ(ierr); 62 nextlevelts->data = (void*)nextlevelrk; 63 ierr = TSSetSolution(nextlevelts,ts->vec_sol);CHKERRQ(ierr); 64 65 /* set up the ts for the fast part */ 66 nextlevelts = currentlevelrk->subts_fast; 67 ierr = PetscNewLog(nextlevelts,&nextlevelrk);CHKERRQ(ierr); 68 nextlevelrk->tableau = rk->tableau; 69 nextlevelrk->work = rk->work; 70 nextlevelrk->Y = rk->Y; 71 nextlevelrk->X0 = rk->X0; 72 nextlevelrk->dtratio = rk->dtratio; 73 ierr = PetscMalloc1(rk->tableau->s,&nextlevelrk->YdotRHS);CHKERRQ(ierr); 74 ierr = TSRHSSplitGetIS(nextlevelts,"slow",&nextlevelrk->is_slow);CHKERRQ(ierr); 75 ierr = TSRHSSplitGetSubTS(nextlevelts,"slow",&nextlevelrk->subts_slow);CHKERRQ(ierr); 76 ierr = TSRHSSplitGetIS(nextlevelts,"fast",&nextlevelrk->is_fast);CHKERRQ(ierr); 77 ierr = TSRHSSplitGetSubTS(nextlevelts,"fast",&nextlevelrk->subts_fast);CHKERRQ(ierr); 78 nextlevelts->data = (void*)nextlevelrk; 79 ierr = TSSetSolution(nextlevelts,ts->vec_sol);CHKERRQ(ierr); 80 81 currentlevelrk = nextlevelrk; 82 } 83 PetscFunctionReturn(0); 84 } 85 86 static PetscErrorCode TSReset_MRKSPLIT(TS ts) 87 { 88 TS_RK *rk = (TS_RK*)ts->data; 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 if (rk->subts_slow) { 93 ierr = TSReset_MRKSPLIT(rk->subts_slow);CHKERRQ(ierr); 94 ierr = PetscFree(rk->subts_slow->data);CHKERRQ(ierr); 95 rk->subts_slow = NULL; 96 } 97 if (rk->subts_fast) { 98 ierr = VecDestroy(&rk->X0);CHKERRQ(ierr); 99 ierr = TSReset_MRKSPLIT(rk->subts_fast); 100 ierr = PetscFree(rk->subts_fast->data);CHKERRQ(ierr); 101 rk->subts_fast = NULL; 102 } 103 ierr = PetscFree(rk->YdotRHS_fast);CHKERRQ(ierr); 104 ierr = PetscFree(rk->YdotRHS_slow);CHKERRQ(ierr); 105 ierr = PetscFree(rk->YdotRHS);CHKERRQ(ierr); 106 PetscFunctionReturn(0); 107 } 108 109 static PetscErrorCode TSSetUp_MRKNONSPLIT(TS ts) 110 { 111 TS_RK *rk = (TS_RK*)ts->data; 112 RKTableau tab = rk->tableau; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 ierr = TSRHSSplitGetIS(ts,"slow",&rk->is_slow);CHKERRQ(ierr); 117 ierr = TSRHSSplitGetIS(ts,"fast",&rk->is_fast);CHKERRQ(ierr); 118 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"); 119 ierr = VecDuplicate(ts->vec_sol,&rk->X0);CHKERRQ(ierr); 120 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS_slow);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 124 static PetscErrorCode TSReset_MRKNONSPLIT(TS ts) 125 { 126 TS_RK *rk = (TS_RK*)ts->data; 127 RKTableau tab = rk->tableau; 128 PetscErrorCode ierr; 129 130 PetscFunctionBegin; 131 ierr = VecDestroy(&rk->X0);CHKERRQ(ierr); 132 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS_slow);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 static PetscErrorCode TSInterpolate_MRKNONSPLIT(TS ts,PetscReal itime,Vec X) 137 { 138 TS_RK *rk = (TS_RK*)ts->data; 139 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 140 PetscReal h; 141 PetscReal tt,t; 142 PetscScalar *b; 143 const PetscReal *B = rk->tableau->binterp; 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 148 149 switch (rk->status) { 150 case TS_STEP_INCOMPLETE: 151 case TS_STEP_PENDING: 152 h = ts->time_step; 153 t = (itime - ts->ptime)/h; 154 break; 155 case TS_STEP_COMPLETE: 156 h = ts->ptime - ts->ptime_prev; 157 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 158 break; 159 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 160 } 161 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 162 for (i=0; i<s; i++) b[i] = 0; 163 for (j=0,tt=t; j<p; j++,tt*=t) { 164 for (i=0; i<s; i++) { 165 b[i] += h * B[i*p+j] * tt; 166 } 167 } 168 ierr = VecCopy(rk->X0,X);CHKERRQ(ierr); 169 ierr = VecMAXPY(X,s,b,rk->YdotRHS_slow);CHKERRQ(ierr); 170 ierr = PetscFree(b);CHKERRQ(ierr); 171 PetscFunctionReturn(0); 172 } 173 174 static PetscErrorCode TSInterpolate_MRKSPLIT(TS ts,PetscReal itime,Vec X) 175 { 176 TS_RK *rk = (TS_RK*)ts->data; 177 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 178 Vec Yslow; /* vector holds the slow components in Y[0] */ 179 PetscReal h; 180 PetscReal tt,t; 181 PetscScalar *b; 182 const PetscReal *B = rk->tableau->binterp; 183 PetscErrorCode ierr; 184 185 PetscFunctionBegin; 186 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 187 188 switch (rk->status) { 189 case TS_STEP_INCOMPLETE: 190 case TS_STEP_PENDING: 191 h = ts->time_step; 192 t = (itime - ts->ptime)/h; 193 break; 194 case TS_STEP_COMPLETE: 195 h = ts->ptime - ts->ptime_prev; 196 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 197 break; 198 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 199 } 200 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 201 for (i=0; i<s; i++) b[i] = 0; 202 for (j=0,tt=t; j<p; j++,tt*=t) { 203 for (i=0; i<s; i++) { 204 b[i] += h * B[i*p+j] * tt; 205 } 206 } 207 ierr = VecGetSubVector(rk->X0,rk->is_slow,&Yslow);CHKERRQ(ierr); 208 ierr = VecCopy(Yslow,X);CHKERRQ(ierr); 209 ierr = VecMAXPY(X,s,b,((TS_RK*)rk->subts_slow->data)->YdotRHS);CHKERRQ(ierr); 210 ierr = VecRestoreSubVector(rk->X0,rk->is_slow,&Yslow);CHKERRQ(ierr); 211 ierr = PetscFree(b);CHKERRQ(ierr); 212 PetscFunctionReturn(0); 213 } 214 215 static PetscErrorCode TSStep_MRKNONSPLIT(TS ts) 216 { 217 TS_RK *rk = (TS_RK*)ts->data; 218 RKTableau tab = rk->tableau; 219 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS,*YdotRHS_slow = rk->YdotRHS_slow; 220 Vec stage_slow,sol_slow; /* vectors store the slow components */ 221 Vec subvec_slow; /* sub vector to store the slow components */ 222 const PetscInt s = tab->s; 223 const PetscReal *A = tab->A,*c = tab->c; 224 PetscScalar *w = rk->work; 225 PetscInt i,j,k,dtratio = rk->dtratio; 226 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 227 PetscErrorCode ierr; 228 229 PetscFunctionBegin; 230 rk->status = TS_STEP_INCOMPLETE; 231 ierr = VecDuplicate(ts->vec_sol,&stage_slow);CHKERRQ(ierr); 232 ierr = VecDuplicate(ts->vec_sol,&sol_slow);CHKERRQ(ierr); 233 ierr = VecCopy(ts->vec_sol,rk->X0);CHKERRQ(ierr); 234 for (i=0; i<s; i++) { 235 rk->stage_time = t + h*c[i]; 236 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 237 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 238 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 239 ierr = VecMAXPY(Y[i],i,w,YdotRHS_slow);CHKERRQ(ierr); 240 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 241 /* compute the stage RHS */ 242 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 243 } 244 /* update the slow components in the solution */ 245 rk->YdotRHS = YdotRHS_slow; 246 rk->dtratio = 1; 247 ierr = TSEvaluateStep(ts,tab->order,sol_slow,NULL);CHKERRQ(ierr); 248 rk->dtratio = dtratio; 249 rk->YdotRHS = YdotRHS; 250 for (k=0; k<rk->dtratio; k++) { 251 for (i=0; i<s; i++) { 252 rk->stage_time = t + h/rk->dtratio*c[i]; 253 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 254 /* update the fast components in the stage value, the slow components will be overwritten, so it is ok to have garbage in the slow components */ 255 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 256 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 257 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 258 ierr = TSInterpolate_MRKNONSPLIT(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],stage_slow);CHKERRQ(ierr); 259 /* update the slow components in the stage value */ 260 ierr = VecGetSubVector(stage_slow,rk->is_slow,&subvec_slow);CHKERRQ(ierr); 261 ierr = VecISCopy(Y[i],rk->is_slow,SCATTER_FORWARD,subvec_slow);CHKERRQ(ierr); 262 ierr = VecRestoreSubVector(stage_slow,rk->is_slow,&subvec_slow);CHKERRQ(ierr); 263 ierr = TSPostStage(ts,rk->stage_time,i,Y);CHKERRQ(ierr); 264 /* compute the stage RHS */ 265 ierr = TSComputeRHSFunction(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 266 } 267 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 268 } 269 /* update the slow components in the solution */ 270 ierr = VecGetSubVector(sol_slow,rk->is_slow,&subvec_slow);CHKERRQ(ierr); 271 ierr = VecISCopy(ts->vec_sol,rk->is_slow,SCATTER_FORWARD,subvec_slow);CHKERRQ(ierr); 272 ierr = VecRestoreSubVector(sol_slow,rk->is_slow,&subvec_slow);CHKERRQ(ierr); 273 274 ts->ptime += ts->time_step; 275 ts->time_step = next_time_step; 276 rk->status = TS_STEP_COMPLETE; 277 /* free memory of work vectors */ 278 ierr = VecDestroy(&stage_slow);CHKERRQ(ierr); 279 ierr = VecDestroy(&sol_slow);CHKERRQ(ierr); 280 PetscFunctionReturn(0); 281 } 282 283 /* 284 This is for partitioned RHS multirate RK method 285 The step completion formula is 286 287 x1 = x0 + h b^T YdotRHS 288 289 */ 290 static PetscErrorCode TSEvaluateStep_MRKSPLIT(TS ts,PetscInt order,Vec X,PetscBool *done) 291 { 292 TS_RK *rk = (TS_RK*)ts->data; 293 RKTableau tab = rk->tableau; 294 Vec Xslow,Xfast; /* subvectors of X which store slow components and fast components respectively */ 295 PetscScalar *w = rk->work; 296 PetscReal h = ts->time_step; 297 PetscInt s = tab->s,j; 298 PetscErrorCode ierr; 299 300 PetscFunctionBegin; 301 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 302 if (rk->slow) { 303 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 304 ierr = VecGetSubVector(ts->vec_sol,rk->is_slow,&Xslow);CHKERRQ(ierr); 305 ierr = VecMAXPY(Xslow,s,w,((TS_RK*)rk->subts_slow->data)->YdotRHS);CHKERRQ(ierr); 306 ierr = VecRestoreSubVector(ts->vec_sol,rk->is_slow,&Xslow);CHKERRQ(ierr);; 307 } else { 308 for (j=0; j<s; j++) w[j] = h/rk->dtratio*tab->b[j]; 309 ierr = VecGetSubVector(X,rk->is_fast,&Xfast);CHKERRQ(ierr); 310 ierr = VecMAXPY(Xfast,s,w,((TS_RK*)rk->subts_fast->data)->YdotRHS);CHKERRQ(ierr); 311 ierr = VecRestoreSubVector(X,rk->is_fast,&Xfast);CHKERRQ(ierr); 312 } 313 PetscFunctionReturn(0); 314 } 315 316 static PetscErrorCode TSStep_MRKSPLIT(TS ts) 317 { 318 TS_RK *rk = (TS_RK*)ts->data; 319 RKTableau tab = rk->tableau; 320 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 321 Vec *YdotRHS_fast,*YdotRHS_slow; 322 Vec Yslow,Yfast; /* subvectors store the stges of slow components and fast components respectively */ 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 next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 328 PetscErrorCode ierr; 329 330 PetscFunctionBegin; 331 rk->status = TS_STEP_INCOMPLETE; 332 YdotRHS_fast = ((TS_RK*)rk->subts_fast->data)->YdotRHS; 333 YdotRHS_slow = ((TS_RK*)rk->subts_slow->data)->YdotRHS; 334 for (i=0; i<s; i++) { 335 ierr = VecGetSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i]);CHKERRQ(ierr); 336 ierr = VecGetSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);CHKERRQ(ierr); 337 } 338 ierr = VecCopy(ts->vec_sol,rk->X0);CHKERRQ(ierr); 339 /* propogate both slow and fast components using large time steps */ 340 for (i=0; i<s; i++) { 341 rk->stage_time = t + h*c[i]; 342 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 343 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 344 /*ierr = VecCopy(ts->vec_sol,Yc[i]);CHKERRQ(ierr);*/ 345 ierr = VecGetSubVector(Y[i],rk->is_fast,&Yfast);CHKERRQ(ierr); 346 ierr = VecGetSubVector(Y[i],rk->is_slow,&Yslow);CHKERRQ(ierr); 347 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 348 ierr = VecMAXPY(Yfast,i,w,YdotRHS_fast);CHKERRQ(ierr); 349 ierr = VecMAXPY(Yslow,i,w,YdotRHS_slow);CHKERRQ(ierr); 350 ierr = VecRestoreSubVector(Y[i],rk->is_fast,&Yfast);CHKERRQ(ierr); 351 ierr = VecRestoreSubVector(Y[i],rk->is_slow,&Yslow);CHKERRQ(ierr); 352 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 353 ierr = TSComputeRHSFunction(rk->subts_slow,t+h*c[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 354 ierr = TSComputeRHSFunction(rk->subts_fast,t+h*c[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 355 } 356 rk->slow = PETSC_TRUE; 357 /* update the value of slow components when using slow time step */ 358 ierr = TSEvaluateStep_MRKSPLIT(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 359 360 if (((TS_RK*)rk->subts_fast->data)->subts_fast) { 361 rk->subts_fast->ptime = t; 362 rk->subts_fast->time_step = h/rk->dtratio; 363 for (k=0; k<rk->dtratio; k++) { 364 ierr = TSStep_MRKSPLIT(rk->subts_fast);CHKERRQ(ierr); 365 } 366 } else { 367 for (k=0; k<rk->dtratio; k++) { 368 /* propogate fast component using small time steps */ 369 for (i=0; i<s; i++) { 370 rk->stage_time = t + h/rk->dtratio*c[i]; 371 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 372 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 373 /* stage value for fast components */ 374 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 375 ierr = VecGetSubVector(Y[i],rk->is_fast,&Yfast);CHKERRQ(ierr); 376 ierr = VecMAXPY(Yfast,i,w,YdotRHS_fast);CHKERRQ(ierr); 377 ierr = VecRestoreSubVector(Y[i],rk->is_fast,&Yfast);CHKERRQ(ierr); 378 /* stage value for slow components */ 379 ierr = VecGetSubVector(Y[i],rk->is_slow,&Yslow);CHKERRQ(ierr); 380 ierr = TSInterpolate_MRKSPLIT(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Yslow);CHKERRQ(ierr); 381 ierr = VecRestoreSubVector(Y[i],rk->is_slow,&Yslow);CHKERRQ(ierr); 382 ierr = TSPostStage(ts,rk->stage_time,i,Y);CHKERRQ(ierr); 383 /* compute the stage RHS for fast components */ 384 ierr = TSComputeRHSFunction(rk->subts_fast,t+k*h*rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 385 } 386 /* update the value of fast components when using fast time step */ 387 rk->slow = PETSC_FALSE; 388 ierr = TSEvaluateStep_MRKSPLIT(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 389 } 390 } 391 for (i=0; i<s; i++) { 392 ierr = VecRestoreSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i]);CHKERRQ(ierr); 393 ierr = VecRestoreSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);CHKERRQ(ierr); 394 } 395 ts->ptime = t + ts->time_step; 396 ts->time_step = next_time_step; 397 rk->status = TS_STEP_COMPLETE; 398 PetscFunctionReturn(0); 399 } 400 401 /*@C 402 TSRKSetMultirateType - Set the type of RK Multirate scheme 403 404 Logically collective 405 406 Input Parameter: 407 + ts - timestepping context 408 - mrktype - type of MRK-scheme 409 410 Options Database: 411 . -ts_rk_multiarte_type - <none,nonsplit,split> 412 413 Level: intermediate 414 @*/ 415 PetscErrorCode TSRKSetMultirateType(TS ts, TSMRKType mrktype) 416 { 417 TS_RK *rk = (TS_RK*)ts->data; 418 PetscErrorCode ierr; 419 420 PetscFunctionBegin; 421 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 422 switch(mrktype){ 423 case TSMRKNONE: 424 break; 425 case TSMRKNONSPLIT: 426 ts->ops->step = TSStep_MRKNONSPLIT; 427 ts->ops->interpolate = TSInterpolate_MRKNONSPLIT; 428 rk->dtratio = 2; 429 ierr = TSRKSetType(ts,TSMRKDefault);CHKERRQ(ierr); 430 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_MRKNONSPLIT_C",TSSetUp_MRKNONSPLIT);CHKERRQ(ierr); 431 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSReset_MRKNONSPLIT_C",TSReset_MRKNONSPLIT);CHKERRQ(ierr); 432 break; 433 case TSMRKSPLIT: 434 ts->ops->step = TSStep_MRKSPLIT; 435 ts->ops->evaluatestep = TSEvaluateStep_MRKSPLIT; 436 ts->ops->interpolate = TSInterpolate_MRKSPLIT; 437 rk->dtratio = 2; 438 ierr = TSRKSetType(ts,TSMRKDefault);CHKERRQ(ierr); 439 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_MRKSPLIT_C",TSSetUp_MRKSPLIT);CHKERRQ(ierr); 440 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSReset_MRKSPLIT_C",TSReset_MRKSPLIT);CHKERRQ(ierr); 441 break; 442 default : 443 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown type '%s'",mrktype); 444 } 445 PetscFunctionReturn(0); 446 } 447 448