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 CHKERRQ(TSGLLEAdaptInitializePackage()); 55 CHKERRQ(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 CHKERRQ(TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None)); 74 CHKERRQ(TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size)); 75 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID)); 110 CHKERRQ(TSGLLEAdaptRegisterAll()); 111 CHKERRQ(PetscRegisterFinalize(TSGLLEAdaptFinalizePackage)); 112 PetscFunctionReturn(0); 113 } 114 115 PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type) 116 { 117 PetscErrorCode (*r)(TSGLLEAdapt); 118 119 PetscFunctionBegin; 120 CHKERRQ(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) CHKERRQ((*adapt->ops->destroy)(adapt)); 123 CHKERRQ((*r)(adapt)); 124 CHKERRQ(PetscObjectChangeTypeName((PetscObject)adapt,type)); 125 PetscFunctionReturn(0); 126 } 127 128 PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[]) 129 { 130 PetscFunctionBegin; 131 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix)); 132 PetscFunctionReturn(0); 133 } 134 135 PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer) 136 { 137 PetscBool iascii; 138 139 PetscFunctionBegin; 140 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 141 if (iascii) { 142 CHKERRQ(PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer)); 143 if (adapt->ops->view) { 144 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 145 CHKERRQ((*adapt->ops->view)(adapt,viewer)); 146 CHKERRQ(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) CHKERRQ((*(*adapt)->ops->destroy)(*adapt)); 159 CHKERRQ(PetscHeaderDestroy(adapt)); 160 PetscFunctionReturn(0); 161 } 162 163 PetscErrorCode TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt) 164 { 165 PetscErrorCode ierr; 166 char type[256] = TSGLLEADAPT_BOTH; 167 PetscBool flg; 168 169 PetscFunctionBegin; 170 /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this 171 * function can only be called from inside TSSetFromOptions_GLLE() */ 172 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"TSGLLE Adaptivity options")); 173 ierr = PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLLEAdaptSetType",TSGLLEAdaptList, 174 ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);CHKERRQ(ierr); 175 if (flg || !((PetscObject)adapt)->type_name) { 176 CHKERRQ(TSGLLEAdaptSetType(adapt,type)); 177 } 178 if (adapt->ops->setfromoptions) CHKERRQ((*adapt->ops->setfromoptions)(PetscOptionsObject,adapt)); 179 CHKERRQ(PetscOptionsTail()); 180 PetscFunctionReturn(0); 181 } 182 183 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) 184 { 185 PetscFunctionBegin; 186 PetscValidHeaderSpecific(adapt,TSGLLEADAPT_CLASSID,1); 187 PetscValidIntPointer(orders,3); 188 PetscValidPointer(errors,4); 189 PetscValidPointer(cost,5); 190 PetscValidIntPointer(next_sc,9); 191 PetscValidPointer(next_h,10); 192 PetscValidBoolPointer(finish,11); 193 CHKERRQ((*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish)); 194 PetscFunctionReturn(0); 195 } 196 197 PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt) 198 { 199 TSGLLEAdapt adapt; 200 201 PetscFunctionBegin; 202 *inadapt = NULL; 203 CHKERRQ(PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView)); 204 *inadapt = adapt; 205 PetscFunctionReturn(0); 206 } 207 208 /* 209 Implementations 210 */ 211 212 static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt) 213 { 214 PetscFunctionBegin; 215 CHKERRQ(PetscFree(adapt->data)); 216 PetscFunctionReturn(0); 217 } 218 219 /* -------------------------------- None ----------------------------------- */ 220 typedef struct { 221 PetscInt scheme; 222 PetscReal h; 223 } TSGLLEAdapt_None; 224 225 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) 226 { 227 PetscFunctionBegin; 228 *next_sc = cur; 229 *next_h = h; 230 if (*next_h > tleft) { 231 *finish = PETSC_TRUE; 232 *next_h = tleft; 233 } else *finish = PETSC_FALSE; 234 PetscFunctionReturn(0); 235 } 236 237 PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt) 238 { 239 TSGLLEAdapt_None *a; 240 241 PetscFunctionBegin; 242 CHKERRQ(PetscNewLog(adapt,&a)); 243 adapt->data = (void*)a; 244 adapt->ops->choose = TSGLLEAdaptChoose_None; 245 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 246 PetscFunctionReturn(0); 247 } 248 249 /* -------------------------------- Size ----------------------------------- */ 250 typedef struct { 251 PetscReal desired_h; 252 } TSGLLEAdapt_Size; 253 254 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) 255 { 256 TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data; 257 PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h; 258 259 PetscFunctionBegin; 260 *next_sc = cur; 261 optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur])); 262 /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the 263 * one that would have been taken (without smoothing) on the last step. */ 264 last_desired_h = sz->desired_h; 265 sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */ 266 267 /* Normally only happens on the first step */ 268 if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h); 269 else *next_h = sz->desired_h; 270 271 if (*next_h > tleft) { 272 *finish = PETSC_TRUE; 273 *next_h = tleft; 274 } else *finish = PETSC_FALSE; 275 PetscFunctionReturn(0); 276 } 277 278 PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt) 279 { 280 TSGLLEAdapt_Size *a; 281 282 PetscFunctionBegin; 283 CHKERRQ(PetscNewLog(adapt,&a)); 284 adapt->data = (void*)a; 285 adapt->ops->choose = TSGLLEAdaptChoose_Size; 286 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 287 PetscFunctionReturn(0); 288 } 289 290 /* -------------------------------- Both ----------------------------------- */ 291 typedef struct { 292 PetscInt count_at_order; 293 PetscReal desired_h; 294 } TSGLLEAdapt_Both; 295 296 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) 297 { 298 TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data; 299 PetscReal dec = 0.2,inc = 5.0,safe = 0.9; 300 struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0}; 301 PetscInt i; 302 303 PetscFunctionBegin; 304 for (i=0; i<n; i++) { 305 PetscReal optimal; 306 trial.id = i; 307 optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i])); 308 trial.h = h*optimal; 309 trial.eff = trial.h/cost[i]; 310 if (trial.eff > best.eff) CHKERRQ(PetscArraycpy(&best,&trial,1)); 311 if (i == cur) CHKERRQ(PetscArraycpy(¤t,&trial,1)); 312 } 313 /* Only switch orders if the scheme offers significant benefits over the current one. 314 When the scheme is not changing, only change step size if it offers significant benefits. */ 315 if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) { 316 PetscReal last_desired_h; 317 *next_sc = current.id; 318 last_desired_h = both->desired_h; 319 both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h)); 320 *next_h = (both->count_at_order > 0) 321 ? PetscSqrtReal(last_desired_h * both->desired_h) 322 : both->desired_h; 323 both->count_at_order++; 324 } else { 325 PetscReal rat = cost[best.id]/cost[cur]; 326 *next_sc = best.id; 327 *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h)); 328 both->count_at_order = 0; 329 both->desired_h = best.h; 330 } 331 332 if (*next_h > tleft) { 333 *finish = PETSC_TRUE; 334 *next_h = tleft; 335 } else *finish = PETSC_FALSE; 336 PetscFunctionReturn(0); 337 } 338 339 PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt) 340 { 341 TSGLLEAdapt_Both *a; 342 343 PetscFunctionBegin; 344 CHKERRQ(PetscNewLog(adapt,&a)); 345 adapt->data = (void*)a; 346 adapt->ops->choose = TSGLLEAdaptChoose_Both; 347 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 348 PetscFunctionReturn(0); 349 } 350