1 #include <petscconvest.h> /*I "petscconvest.h" I*/ 2 #include <petscts.h> 3 4 #include <petsc/private/petscconvestimpl.h> 5 6 static PetscErrorCode PetscConvEstSetTS_Private(PetscConvEst ce, PetscObject solver) 7 { 8 PetscClassId id; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 ierr = PetscObjectGetClassId(ce->solver, &id);CHKERRQ(ierr); 13 if (id != TS_CLASSID) SETERRQ(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "Solver was not a TS"); 14 ierr = TSGetDM((TS) ce->solver, &ce->idm);CHKERRQ(ierr); 15 PetscFunctionReturn(0); 16 } 17 18 static PetscErrorCode PetscConvEstInitGuessTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u) 19 { 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 ierr = TSComputeInitialCondition((TS) ce->solver, u);CHKERRQ(ierr); 24 PetscFunctionReturn(0); 25 } 26 27 static PetscErrorCode PetscConvEstComputeErrorTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[]) 28 { 29 Vec e; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 ierr = VecDuplicate(u, &e);CHKERRQ(ierr); 34 ierr = TSComputeExactError((TS) ce->solver, u, e);CHKERRQ(ierr); 35 ierr = VecNorm(e, NORM_2, errors);CHKERRQ(ierr); 36 ierr = VecDestroy(&e);CHKERRQ(ierr); 37 PetscFunctionReturn(0); 38 } 39 40 static PetscErrorCode PetscConvEstGetConvRateTS_Private(PetscConvEst ce, PetscReal alpha[]) 41 { 42 TS ts = (TS) ce->solver; 43 Vec u; 44 PetscReal *dt, *x, *y, slope, intercept; 45 PetscInt Ns, oNs, Nf = ce->Nf, f, Nr = ce->Nr, r; 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 ierr = TSGetSolution(ts, &u);CHKERRQ(ierr); 50 ierr = PetscMalloc1(Nr+1, &dt);CHKERRQ(ierr); 51 ierr = TSGetTimeStep(ts, &dt[0]);CHKERRQ(ierr); 52 ierr = TSGetMaxSteps(ts, &oNs);CHKERRQ(ierr); 53 Ns = oNs; 54 for (r = 0; r <= Nr; ++r) { 55 if (r > 0) { 56 dt[r] = dt[r-1]/ce->r; 57 Ns = PetscCeilReal(Ns*ce->r); 58 } 59 ierr = TSSetTime(ts, 0.0);CHKERRQ(ierr); 60 ierr = TSSetStepNumber(ts, 0);CHKERRQ(ierr); 61 ierr = TSSetTimeStep(ts, dt[r]);CHKERRQ(ierr); 62 ierr = TSSetMaxSteps(ts, Ns);CHKERRQ(ierr); 63 ierr = PetscConvEstComputeInitialGuess(ce, r, NULL, u);CHKERRQ(ierr); 64 ierr = TSSolve(ts, u);CHKERRQ(ierr); 65 ierr = PetscLogEventBegin(ce->event, ce, 0, 0, 0);CHKERRQ(ierr); 66 ierr = PetscConvEstComputeError(ce, r, NULL, u, &ce->errors[r*Nf]);CHKERRQ(ierr); 67 ierr = PetscLogEventEnd(ce->event, ce, 0, 0, 0);CHKERRQ(ierr); 68 for (f = 0; f < ce->Nf; ++f) { 69 ierr = PetscLogEventSetDof(ce->event, f, 1.0/dt[r]);CHKERRQ(ierr); 70 ierr = PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);CHKERRQ(ierr); 71 } 72 } 73 /* Fit convergence rate */ 74 if (Nr) { 75 ierr = PetscMalloc2(Nr+1, &x, Nr+1, &y);CHKERRQ(ierr); 76 for (f = 0; f < Nf; ++f) { 77 for (r = 0; r <= Nr; ++r) { 78 x[r] = PetscLog10Real(dt[r*Nf+f]); 79 y[r] = PetscLog10Real(ce->errors[r*Nf+f]); 80 } 81 ierr = PetscLinearRegression(Nr+1, x, y, &slope, &intercept);CHKERRQ(ierr); 82 /* Since lg err = s lg dt + b */ 83 alpha[f] = slope; 84 } 85 ierr = PetscFree2(x, y);CHKERRQ(ierr); 86 } 87 /* Reset solver */ 88 ierr = TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);CHKERRQ(ierr); 89 ierr = TSSetTime(ts, 0.0);CHKERRQ(ierr); 90 ierr = TSSetStepNumber(ts, 0);CHKERRQ(ierr); 91 ierr = TSSetTimeStep(ts, dt[0]);CHKERRQ(ierr); 92 ierr = TSSetMaxSteps(ts, oNs);CHKERRQ(ierr); 93 ierr = PetscConvEstComputeInitialGuess(ce, 0, NULL, u);CHKERRQ(ierr); 94 ierr = PetscFree(dt);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 98 PetscErrorCode PetscConvEstUseTS(PetscConvEst ce) 99 { 100 PetscFunctionBegin; 101 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1); 102 ce->ops->setsolver = PetscConvEstSetTS_Private; 103 ce->ops->initguess = PetscConvEstInitGuessTS_Private; 104 ce->ops->computeerror = PetscConvEstComputeErrorTS_Private; 105 ce->ops->getconvrate = PetscConvEstGetConvRateTS_Private; 106 PetscFunctionReturn(0); 107 } 108