Actual source code: glleadapt.c

petsc-3.11.4 2019-09-28
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: .keywords: TSGLLEAdapt, register

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

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

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

 66:   Not Collective

 68:   Level: advanced

 70: .keywords: TSGLLEAdapt, register, all

 72: .seealso: TSGLLEAdaptRegisterDestroy()
 73: @*/
 74: PetscErrorCode  TSGLLEAdaptRegisterAll(void)
 75: {

 79:   if (TSGLLEAdaptRegisterAllCalled) return(0);
 80:   TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
 81:   TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None);
 82:   TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size);
 83:   TSGLLEAdaptRegister(TSGLLEADAPT_BOTH,TSGLLEAdaptCreate_Both);
 84:   return(0);
 85: }

 87: /*@C
 88:   TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is
 89:   called from PetscFinalize().

 91:   Level: developer

 93: .keywords: Petsc, destroy, package
 94: .seealso: PetscFinalize()
 95: @*/
 96: PetscErrorCode  TSGLLEAdaptFinalizePackage(void)
 97: {

101:   PetscFunctionListDestroy(&TSGLLEAdaptList);
102:   TSGLLEAdaptPackageInitialized = PETSC_FALSE;
103:   TSGLLEAdaptRegisterAllCalled  = PETSC_FALSE;
104:   return(0);
105: }

107: /*@C
108:   TSGLLEAdaptInitializePackage - This function initializes everything in the TSGLLEAdapt package. It is
109:   called from TSInitializePackage().

111:   Level: developer

113: .keywords: TSGLLEAdapt, initialize, package
114: .seealso: PetscInitialize()
115: @*/
116: PetscErrorCode  TSGLLEAdaptInitializePackage(void)
117: {

121:   if (TSGLLEAdaptPackageInitialized) return(0);
122:   TSGLLEAdaptPackageInitialized = PETSC_TRUE;
123:   PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID);
124:   TSGLLEAdaptRegisterAll();
125:   PetscRegisterFinalize(TSGLLEAdaptFinalizePackage);
126:   return(0);
127: }

129: PetscErrorCode  TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)
130: {
131:   PetscErrorCode ierr,(*r)(TSGLLEAdapt);

134:   PetscFunctionListFind(TSGLLEAdaptList,type,&r);
135:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAdapt type \"%s\" given",type);
136:   if (((PetscObject)adapt)->type_name) {(*adapt->ops->destroy)(adapt);}
137:   (*r)(adapt);
138:   PetscObjectChangeTypeName((PetscObject)adapt,type);
139:   return(0);
140: }

142: PetscErrorCode  TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])
143: {

147:   PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
148:   return(0);
149: }

151: PetscErrorCode  TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)
152: {
154:   PetscBool      iascii;

157:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
158:   if (iascii) {
159:     PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
160:     if (adapt->ops->view) {
161:       PetscViewerASCIIPushTab(viewer);
162:       (*adapt->ops->view)(adapt,viewer);
163:       PetscViewerASCIIPopTab(viewer);
164:     }
165:   }
166:   return(0);
167: }

169: PetscErrorCode  TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
170: {

174:   if (!*adapt) return(0);
176:   if (--((PetscObject)(*adapt))->refct > 0) {*adapt = 0; return(0);}
177:   if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
178:   PetscHeaderDestroy(adapt);
179:   return(0);
180: }

182: PetscErrorCode  TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt)
183: {
185:   char           type[256] = TSGLLEADAPT_BOTH;
186:   PetscBool      flg;

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

202: 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)
203: {

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

218: PetscErrorCode  TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt)
219: {
221:   TSGLLEAdapt      adapt;

224:   *inadapt = NULL;
225:   PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView);
226:   *inadapt = adapt;
227:   return(0);
228: }


231: /*
232: *  Implementations
233: */

235: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
236: {

240:   PetscFree(adapt->data);
241:   return(0);
242: }

244: /* -------------------------------- None ----------------------------------- */
245: typedef struct {
246:   PetscInt  scheme;
247:   PetscReal h;
248: } TSGLLEAdapt_None;

250: 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)
251: {

254:   *next_sc = cur;
255:   *next_h  = h;
256:   if (*next_h > tleft) {
257:     *finish = PETSC_TRUE;
258:     *next_h = tleft;
259:   } else *finish = PETSC_FALSE;
260:   return(0);
261: }

263: PetscErrorCode  TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
264: {
266:   TSGLLEAdapt_None *a;

269:   PetscNewLog(adapt,&a);
270:   adapt->data         = (void*)a;
271:   adapt->ops->choose  = TSGLLEAdaptChoose_None;
272:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
273:   return(0);
274: }

276: /* -------------------------------- Size ----------------------------------- */
277: typedef struct {
278:   PetscReal desired_h;
279: } TSGLLEAdapt_Size;


282: 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)
283: {
284:   TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data;
285:   PetscReal      dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;

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

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

299:   if (*next_h > tleft) {
300:     *finish = PETSC_TRUE;
301:     *next_h = tleft;
302:   } else *finish = PETSC_FALSE;
303:   return(0);
304: }

306: PetscErrorCode  TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
307: {
309:   TSGLLEAdapt_Size *a;

312:   PetscNewLog(adapt,&a);
313:   adapt->data         = (void*)a;
314:   adapt->ops->choose  = TSGLLEAdaptChoose_Size;
315:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
316:   return(0);
317: }

319: /* -------------------------------- Both ----------------------------------- */
320: typedef struct {
321:   PetscInt  count_at_order;
322:   PetscReal desired_h;
323: } TSGLLEAdapt_Both;


326: 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)
327: {
328:   TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data;
330:   PetscReal      dec = 0.2,inc = 5.0,safe = 0.9;
331:   struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
332:   PetscInt i;

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

363:   if (*next_h > tleft) {
364:     *finish = PETSC_TRUE;
365:     *next_h = tleft;
366:   } else *finish = PETSC_FALSE;
367:   return(0);
368: }

370: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
371: {
373:   TSGLLEAdapt_Both *a;

376:   PetscNewLog(adapt,&a);
377:   adapt->data         = (void*)a;
378:   adapt->ops->choose  = TSGLLEAdaptChoose_Both;
379:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
380:   return(0);
381: }