Actual source code: glleadapt.c


  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: {
 53:   TSGLLEAdaptInitializePackage();
 54:   PetscFunctionListAdd(&TSGLLEAdaptList,sname,function);
 55:   return 0;
 56: }

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

 61:   Not Collective

 63:   Level: advanced

 65: .seealso: TSGLLEAdaptRegisterDestroy()
 66: @*/
 67: PetscErrorCode  TSGLLEAdaptRegisterAll(void)
 68: {
 69:   if (TSGLLEAdaptRegisterAllCalled) return 0;
 70:   TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
 71:   TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None);
 72:   TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size);
 73:   TSGLLEAdaptRegister(TSGLLEADAPT_BOTH,TSGLLEAdaptCreate_Both);
 74:   return 0;
 75: }

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

 81:   Level: developer

 83: .seealso: PetscFinalize()
 84: @*/
 85: PetscErrorCode  TSGLLEAdaptFinalizePackage(void)
 86: {
 87:   PetscFunctionListDestroy(&TSGLLEAdaptList);
 88:   TSGLLEAdaptPackageInitialized = PETSC_FALSE;
 89:   TSGLLEAdaptRegisterAllCalled  = PETSC_FALSE;
 90:   return 0;
 91: }

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

 97:   Level: developer

 99: .seealso: PetscInitialize()
100: @*/
101: PetscErrorCode  TSGLLEAdaptInitializePackage(void)
102: {
103:   if (TSGLLEAdaptPackageInitialized) return 0;
104:   TSGLLEAdaptPackageInitialized = PETSC_TRUE;
105:   PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID);
106:   TSGLLEAdaptRegisterAll();
107:   PetscRegisterFinalize(TSGLLEAdaptFinalizePackage);
108:   return 0;
109: }

111: PetscErrorCode  TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)
112: {
113:   PetscErrorCode (*r)(TSGLLEAdapt);

115:   PetscFunctionListFind(TSGLLEAdaptList,type,&r);
117:   if (((PetscObject)adapt)->type_name) (*adapt->ops->destroy)(adapt);
118:   (*r)(adapt);
119:   PetscObjectChangeTypeName((PetscObject)adapt,type);
120:   return 0;
121: }

123: PetscErrorCode  TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])
124: {
125:   PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
126:   return 0;
127: }

129: PetscErrorCode  TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)
130: {
131:   PetscBool      iascii;

133:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
134:   if (iascii) {
135:     PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
136:     if (adapt->ops->view) {
137:       PetscViewerASCIIPushTab(viewer);
138:       (*adapt->ops->view)(adapt,viewer);
139:       PetscViewerASCIIPopTab(viewer);
140:     }
141:   }
142:   return 0;
143: }

145: PetscErrorCode  TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
146: {
147:   if (!*adapt) return 0;
149:   if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; return 0;}
150:   if ((*adapt)->ops->destroy) (*(*adapt)->ops->destroy)(*adapt);
151:   PetscHeaderDestroy(adapt);
152:   return 0;
153: }

155: PetscErrorCode  TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt)
156: {
157:   char      type[256] = TSGLLEADAPT_BOTH;
158:   PetscBool flg;

160:   /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
161:   * function can only be called from inside TSSetFromOptions_GLLE()  */
162:   PetscOptionsHead(PetscOptionsObject,"TSGLLE Adaptivity options");
163:   PetscCall(PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLLEAdaptSetType",TSGLLEAdaptList,
164:                             ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg));
165:   if (flg || !((PetscObject)adapt)->type_name) TSGLLEAdaptSetType(adapt,type);
166:   if (adapt->ops->setfromoptions) (*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);
167:   PetscOptionsTail();
168:   return 0;
169: }

171: 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)
172: {
180:   (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);
181:   return 0;
182: }

184: PetscErrorCode  TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt)
185: {
186:   TSGLLEAdapt      adapt;

188:   *inadapt = NULL;
189:   PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView);
190:   *inadapt = adapt;
191:   return 0;
192: }

194: /*
195:    Implementations
196: */

198: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
199: {
200:   PetscFree(adapt->data);
201:   return 0;
202: }

204: /* -------------------------------- None ----------------------------------- */
205: typedef struct {
206:   PetscInt  scheme;
207:   PetscReal h;
208: } TSGLLEAdapt_None;

210: 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)
211: {
212:   *next_sc = cur;
213:   *next_h  = h;
214:   if (*next_h > tleft) {
215:     *finish = PETSC_TRUE;
216:     *next_h = tleft;
217:   } else *finish = PETSC_FALSE;
218:   return 0;
219: }

221: PetscErrorCode  TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
222: {
223:   TSGLLEAdapt_None *a;

225:   PetscNewLog(adapt,&a);
226:   adapt->data         = (void*)a;
227:   adapt->ops->choose  = TSGLLEAdaptChoose_None;
228:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
229:   return 0;
230: }

232: /* -------------------------------- Size ----------------------------------- */
233: typedef struct {
234:   PetscReal desired_h;
235: } TSGLLEAdapt_Size;

237: 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)
238: {
239:   TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data;
240:   PetscReal      dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;

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

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

253:   if (*next_h > tleft) {
254:     *finish = PETSC_TRUE;
255:     *next_h = tleft;
256:   } else *finish = PETSC_FALSE;
257:   return 0;
258: }

260: PetscErrorCode  TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
261: {
262:   TSGLLEAdapt_Size *a;

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

271: /* -------------------------------- Both ----------------------------------- */
272: typedef struct {
273:   PetscInt  count_at_order;
274:   PetscReal desired_h;
275: } TSGLLEAdapt_Both;

277: 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)
278: {
279:   TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data;
280:   PetscReal        dec = 0.2,inc = 5.0,safe = 0.9;
281:   struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
282:   PetscInt        i;

284:   for (i=0; i<n; i++) {
285:     PetscReal optimal;
286:     trial.id  = i;
287:     optimal   = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
288:     trial.h   = h*optimal;
289:     trial.eff = trial.h/cost[i];
290:     if (trial.eff > best.eff) PetscArraycpy(&best,&trial,1);
291:     if (i == cur) PetscArraycpy(&current,&trial,1);
292:   }
293:   /* Only switch orders if the scheme offers significant benefits over the current one.
294:   When the scheme is not changing, only change step size if it offers significant benefits. */
295:   if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
296:     PetscReal last_desired_h;
297:     *next_sc        = current.id;
298:     last_desired_h  = both->desired_h;
299:     both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
300:     *next_h         = (both->count_at_order > 0)
301:                       ? PetscSqrtReal(last_desired_h * both->desired_h)
302:                       : both->desired_h;
303:     both->count_at_order++;
304:   } else {
305:     PetscReal rat = cost[best.id]/cost[cur];
306:     *next_sc = best.id;
307:     *next_h  = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
308:     both->count_at_order = 0;
309:     both->desired_h      = best.h;
310:   }

312:   if (*next_h > tleft) {
313:     *finish = PETSC_TRUE;
314:     *next_h = tleft;
315:   } else *finish = PETSC_FALSE;
316:   return 0;
317: }

319: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
320: {
321:   TSGLLEAdapt_Both *a;

323:   PetscNewLog(adapt,&a);
324:   adapt->data         = (void*)a;
325:   adapt->ops->choose  = TSGLLEAdaptChoose_Both;
326:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
327:   return 0;
328: }