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