1 2 #include <../src/ts/impls/implicit/glle/glle.h> /*I "petscts.h" I*/ 3 4 static PetscFunctionList TSGLLEAdaptList; 5 static PetscBool TSGLLEAdaptPackageInitialized; 6 static PetscBool TSGLLEAdaptRegisterAllCalled; 7 static PetscClassId TSGLLEADAPT_CLASSID; 8 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 }; 15 16 struct _p_TSGLLEAdapt { 17 PETSCHEADER(struct _TSGLLEAdaptOps); 18 void *data; 19 }; 20 21 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt); 22 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt); 23 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt); 24 25 /*@C 26 TSGLLEAdaptRegister - adds a TSGLLEAdapt implementation 27 28 Not Collective 29 30 Input Parameters: 31 + name_scheme - name of user-defined adaptivity scheme 32 - routine_create - routine to create method context 33 34 Notes: 35 TSGLLEAdaptRegister() may be called multiple times to add several user-defined families. 36 37 Sample usage: 38 .vb 39 TSGLLEAdaptRegister("my_scheme",MySchemeCreate); 40 .ve 41 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 46 47 Level: advanced 48 49 .seealso: TSGLLEAdaptRegisterAll() 50 @*/ 51 PetscErrorCode TSGLLEAdaptRegister(const char sname[],PetscErrorCode (*function)(TSGLLEAdapt)) 52 { 53 PetscFunctionBegin; 54 PetscCall(TSGLLEAdaptInitializePackage()); 55 PetscCall(PetscFunctionListAdd(&TSGLLEAdaptList,sname,function)); 56 PetscFunctionReturn(0); 57 } 58 59 /*@C 60 TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in TSGLLEAdapt 61 62 Not Collective 63 64 Level: advanced 65 66 .seealso: TSGLLEAdaptRegisterDestroy() 67 @*/ 68 PetscErrorCode TSGLLEAdaptRegisterAll(void) 69 { 70 PetscFunctionBegin; 71 if (TSGLLEAdaptRegisterAllCalled) PetscFunctionReturn(0); 72 TSGLLEAdaptRegisterAllCalled = PETSC_TRUE; 73 PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None)); 74 PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size)); 75 PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_BOTH,TSGLLEAdaptCreate_Both)); 76 PetscFunctionReturn(0); 77 } 78 79 /*@C 80 TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is 81 called from PetscFinalize(). 82 83 Level: developer 84 85 .seealso: PetscFinalize() 86 @*/ 87 PetscErrorCode TSGLLEAdaptFinalizePackage(void) 88 { 89 PetscFunctionBegin; 90 PetscCall(PetscFunctionListDestroy(&TSGLLEAdaptList)); 91 TSGLLEAdaptPackageInitialized = PETSC_FALSE; 92 TSGLLEAdaptRegisterAllCalled = PETSC_FALSE; 93 PetscFunctionReturn(0); 94 } 95 96 /*@C 97 TSGLLEAdaptInitializePackage - This function initializes everything in the TSGLLEAdapt package. It is 98 called from TSInitializePackage(). 99 100 Level: developer 101 102 .seealso: PetscInitialize() 103 @*/ 104 PetscErrorCode TSGLLEAdaptInitializePackage(void) 105 { 106 PetscFunctionBegin; 107 if (TSGLLEAdaptPackageInitialized) PetscFunctionReturn(0); 108 TSGLLEAdaptPackageInitialized = PETSC_TRUE; 109 PetscCall(PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID)); 110 PetscCall(TSGLLEAdaptRegisterAll()); 111 PetscCall(PetscRegisterFinalize(TSGLLEAdaptFinalizePackage)); 112 PetscFunctionReturn(0); 113 } 114 115 PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type) 116 { 117 PetscErrorCode (*r)(TSGLLEAdapt); 118 119 PetscFunctionBegin; 120 PetscCall(PetscFunctionListFind(TSGLLEAdaptList,type,&r)); 121 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAdapt type \"%s\" given",type); 122 if (((PetscObject)adapt)->type_name) PetscCall((*adapt->ops->destroy)(adapt)); 123 PetscCall((*r)(adapt)); 124 PetscCall(PetscObjectChangeTypeName((PetscObject)adapt,type)); 125 PetscFunctionReturn(0); 126 } 127 128 PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[]) 129 { 130 PetscFunctionBegin; 131 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix)); 132 PetscFunctionReturn(0); 133 } 134 135 PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer) 136 { 137 PetscBool iascii; 138 139 PetscFunctionBegin; 140 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 141 if (iascii) { 142 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer)); 143 if (adapt->ops->view) { 144 PetscCall(PetscViewerASCIIPushTab(viewer)); 145 PetscCall((*adapt->ops->view)(adapt,viewer)); 146 PetscCall(PetscViewerASCIIPopTab(viewer)); 147 } 148 } 149 PetscFunctionReturn(0); 150 } 151 152 PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt) 153 { 154 PetscFunctionBegin; 155 if (!*adapt) PetscFunctionReturn(0); 156 PetscValidHeaderSpecific(*adapt,TSGLLEADAPT_CLASSID,1); 157 if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; PetscFunctionReturn(0);} 158 if ((*adapt)->ops->destroy) PetscCall((*(*adapt)->ops->destroy)(*adapt)); 159 PetscCall(PetscHeaderDestroy(adapt)); 160 PetscFunctionReturn(0); 161 } 162 163 PetscErrorCode TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt) 164 { 165 char type[256] = TSGLLEADAPT_BOTH; 166 PetscBool flg; 167 168 PetscFunctionBegin; 169 /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this 170 * function can only be called from inside TSSetFromOptions_GLLE() */ 171 PetscCall(PetscOptionsHead(PetscOptionsObject,"TSGLLE Adaptivity options")); 172 PetscCall(PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLLEAdaptSetType",TSGLLEAdaptList, 173 ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg)); 174 if (flg || !((PetscObject)adapt)->type_name) PetscCall(TSGLLEAdaptSetType(adapt,type)); 175 if (adapt->ops->setfromoptions) PetscCall((*adapt->ops->setfromoptions)(PetscOptionsObject,adapt)); 176 PetscCall(PetscOptionsTail()); 177 PetscFunctionReturn(0); 178 } 179 180 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) 181 { 182 PetscFunctionBegin; 183 PetscValidHeaderSpecific(adapt,TSGLLEADAPT_CLASSID,1); 184 PetscValidIntPointer(orders,3); 185 PetscValidRealPointer(errors,4); 186 PetscValidRealPointer(cost,5); 187 PetscValidIntPointer(next_sc,9); 188 PetscValidRealPointer(next_h,10); 189 PetscValidBoolPointer(finish,11); 190 PetscCall((*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish)); 191 PetscFunctionReturn(0); 192 } 193 194 PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt) 195 { 196 TSGLLEAdapt adapt; 197 198 PetscFunctionBegin; 199 *inadapt = NULL; 200 PetscCall(PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView)); 201 *inadapt = adapt; 202 PetscFunctionReturn(0); 203 } 204 205 /* 206 Implementations 207 */ 208 209 static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt) 210 { 211 PetscFunctionBegin; 212 PetscCall(PetscFree(adapt->data)); 213 PetscFunctionReturn(0); 214 } 215 216 /* -------------------------------- None ----------------------------------- */ 217 typedef struct { 218 PetscInt scheme; 219 PetscReal h; 220 } TSGLLEAdapt_None; 221 222 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) 223 { 224 PetscFunctionBegin; 225 *next_sc = cur; 226 *next_h = h; 227 if (*next_h > tleft) { 228 *finish = PETSC_TRUE; 229 *next_h = tleft; 230 } else *finish = PETSC_FALSE; 231 PetscFunctionReturn(0); 232 } 233 234 PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt) 235 { 236 TSGLLEAdapt_None *a; 237 238 PetscFunctionBegin; 239 PetscCall(PetscNewLog(adapt,&a)); 240 adapt->data = (void*)a; 241 adapt->ops->choose = TSGLLEAdaptChoose_None; 242 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 243 PetscFunctionReturn(0); 244 } 245 246 /* -------------------------------- Size ----------------------------------- */ 247 typedef struct { 248 PetscReal desired_h; 249 } TSGLLEAdapt_Size; 250 251 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) 252 { 253 TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data; 254 PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h; 255 256 PetscFunctionBegin; 257 *next_sc = cur; 258 optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur])); 259 /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the 260 * one that would have been taken (without smoothing) on the last step. */ 261 last_desired_h = sz->desired_h; 262 sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */ 263 264 /* Normally only happens on the first step */ 265 if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h); 266 else *next_h = sz->desired_h; 267 268 if (*next_h > tleft) { 269 *finish = PETSC_TRUE; 270 *next_h = tleft; 271 } else *finish = PETSC_FALSE; 272 PetscFunctionReturn(0); 273 } 274 275 PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt) 276 { 277 TSGLLEAdapt_Size *a; 278 279 PetscFunctionBegin; 280 PetscCall(PetscNewLog(adapt,&a)); 281 adapt->data = (void*)a; 282 adapt->ops->choose = TSGLLEAdaptChoose_Size; 283 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 284 PetscFunctionReturn(0); 285 } 286 287 /* -------------------------------- Both ----------------------------------- */ 288 typedef struct { 289 PetscInt count_at_order; 290 PetscReal desired_h; 291 } TSGLLEAdapt_Both; 292 293 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) 294 { 295 TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data; 296 PetscReal dec = 0.2,inc = 5.0,safe = 0.9; 297 struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0}; 298 PetscInt i; 299 300 PetscFunctionBegin; 301 for (i=0; i<n; i++) { 302 PetscReal optimal; 303 trial.id = i; 304 optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i])); 305 trial.h = h*optimal; 306 trial.eff = trial.h/cost[i]; 307 if (trial.eff > best.eff) PetscCall(PetscArraycpy(&best,&trial,1)); 308 if (i == cur) PetscCall(PetscArraycpy(¤t,&trial,1)); 309 } 310 /* Only switch orders if the scheme offers significant benefits over the current one. 311 When the scheme is not changing, only change step size if it offers significant benefits. */ 312 if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) { 313 PetscReal last_desired_h; 314 *next_sc = current.id; 315 last_desired_h = both->desired_h; 316 both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h)); 317 *next_h = (both->count_at_order > 0) 318 ? PetscSqrtReal(last_desired_h * both->desired_h) 319 : both->desired_h; 320 both->count_at_order++; 321 } else { 322 PetscReal rat = cost[best.id]/cost[cur]; 323 *next_sc = best.id; 324 *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h)); 325 both->count_at_order = 0; 326 both->desired_h = best.h; 327 } 328 329 if (*next_h > tleft) { 330 *finish = PETSC_TRUE; 331 *next_h = tleft; 332 } else *finish = PETSC_FALSE; 333 PetscFunctionReturn(0); 334 } 335 336 PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt) 337 { 338 TSGLLEAdapt_Both *a; 339 340 PetscFunctionBegin; 341 PetscCall(PetscNewLog(adapt,&a)); 342 adapt->data = (void*)a; 343 adapt->ops->choose = TSGLLEAdaptChoose_Both; 344 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 345 PetscFunctionReturn(0); 346 } 347