Actual source code: glleadapt.c

petsc-3.13.6 2020-09-29
Report Typos and Errors

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

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

  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)(PetscOptionItems*,TSGLLEAdapt);
 14: };

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

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

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

 28:    Not Collective

 30:    Input Parameters:
 31: +  name_scheme - name of user-defined adaptivity scheme
 32: -  routine_create - routine to create method context

 34:    Notes:
 35:    TSGLLEAdaptRegister() may be called multiple times to add several user-defined families.

 37:    Sample usage:
 38: .vb
 39:    TSGLLEAdaptRegister("my_scheme",MySchemeCreate);
 40: .ve

 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

 47:    Level: advanced

 49: .seealso: TSGLLEAdaptRegisterAll()
 50: @*/
 51: PetscErrorCode  TSGLLEAdaptRegister(const char sname[],PetscErrorCode (*function)(TSGLLEAdapt))
 52: {

 56:   TSGLLEAdaptInitializePackage();
 57:   PetscFunctionListAdd(&TSGLLEAdaptList,sname,function);
 58:   return(0);
 59: }

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

 64:   Not Collective

 66:   Level: advanced

 68: .seealso: TSGLLEAdaptRegisterDestroy()
 69: @*/
 70: PetscErrorCode  TSGLLEAdaptRegisterAll(void)
 71: {

 75:   if (TSGLLEAdaptRegisterAllCalled) return(0);
 76:   TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
 77:   TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None);
 78:   TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size);
 79:   TSGLLEAdaptRegister(TSGLLEADAPT_BOTH,TSGLLEAdaptCreate_Both);
 80:   return(0);
 81: }

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

 87:   Level: developer

 89: .seealso: PetscFinalize()
 90: @*/
 91: PetscErrorCode  TSGLLEAdaptFinalizePackage(void)
 92: {

 96:   PetscFunctionListDestroy(&TSGLLEAdaptList);
 97:   TSGLLEAdaptPackageInitialized = PETSC_FALSE;
 98:   TSGLLEAdaptRegisterAllCalled  = PETSC_FALSE;
 99:   return(0);
100: }

102: /*@C
103:   TSGLLEAdaptInitializePackage - This function initializes everything in the TSGLLEAdapt package. It is
104:   called from TSInitializePackage().

106:   Level: developer

108: .seealso: PetscInitialize()
109: @*/
110: PetscErrorCode  TSGLLEAdaptInitializePackage(void)
111: {

115:   if (TSGLLEAdaptPackageInitialized) return(0);
116:   TSGLLEAdaptPackageInitialized = PETSC_TRUE;
117:   PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID);
118:   TSGLLEAdaptRegisterAll();
119:   PetscRegisterFinalize(TSGLLEAdaptFinalizePackage);
120:   return(0);
121: }

123: PetscErrorCode  TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)
124: {
125:   PetscErrorCode ierr,(*r)(TSGLLEAdapt);

128:   PetscFunctionListFind(TSGLLEAdaptList,type,&r);
129:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAdapt type \"%s\" given",type);
130:   if (((PetscObject)adapt)->type_name) {(*adapt->ops->destroy)(adapt);}
131:   (*r)(adapt);
132:   PetscObjectChangeTypeName((PetscObject)adapt,type);
133:   return(0);
134: }

136: PetscErrorCode  TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])
137: {

141:   PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
142:   return(0);
143: }

145: PetscErrorCode  TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)
146: {
148:   PetscBool      iascii;

151:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
152:   if (iascii) {
153:     PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
154:     if (adapt->ops->view) {
155:       PetscViewerASCIIPushTab(viewer);
156:       (*adapt->ops->view)(adapt,viewer);
157:       PetscViewerASCIIPopTab(viewer);
158:     }
159:   }
160:   return(0);
161: }

163: PetscErrorCode  TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
164: {

168:   if (!*adapt) return(0);
170:   if (--((PetscObject)(*adapt))->refct > 0) {*adapt = 0; return(0);}
171:   if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
172:   PetscHeaderDestroy(adapt);
173:   return(0);
174: }

176: PetscErrorCode  TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt)
177: {
179:   char           type[256] = TSGLLEADAPT_BOTH;
180:   PetscBool      flg;

183:   /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
184:   * function can only be called from inside TSSetFromOptions_GLLE()  */
185:   PetscOptionsHead(PetscOptionsObject,"TSGLLE Adaptivity options");
186:   PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLLEAdaptSetType",TSGLLEAdaptList,
187:                           ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
188:   if (flg || !((PetscObject)adapt)->type_name) {
189:     TSGLLEAdaptSetType(adapt,type);
190:   }
191:   if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
192:   PetscOptionsTail();
193:   return(0);
194: }

196: 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)
197: {

208:   (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);
209:   return(0);
210: }

212: PetscErrorCode  TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt)
213: {
215:   TSGLLEAdapt      adapt;

218:   *inadapt = NULL;
219:   PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView);
220:   *inadapt = adapt;
221:   return(0);
222: }


225: /*
226: *  Implementations
227: */

229: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
230: {

234:   PetscFree(adapt->data);
235:   return(0);
236: }

238: /* -------------------------------- None ----------------------------------- */
239: typedef struct {
240:   PetscInt  scheme;
241:   PetscReal h;
242: } TSGLLEAdapt_None;

244: 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)
245: {

248:   *next_sc = cur;
249:   *next_h  = h;
250:   if (*next_h > tleft) {
251:     *finish = PETSC_TRUE;
252:     *next_h = tleft;
253:   } else *finish = PETSC_FALSE;
254:   return(0);
255: }

257: PetscErrorCode  TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
258: {
260:   TSGLLEAdapt_None *a;

263:   PetscNewLog(adapt,&a);
264:   adapt->data         = (void*)a;
265:   adapt->ops->choose  = TSGLLEAdaptChoose_None;
266:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
267:   return(0);
268: }

270: /* -------------------------------- Size ----------------------------------- */
271: typedef struct {
272:   PetscReal desired_h;
273: } TSGLLEAdapt_Size;


276: 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)
277: {
278:   TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data;
279:   PetscReal      dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;

282:   *next_sc = cur;
283:   optimal  = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur]));
284:   /* Step sizes oscillate when there is no smoothing.  Here we use a geometric mean of the current step size and the
285:   * one that would have been taken (without smoothing) on the last step. */
286:   last_desired_h = sz->desired_h;
287:   sz->desired_h  = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */

289:   /* Normally only happens on the first step */
290:   if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
291:   else *next_h = sz->desired_h;

293:   if (*next_h > tleft) {
294:     *finish = PETSC_TRUE;
295:     *next_h = tleft;
296:   } else *finish = PETSC_FALSE;
297:   return(0);
298: }

300: PetscErrorCode  TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
301: {
303:   TSGLLEAdapt_Size *a;

306:   PetscNewLog(adapt,&a);
307:   adapt->data         = (void*)a;
308:   adapt->ops->choose  = TSGLLEAdaptChoose_Size;
309:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
310:   return(0);
311: }

313: /* -------------------------------- Both ----------------------------------- */
314: typedef struct {
315:   PetscInt  count_at_order;
316:   PetscReal desired_h;
317: } TSGLLEAdapt_Both;


320: 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)
321: {
322:   TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data;
323:   PetscErrorCode   ierr;
324:   PetscReal        dec = 0.2,inc = 5.0,safe = 0.9;
325:   struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
326:   PetscInt        i;

329:   for (i=0; i<n; i++) {
330:     PetscReal optimal;
331:     trial.id  = i;
332:     optimal   = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
333:     trial.h   = h*optimal;
334:     trial.eff = trial.h/cost[i];
335:     if (trial.eff > best.eff) {PetscArraycpy(&best,&trial,1);}
336:     if (i == cur) {PetscArraycpy(&current,&trial,1);}
337:   }
338:   /* Only switch orders if the scheme offers significant benefits over the current one.
339:   When the scheme is not changing, only change step size if it offers significant benefits. */
340:   if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
341:     PetscReal last_desired_h;
342:     *next_sc        = current.id;
343:     last_desired_h  = both->desired_h;
344:     both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
345:     *next_h         = (both->count_at_order > 0)
346:                       ? PetscSqrtReal(last_desired_h * both->desired_h)
347:                       : both->desired_h;
348:     both->count_at_order++;
349:   } else {
350:     PetscReal rat = cost[best.id]/cost[cur];
351:     *next_sc = best.id;
352:     *next_h  = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
353:     both->count_at_order = 0;
354:     both->desired_h      = best.h;
355:   }

357:   if (*next_h > tleft) {
358:     *finish = PETSC_TRUE;
359:     *next_h = tleft;
360:   } else *finish = PETSC_FALSE;
361:   return(0);
362: }

364: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
365: {
367:   TSGLLEAdapt_Both *a;

370:   PetscNewLog(adapt,&a);
371:   adapt->data         = (void*)a;
372:   adapt->ops->choose  = TSGLLEAdaptChoose_Both;
373:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
374:   return(0);
375: }