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

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

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

     65:   Not Collective

     67:   Level: advanced

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

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

     86:   Level: developer

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

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

    103:   Level: developer

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

    118: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt, TSGLLEAdaptType type)
    119: {
    120:   PetscErrorCode (*r)(TSGLLEAdapt);

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

    131: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt, const char prefix[])
    132: {
    133:   PetscFunctionBegin;
    134:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adapt, prefix));
    135:   PetscFunctionReturn(PETSC_SUCCESS);
    136: }

    138: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt, PetscViewer viewer)
    139: {
    140:   PetscBool isascii;

    142:   PetscFunctionBegin;
    143:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    144:   if (isascii) {
    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: }

    155: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
    156: {
    157:   PetscFunctionBegin;
    158:   if (!*adapt) PetscFunctionReturn(PETSC_SUCCESS);
    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: }

    169: PetscErrorCode TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt, PetscOptionItems PetscOptionsObject)
    170: {
    171:   char      type[256] = TSGLLEADAPT_BOTH;
    172:   PetscBool flg;

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

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

    199: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm, TSGLLEAdapt *inadapt)
    200: {
    201:   TSGLLEAdapt adapt;

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

    210: /*
    211:    Implementations
    212: */

    214: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
    215: {
    216:   PetscFunctionBegin;
    217:   PetscCall(PetscFree(adapt->data));
    218:   PetscFunctionReturn(PETSC_SUCCESS);
    219: }

    221: /* -------------------------------- None ----------------------------------- */
    222: typedef struct {
    223:   PetscInt  scheme;
    224:   PetscReal h;
    225: } TSGLLEAdapt_None;

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

    239: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
    240: {
    241:   TSGLLEAdapt_None *a;

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

    251: /* -------------------------------- Size ----------------------------------- */
    252: typedef struct {
    253:   PetscReal desired_h;
    254: } TSGLLEAdapt_Size;

    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;

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

    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;

    273:   if (*next_h > tleft) {
    274:     *finish = PETSC_TRUE;
    275:     *next_h = tleft;
    276:   } else *finish = PETSC_FALSE;
    277:   PetscFunctionReturn(PETSC_SUCCESS);
    278: }

    280: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
    281: {
    282:   TSGLLEAdapt_Size *a;

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

    292: /* -------------------------------- Both ----------------------------------- */
    293: typedef struct {
    294:   PetscInt  count_at_order;
    295:   PetscReal desired_h;
    296: } TSGLLEAdapt_Both;

    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;

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

    335:   if (*next_h > tleft) {
    336:     *finish = PETSC_TRUE;
    337:     *next_h = tleft;
    338:   } else *finish = PETSC_FALSE;
    339:   PetscFunctionReturn(PETSC_SUCCESS);
    340: }

    342: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
    343: {
    344:   TSGLLEAdapt_Both *a;

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