Actual source code: glleadapt.c

  1: #include <../src/ts/impls/implicit/glle/glle.h>

  3: static PetscFunctionList TSGLLEAdaptList;
  4: static PetscBool         TSGLLEAdaptPackageInitialized;
  5: static PetscBool         TSGLLEAdaptRegisterAllCalled;
  6: static PetscClassId      TSGLLEADAPT_CLASSID;

  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: };

 15: struct _p_TSGLLEAdapt {
 16:   PETSCHEADER(struct _TSGLLEAdaptOps);
 17:   void *data;
 18: };

 20: PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt);
 21: PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt);
 22: PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt);

 24: /*@C
 25:   TSGLLEAdaptRegister -  adds a `TSGLLEAdapt` implementation

 27:   Not Collective, No Fortran Support

 29:   Input Parameters:
 30: + sname    - name of user-defined adaptivity scheme
 31: - function - routine to create method context

 33:   Level: advanced

 35:   Note:
 36:   `TSGLLEAdaptRegister()` may be called multiple times to add several user-defined families.

 38:   Example Usage:
 39: .vb
 40:    TSGLLEAdaptRegister("my_scheme", MySchemeCreate);
 41: .ve

 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

 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: }

 58: /*@C
 59:   TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in `TSGLLEAdapt`

 61:   Not Collective

 63:   Level: advanced

 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: }

 78: /*@C
 79:   TSGLLEAdaptFinalizePackage - This function destroys everything in the `TSGLLE` package. It is
 80:   called from `PetscFinalize()`.

 82:   Level: developer

 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: }

 95: /*@C
 96:   TSGLLEAdaptInitializePackage - This function initializes everything in the `TSGLLEAdapt` package. It is
 97:   called from `TSInitializePackage()`.

 99:   Level: developer

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: }

114: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt, TSGLLEAdaptType type)
115: {
116:   PetscErrorCode (*r)(TSGLLEAdapt);

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: }

127: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt, const char prefix[])
128: {
129:   PetscFunctionBegin;
130:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adapt, prefix));
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }

134: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt, PetscViewer viewer)
135: {
136:   PetscBool iascii;

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: }

151: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
152: {
153:   PetscFunctionBegin;
154:   if (!*adapt) PetscFunctionReturn(PETSC_SUCCESS);
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: }

165: PetscErrorCode TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt, PetscOptionItems *PetscOptionsObject)
166: {
167:   char      type[256] = TSGLLEADAPT_BOTH;
168:   PetscBool flg;

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: }

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;
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: }

195: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm, TSGLLEAdapt *inadapt)
196: {
197:   TSGLLEAdapt adapt;

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: }

206: /*
207:    Implementations
208: */

210: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
211: {
212:   PetscFunctionBegin;
213:   PetscCall(PetscFree(adapt->data));
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: /* -------------------------------- None ----------------------------------- */
218: typedef struct {
219:   PetscInt  scheme;
220:   PetscReal h;
221: } TSGLLEAdapt_None;

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: }

235: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
236: {
237:   TSGLLEAdapt_None *a;

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: }

247: /* -------------------------------- Size ----------------------------------- */
248: typedef struct {
249:   PetscReal desired_h;
250: } TSGLLEAdapt_Size;

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;

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] */

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;

269:   if (*next_h > tleft) {
270:     *finish = PETSC_TRUE;
271:     *next_h = tleft;
272:   } else *finish = PETSC_FALSE;
273:   PetscFunctionReturn(PETSC_SUCCESS);
274: }

276: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
277: {
278:   TSGLLEAdapt_Size *a;

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: }

288: /* -------------------------------- Both ----------------------------------- */
289: typedef struct {
290:   PetscInt  count_at_order;
291:   PetscReal desired_h;
292: } TSGLLEAdapt_Both;

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;

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(&current, &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:   }

331:   if (*next_h > tleft) {
332:     *finish = PETSC_TRUE;
333:     *next_h = tleft;
334:   } else *finish = PETSC_FALSE;
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
339: {
340:   TSGLLEAdapt_Both *a;

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: }