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