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