1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscdm.h> 3 static PetscErrorCode TSRHSSplitGetRHSSplit(TS ts,const char splitname[],TS_RHSSplitLink *isplit) 4 { 5 PetscBool found = PETSC_FALSE; 6 PetscErrorCode ierr; 7 8 PetscFunctionBegin; 9 *isplit = ts->tsrhssplit; 10 /* look up the split */ 11 while (*isplit) { 12 ierr = PetscStrcmp((*isplit)->splitname,splitname,&found);CHKERRQ(ierr); 13 if (found) break; 14 *isplit = (*isplit)->next; 15 } 16 PetscFunctionReturn(0); 17 } 18 19 /*@C 20 TSRHSSplitSetIS - Set the index set for the specified split 21 22 Logically Collective on TS 23 24 Input Parameters: 25 + ts - the TS context obtained from TSCreate() 26 . splitname - name of this split, if NULL the number of the split is used 27 - is - the index set for part of the solution vector 28 29 Level: intermediate 30 31 .seealso: TSRHSSplitGetIS() 32 33 .keywords: TS, TSRHSSplit 34 @*/ 35 PetscErrorCode TSRHSSplitSetIS(TS ts,const char splitname[],IS is) 36 { 37 TS_RHSSplitLink newsplit,next = ts->tsrhssplit; 38 char prefix[128]; 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 43 PetscValidHeaderSpecific(is,IS_CLASSID,3); 44 45 ierr = PetscNew(&newsplit);CHKERRQ(ierr); 46 if (splitname) { 47 ierr = PetscStrallocpy(splitname,&newsplit->splitname);CHKERRQ(ierr); 48 } else { 49 ierr = PetscMalloc1(8,&newsplit->splitname);CHKERRQ(ierr); 50 ierr = PetscSNPrintf(newsplit->splitname,7,"%D",ts->num_rhs_splits);CHKERRQ(ierr); 51 } 52 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 53 newsplit->is = is; 54 ierr = TSCreate(PetscObjectComm((PetscObject)ts),&newsplit->ts);CHKERRQ(ierr); 55 ierr = PetscObjectIncrementTabLevel((PetscObject)newsplit->ts,(PetscObject)ts,1);CHKERRQ(ierr); 56 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)newsplit->ts);CHKERRQ(ierr); 57 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%srhsplit_%s_",((PetscObject)ts)->prefix ? ((PetscObject)ts)->prefix : "",newsplit->splitname);CHKERRQ(ierr); 58 ierr = TSSetOptionsPrefix(newsplit->ts,prefix);CHKERRQ(ierr); 59 if (!next) ts->tsrhssplit = newsplit; 60 else { 61 while (next->next) next = next->next; 62 next->next = newsplit; 63 } 64 ts->num_rhs_splits++; 65 PetscFunctionReturn(0); 66 } 67 68 /*@C 69 TSRHSSplitGetIS - Retrieves the elements for a split as an IS 70 71 Logically Collective on TS 72 73 Input Parameters: 74 + ts - the TS context obtained from TSCreate() 75 - splitname - name of this split 76 77 Output Parameters: 78 - is - the index set for part of the solution vector 79 80 Level: intermediate 81 82 .seealso: TSRHSSplitSetIS() 83 84 .keywords: TS, TSRHSSplit 85 @*/ 86 PetscErrorCode TSRHSSplitGetIS(TS ts,const char splitname[],IS *is) 87 { 88 TS_RHSSplitLink isplit; 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 93 *is = NULL; 94 /* look up the split */ 95 ierr = TSRHSSplitGetRHSSplit(ts,splitname,&isplit);CHKERRQ(ierr); 96 if (isplit) *is = isplit->is; 97 PetscFunctionReturn(0); 98 } 99 100 /*@C 101 TSRHSSplitSetRHSFunction - Set the split right-hand-side functions. 102 103 Logically Collective on TS 104 105 Input Parameters: 106 + ts - the TS context obtained from TSCreate() 107 . splitname - name of this split 108 . r - vector to hold the residual (or NULL to have it created internally) 109 . rhsfunc - the RHS function evaluation routine 110 - ctx - user-defined context for private data for the split function evaluation routine (may be NULL) 111 112 Calling sequence of fun: 113 $ rhsfunc(TS ts,PetscReal t,Vec u,Vec f,ctx); 114 115 + t - time at step/stage being solved 116 . u - state vector 117 . f - function vector 118 - ctx - [optional] user-defined context for matrix evaluation routine (may be NULL) 119 120 Level: beginner 121 122 .keywords: TS, timestep, set, ODE, Hamiltonian, Function 123 @*/ 124 PetscErrorCode TSRHSSplitSetRHSFunction(TS ts,const char splitname[],Vec r,TSRHSFunction rhsfunc,void *ctx) 125 { 126 TS_RHSSplitLink isplit; 127 DM dmc; 128 Vec subvec,ralloc = NULL; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 133 if (r) PetscValidHeaderSpecific(r,VEC_CLASSID,2); 134 135 /* look up the split */ 136 ierr = TSRHSSplitGetRHSSplit(ts,splitname,&isplit);CHKERRQ(ierr); 137 if (!isplit) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"The split %s is not created, check the split name or call TSRHSSplitSetIS() to create one",splitname); 138 139 if (!r && ts->vec_sol) { 140 ierr = VecGetSubVector(ts->vec_sol,isplit->is,&subvec);CHKERRQ(ierr); 141 ierr = VecDuplicate(subvec,&ralloc);CHKERRQ(ierr); 142 r = ralloc; 143 ierr = VecRestoreSubVector(ts->vec_sol,isplit->is,&subvec);CHKERRQ(ierr); 144 } 145 146 if(ts->dm){ 147 ierr = DMClone(ts->dm, &dmc);CHKERRQ(ierr); 148 ierr = TSSetDM(isplit->ts, dmc);CHKERRQ(ierr); 149 ierr = DMDestroy(&dmc);CHKERRQ(ierr); 150 } 151 152 ierr = TSSetRHSFunction(isplit->ts,r,rhsfunc,ctx);CHKERRQ(ierr); 153 ierr = VecDestroy(&ralloc);CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 157 /*@C 158 TSRHSSplitGetSubTS - Get the sub-TS by split name. 159 160 Logically Collective on TS 161 162 Output Parameters: 163 + splitname - the number of the split 164 - subts - the array of TS contexts 165 166 Level: advanced 167 168 .seealso: TSGetRHSSplitFunction() 169 @*/ 170 PetscErrorCode TSRHSSplitGetSubTS(TS ts,const char splitname[],TS *subts) 171 { 172 TS_RHSSplitLink isplit; 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 177 PetscValidPointer(subts,3); 178 *subts = NULL; 179 /* look up the split */ 180 ierr = TSRHSSplitGetRHSSplit(ts,splitname,&isplit);CHKERRQ(ierr); 181 if (isplit) *subts = isplit->ts; 182 PetscFunctionReturn(0); 183 } 184 185 /*@C 186 TSRHSSplitGetSubTSs - Get an array of all sub-TS contexts. 187 188 Logically Collective on TS 189 190 Output Parameters: 191 + n - the number of splits 192 - subksp - the array of TS contexts 193 194 Note: 195 After TSRHSSplitGetSubTS() the array of TSs is to be freed by the user with PetscFree() 196 (not the TS just the array that contains them). 197 198 Level: advanced 199 200 .seealso: TSGetRHSSplitFunction() 201 @*/ 202 PetscErrorCode TSRHSSplitGetSubTSs(TS ts,PetscInt *n,TS *subts[]) 203 { 204 TS_RHSSplitLink ilink = ts->tsrhssplit; 205 PetscInt i = 0; 206 PetscErrorCode ierr; 207 208 PetscFunctionBegin; 209 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 210 if (subts) { 211 ierr = PetscMalloc1(ts->num_rhs_splits,subts);CHKERRQ(ierr); 212 while (ilink) { 213 (*subts)[i++] = ilink->ts; 214 ilink = ilink->next; 215 } 216 } 217 if (n) *n = ts->num_rhs_splits; 218 PetscFunctionReturn(0); 219 } 220