1 #include <../src/ts/impls/implicit/glle/glle.h> /*I "petscts.h" I*/ 2 3 static PetscFunctionList TSGLLEAdaptList; 4 static PetscBool TSGLLEAdaptPackageInitialized; 5 static PetscBool TSGLLEAdaptRegisterAllCalled; 6 static PetscClassId TSGLLEADAPT_CLASSID; 7 8 struct _TSGLLEAdaptOps { 9 PetscErrorCode (*choose)(TSGLLEAdapt, PetscInt, const PetscInt[], const PetscReal[], const PetscReal[], PetscInt, PetscReal, PetscReal, PetscInt *, PetscReal *, PetscBool *); 10 PetscErrorCode (*destroy)(TSGLLEAdapt); 11 PetscErrorCode (*view)(TSGLLEAdapt, PetscViewer); 12 PetscErrorCode (*setfromoptions)(TSGLLEAdapt, PetscOptionItems); 13 }; 14 15 struct _p_TSGLLEAdapt { 16 PETSCHEADER(struct _TSGLLEAdaptOps); 17 void *data; 18 }; 19 20 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt); 21 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt); 22 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt); 23 24 /*@C 25 TSGLLEAdaptRegister - adds a `TSGLLEAdapt` implementation 26 27 Not Collective, No Fortran Support 28 29 Input Parameters: 30 + sname - name of user-defined adaptivity scheme 31 - function - routine to create method context 32 33 Level: advanced 34 35 Note: 36 `TSGLLEAdaptRegister()` may be called multiple times to add several user-defined families. 37 38 Example Usage: 39 .vb 40 TSGLLEAdaptRegister("my_scheme", MySchemeCreate); 41 .ve 42 43 Then, your scheme can be chosen with the procedural interface via 44 $ TSGLLEAdaptSetType(ts, "my_scheme") 45 or at runtime via the option 46 $ -ts_adapt_type my_scheme 47 48 .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegisterAll()` 49 @*/ 50 PetscErrorCode TSGLLEAdaptRegister(const char sname[], PetscErrorCode (*function)(TSGLLEAdapt)) 51 { 52 PetscFunctionBegin; 53 PetscCall(TSGLLEAdaptInitializePackage()); 54 PetscCall(PetscFunctionListAdd(&TSGLLEAdaptList, sname, function)); 55 PetscFunctionReturn(PETSC_SUCCESS); 56 } 57 58 /*@C 59 TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in `TSGLLEAdapt` 60 61 Not Collective 62 63 Level: advanced 64 65 .seealso: [](ch_ts), `TSGLLEAdapt`, `TSGLLE`, `TSGLLEAdaptRegisterDestroy()` 66 @*/ 67 PetscErrorCode TSGLLEAdaptRegisterAll(void) 68 { 69 PetscFunctionBegin; 70 if (TSGLLEAdaptRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 71 TSGLLEAdaptRegisterAllCalled = PETSC_TRUE; 72 PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_NONE, TSGLLEAdaptCreate_None)); 73 PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_SIZE, TSGLLEAdaptCreate_Size)); 74 PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_BOTH, TSGLLEAdaptCreate_Both)); 75 PetscFunctionReturn(PETSC_SUCCESS); 76 } 77 78 /*@C 79 TSGLLEAdaptFinalizePackage - This function destroys everything in the `TSGLLE` package. It is 80 called from `PetscFinalize()`. 81 82 Level: developer 83 84 .seealso: [](ch_ts), `PetscFinalize()`, `TSGLLEAdapt`, `TSGLLEAdaptInitializePackage()` 85 @*/ 86 PetscErrorCode TSGLLEAdaptFinalizePackage(void) 87 { 88 PetscFunctionBegin; 89 PetscCall(PetscFunctionListDestroy(&TSGLLEAdaptList)); 90 TSGLLEAdaptPackageInitialized = PETSC_FALSE; 91 TSGLLEAdaptRegisterAllCalled = PETSC_FALSE; 92 PetscFunctionReturn(PETSC_SUCCESS); 93 } 94 95 /*@C 96 TSGLLEAdaptInitializePackage - This function initializes everything in the `TSGLLEAdapt` package. It is 97 called from `TSInitializePackage()`. 98 99 Level: developer 100 101 .seealso: [](ch_ts), `PetscInitialize()`, `TSGLLEAdapt`, `TSGLLEAdaptFinalizePackage()` 102 @*/ 103 PetscErrorCode TSGLLEAdaptInitializePackage(void) 104 { 105 PetscFunctionBegin; 106 if (TSGLLEAdaptPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 107 TSGLLEAdaptPackageInitialized = PETSC_TRUE; 108 PetscCall(PetscClassIdRegister("TSGLLEAdapt", &TSGLLEADAPT_CLASSID)); 109 PetscCall(TSGLLEAdaptRegisterAll()); 110 PetscCall(PetscRegisterFinalize(TSGLLEAdaptFinalizePackage)); 111 PetscFunctionReturn(PETSC_SUCCESS); 112 } 113 114 PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt, TSGLLEAdaptType type) 115 { 116 PetscErrorCode (*r)(TSGLLEAdapt); 117 118 PetscFunctionBegin; 119 PetscCall(PetscFunctionListFind(TSGLLEAdaptList, type, &r)); 120 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLEAdapt type \"%s\" given", type); 121 if (((PetscObject)adapt)->type_name) PetscUseTypeMethod(adapt, destroy); 122 PetscCall((*r)(adapt)); 123 PetscCall(PetscObjectChangeTypeName((PetscObject)adapt, type)); 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt, const char prefix[]) 128 { 129 PetscFunctionBegin; 130 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adapt, prefix)); 131 PetscFunctionReturn(PETSC_SUCCESS); 132 } 133 134 PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt, PetscViewer viewer) 135 { 136 PetscBool iascii; 137 138 PetscFunctionBegin; 139 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 140 if (iascii) { 141 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adapt, viewer)); 142 if (adapt->ops->view) { 143 PetscCall(PetscViewerASCIIPushTab(viewer)); 144 PetscUseTypeMethod(adapt, view, viewer); 145 PetscCall(PetscViewerASCIIPopTab(viewer)); 146 } 147 } 148 PetscFunctionReturn(PETSC_SUCCESS); 149 } 150 151 PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt) 152 { 153 PetscFunctionBegin; 154 if (!*adapt) PetscFunctionReturn(PETSC_SUCCESS); 155 PetscValidHeaderSpecific(*adapt, TSGLLEADAPT_CLASSID, 1); 156 if (--((PetscObject)*adapt)->refct > 0) { 157 *adapt = NULL; 158 PetscFunctionReturn(PETSC_SUCCESS); 159 } 160 PetscTryTypeMethod(*adapt, destroy); 161 PetscCall(PetscHeaderDestroy(adapt)); 162 PetscFunctionReturn(PETSC_SUCCESS); 163 } 164 165 PetscErrorCode TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt, PetscOptionItems PetscOptionsObject) 166 { 167 char type[256] = TSGLLEADAPT_BOTH; 168 PetscBool flg; 169 170 PetscFunctionBegin; 171 /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this 172 * function can only be called from inside TSSetFromOptions_GLLE() */ 173 PetscOptionsHeadBegin(PetscOptionsObject, "TSGLLE Adaptivity options"); 174 PetscCall(PetscOptionsFList("-ts_adapt_type", "Algorithm to use for adaptivity", "TSGLLEAdaptSetType", TSGLLEAdaptList, ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type, type, sizeof(type), &flg)); 175 if (flg || !((PetscObject)adapt)->type_name) PetscCall(TSGLLEAdaptSetType(adapt, type)); 176 PetscTryTypeMethod(adapt, setfromoptions, PetscOptionsObject); 177 PetscOptionsHeadEnd(); 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180 181 PetscErrorCode TSGLLEAdaptChoose(TSGLLEAdapt adapt, PetscInt n, const PetscInt orders[], const PetscReal errors[], const PetscReal cost[], PetscInt cur, PetscReal h, PetscReal tleft, PetscInt *next_sc, PetscReal *next_h, PetscBool *finish) 182 { 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(adapt, TSGLLEADAPT_CLASSID, 1); 185 PetscAssertPointer(orders, 3); 186 PetscAssertPointer(errors, 4); 187 PetscAssertPointer(cost, 5); 188 PetscAssertPointer(next_sc, 9); 189 PetscAssertPointer(next_h, 10); 190 PetscAssertPointer(finish, 11); 191 PetscUseTypeMethod(adapt, choose, n, orders, errors, cost, cur, h, tleft, next_sc, next_h, finish); 192 PetscFunctionReturn(PETSC_SUCCESS); 193 } 194 195 PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm, TSGLLEAdapt *inadapt) 196 { 197 TSGLLEAdapt adapt; 198 199 PetscFunctionBegin; 200 *inadapt = NULL; 201 PetscCall(PetscHeaderCreate(adapt, TSGLLEADAPT_CLASSID, "TSGLLEAdapt", "General Linear adaptivity", "TS", comm, TSGLLEAdaptDestroy, TSGLLEAdaptView)); 202 *inadapt = adapt; 203 PetscFunctionReturn(PETSC_SUCCESS); 204 } 205 206 /* 207 Implementations 208 */ 209 210 static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt) 211 { 212 PetscFunctionBegin; 213 PetscCall(PetscFree(adapt->data)); 214 PetscFunctionReturn(PETSC_SUCCESS); 215 } 216 217 /* -------------------------------- None ----------------------------------- */ 218 typedef struct { 219 PetscInt scheme; 220 PetscReal h; 221 } TSGLLEAdapt_None; 222 223 static PetscErrorCode TSGLLEAdaptChoose_None(TSGLLEAdapt adapt, PetscInt n, const PetscInt orders[], const PetscReal errors[], const PetscReal cost[], PetscInt cur, PetscReal h, PetscReal tleft, PetscInt *next_sc, PetscReal *next_h, PetscBool *finish) 224 { 225 PetscFunctionBegin; 226 *next_sc = cur; 227 *next_h = h; 228 if (*next_h > tleft) { 229 *finish = PETSC_TRUE; 230 *next_h = tleft; 231 } else *finish = PETSC_FALSE; 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 235 PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt) 236 { 237 TSGLLEAdapt_None *a; 238 239 PetscFunctionBegin; 240 PetscCall(PetscNew(&a)); 241 adapt->data = (void *)a; 242 adapt->ops->choose = TSGLLEAdaptChoose_None; 243 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 244 PetscFunctionReturn(PETSC_SUCCESS); 245 } 246 247 /* -------------------------------- Size ----------------------------------- */ 248 typedef struct { 249 PetscReal desired_h; 250 } TSGLLEAdapt_Size; 251 252 static PetscErrorCode TSGLLEAdaptChoose_Size(TSGLLEAdapt adapt, PetscInt n, const PetscInt orders[], const PetscReal errors[], const PetscReal cost[], PetscInt cur, PetscReal h, PetscReal tleft, PetscInt *next_sc, PetscReal *next_h, PetscBool *finish) 253 { 254 TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size *)adapt->data; 255 PetscReal dec = 0.2, inc = 5.0, safe = 0.9, optimal, last_desired_h; 256 257 PetscFunctionBegin; 258 *next_sc = cur; 259 optimal = PetscPowReal((PetscReal)errors[cur], (PetscReal)-1. / (safe * orders[cur])); 260 /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the 261 * one that would have been taken (without smoothing) on the last step. */ 262 last_desired_h = sz->desired_h; 263 sz->desired_h = h * PetscMax(dec, PetscMin(inc, optimal)); /* Trim to [dec,inc] */ 264 265 /* Normally only happens on the first step */ 266 if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h); 267 else *next_h = sz->desired_h; 268 269 if (*next_h > tleft) { 270 *finish = PETSC_TRUE; 271 *next_h = tleft; 272 } else *finish = PETSC_FALSE; 273 PetscFunctionReturn(PETSC_SUCCESS); 274 } 275 276 PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt) 277 { 278 TSGLLEAdapt_Size *a; 279 280 PetscFunctionBegin; 281 PetscCall(PetscNew(&a)); 282 adapt->data = (void *)a; 283 adapt->ops->choose = TSGLLEAdaptChoose_Size; 284 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 285 PetscFunctionReturn(PETSC_SUCCESS); 286 } 287 288 /* -------------------------------- Both ----------------------------------- */ 289 typedef struct { 290 PetscInt count_at_order; 291 PetscReal desired_h; 292 } TSGLLEAdapt_Both; 293 294 static PetscErrorCode TSGLLEAdaptChoose_Both(TSGLLEAdapt adapt, PetscInt n, const PetscInt orders[], const PetscReal errors[], const PetscReal cost[], PetscInt cur, PetscReal h, PetscReal tleft, PetscInt *next_sc, PetscReal *next_h, PetscBool *finish) 295 { 296 TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both *)adapt->data; 297 PetscReal dec = 0.2, inc = 5.0, safe = 0.9; 298 struct { 299 PetscInt id; 300 PetscReal h, eff; 301 } best = {-1, 0, 0}, trial = {-1, 0, 0}, current = {-1, 0, 0}; 302 PetscInt i; 303 304 PetscFunctionBegin; 305 for (i = 0; i < n; i++) { 306 PetscReal optimal; 307 trial.id = i; 308 optimal = PetscPowReal((PetscReal)errors[i], (PetscReal)-1. / (safe * orders[i])); 309 trial.h = h * optimal; 310 trial.eff = trial.h / cost[i]; 311 if (trial.eff > best.eff) PetscCall(PetscArraycpy(&best, &trial, 1)); 312 if (i == cur) PetscCall(PetscArraycpy(¤t, &trial, 1)); 313 } 314 /* Only switch orders if the scheme offers significant benefits over the current one. 315 When the scheme is not changing, only change step size if it offers significant benefits. */ 316 if (best.eff < 1.2 * current.eff || both->count_at_order < orders[cur] + 2) { 317 PetscReal last_desired_h; 318 *next_sc = current.id; 319 last_desired_h = both->desired_h; 320 both->desired_h = PetscMax(h * dec, PetscMin(h * inc, current.h)); 321 *next_h = (both->count_at_order > 0) ? PetscSqrtReal(last_desired_h * both->desired_h) : both->desired_h; 322 both->count_at_order++; 323 } else { 324 PetscReal rat = cost[best.id] / cost[cur]; 325 *next_sc = best.id; 326 *next_h = PetscMax(h * rat * dec, PetscMin(h * rat * inc, best.h)); 327 both->count_at_order = 0; 328 both->desired_h = best.h; 329 } 330 331 if (*next_h > tleft) { 332 *finish = PETSC_TRUE; 333 *next_h = tleft; 334 } else *finish = PETSC_FALSE; 335 PetscFunctionReturn(PETSC_SUCCESS); 336 } 337 338 PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt) 339 { 340 TSGLLEAdapt_Both *a; 341 342 PetscFunctionBegin; 343 PetscCall(PetscNew(&a)); 344 adapt->data = (void *)a; 345 adapt->ops->choose = TSGLLEAdaptChoose_Both; 346 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 347 PetscFunctionReturn(PETSC_SUCCESS); 348 } 349