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