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