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