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