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 PetscValidBoolPointer(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 Implementations 226 */ 227 228 static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt) 229 { 230 PetscErrorCode ierr; 231 232 PetscFunctionBegin; 233 ierr = PetscFree(adapt->data);CHKERRQ(ierr); 234 PetscFunctionReturn(0); 235 } 236 237 /* -------------------------------- None ----------------------------------- */ 238 typedef struct { 239 PetscInt scheme; 240 PetscReal h; 241 } TSGLLEAdapt_None; 242 243 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) 244 { 245 246 PetscFunctionBegin; 247 *next_sc = cur; 248 *next_h = h; 249 if (*next_h > tleft) { 250 *finish = PETSC_TRUE; 251 *next_h = tleft; 252 } else *finish = PETSC_FALSE; 253 PetscFunctionReturn(0); 254 } 255 256 PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt) 257 { 258 PetscErrorCode ierr; 259 TSGLLEAdapt_None *a; 260 261 PetscFunctionBegin; 262 ierr = PetscNewLog(adapt,&a);CHKERRQ(ierr); 263 adapt->data = (void*)a; 264 adapt->ops->choose = TSGLLEAdaptChoose_None; 265 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 266 PetscFunctionReturn(0); 267 } 268 269 /* -------------------------------- Size ----------------------------------- */ 270 typedef struct { 271 PetscReal desired_h; 272 } TSGLLEAdapt_Size; 273 274 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) 275 { 276 TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data; 277 PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h; 278 279 PetscFunctionBegin; 280 *next_sc = cur; 281 optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur])); 282 /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the 283 * one that would have been taken (without smoothing) on the last step. */ 284 last_desired_h = sz->desired_h; 285 sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */ 286 287 /* Normally only happens on the first step */ 288 if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h); 289 else *next_h = sz->desired_h; 290 291 if (*next_h > tleft) { 292 *finish = PETSC_TRUE; 293 *next_h = tleft; 294 } else *finish = PETSC_FALSE; 295 PetscFunctionReturn(0); 296 } 297 298 PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt) 299 { 300 PetscErrorCode ierr; 301 TSGLLEAdapt_Size *a; 302 303 PetscFunctionBegin; 304 ierr = PetscNewLog(adapt,&a);CHKERRQ(ierr); 305 adapt->data = (void*)a; 306 adapt->ops->choose = TSGLLEAdaptChoose_Size; 307 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 308 PetscFunctionReturn(0); 309 } 310 311 /* -------------------------------- Both ----------------------------------- */ 312 typedef struct { 313 PetscInt count_at_order; 314 PetscReal desired_h; 315 } TSGLLEAdapt_Both; 316 317 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) 318 { 319 TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data; 320 PetscErrorCode ierr; 321 PetscReal dec = 0.2,inc = 5.0,safe = 0.9; 322 struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0}; 323 PetscInt i; 324 325 PetscFunctionBegin; 326 for (i=0; i<n; i++) { 327 PetscReal optimal; 328 trial.id = i; 329 optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i])); 330 trial.h = h*optimal; 331 trial.eff = trial.h/cost[i]; 332 if (trial.eff > best.eff) {ierr = PetscArraycpy(&best,&trial,1);CHKERRQ(ierr);} 333 if (i == cur) {ierr = PetscArraycpy(¤t,&trial,1);CHKERRQ(ierr);} 334 } 335 /* Only switch orders if the scheme offers significant benefits over the current one. 336 When the scheme is not changing, only change step size if it offers significant benefits. */ 337 if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) { 338 PetscReal last_desired_h; 339 *next_sc = current.id; 340 last_desired_h = both->desired_h; 341 both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h)); 342 *next_h = (both->count_at_order > 0) 343 ? PetscSqrtReal(last_desired_h * both->desired_h) 344 : both->desired_h; 345 both->count_at_order++; 346 } else { 347 PetscReal rat = cost[best.id]/cost[cur]; 348 *next_sc = best.id; 349 *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h)); 350 both->count_at_order = 0; 351 both->desired_h = best.h; 352 } 353 354 if (*next_h > tleft) { 355 *finish = PETSC_TRUE; 356 *next_h = tleft; 357 } else *finish = PETSC_FALSE; 358 PetscFunctionReturn(0); 359 } 360 361 PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt) 362 { 363 PetscErrorCode ierr; 364 TSGLLEAdapt_Both *a; 365 366 PetscFunctionBegin; 367 ierr = PetscNewLog(adapt,&a);CHKERRQ(ierr); 368 adapt->data = (void*)a; 369 adapt->ops->choose = TSGLLEAdaptChoose_Both; 370 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree; 371 PetscFunctionReturn(0); 372 } 373