1 #define PETSCTS_DLL 2 3 /* 4 Code for timestepping with implicit Theta method 5 6 Notes: 7 This method can be applied to DAE. 8 9 This method is cast as a 1-stage implicit Runge-Kutta method. 10 11 Theta | Theta 12 ------------- 13 | 1 14 15 To apply a diagonally implicit RK method to DAE, the stage formula 16 17 X_i = x + h sum_j a_ij X'_j 18 19 is interpreted as a formula for X'_i in terms of X_i and known stuff (X'_j, j<i) 20 */ 21 #include "private/tsimpl.h" /*I "petscts.h" I*/ 22 23 typedef struct { 24 Vec Xold; 25 Vec X,Xdot; /* Storage for one stage */ 26 Vec res; /* DAE residuals */ 27 PetscReal Theta; 28 PetscReal shift; 29 PetscReal stage_time; 30 } TS_Theta; 31 32 #undef __FUNCT__ 33 #define __FUNCT__ "TSStep_Theta" 34 static PetscErrorCode TSStep_Theta(TS ts,PetscInt *steps,PetscReal *ptime) 35 { 36 Vec sol = ts->vec_sol; 37 PetscErrorCode ierr; 38 PetscInt i,max_steps = ts->max_steps,its,lits; 39 TS_Theta *th = (TS_Theta*)ts->data; 40 41 PetscFunctionBegin; 42 *steps = -ts->steps; 43 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 44 45 for (i=0; i<max_steps; i++) { 46 if (ts->ptime + ts->time_step > ts->max_time) break; 47 th->stage_time = ts->ptime + th->Theta*ts->time_step; 48 th->shift = 1./(th->Theta*ts->time_step); 49 ts->ptime += ts->time_step; 50 51 ierr = VecCopy(sol,th->Xold);CHKERRQ(ierr); /* Used within function evalutaion */ 52 ierr = VecWAXPY(th->X,1./th->shift,th->Xdot,sol);CHKERRQ(ierr); 53 ierr = SNESSolve(ts->snes,PETSC_NULL,th->X);CHKERRQ(ierr); 54 ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); 55 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 56 ts->nonlinear_its += its; ts->linear_its += lits; 57 ierr = VecAXPY(sol,ts->time_step,th->Xdot);CHKERRQ(ierr); 58 ts->steps++; 59 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 60 } 61 62 *steps += ts->steps; 63 *ptime = ts->ptime; 64 PetscFunctionReturn(0); 65 } 66 67 /*------------------------------------------------------------*/ 68 #undef __FUNCT__ 69 #define __FUNCT__ "TSDestroy_Theta" 70 static PetscErrorCode TSDestroy_Theta(TS ts) 71 { 72 TS_Theta *th = (TS_Theta*)ts->data; 73 PetscErrorCode ierr; 74 75 PetscFunctionBegin; 76 ierr = VecDestroy(th->Xold);CHKERRQ(ierr); 77 ierr = VecDestroy(th->X);CHKERRQ(ierr); 78 ierr = VecDestroy(th->Xdot);CHKERRQ(ierr); 79 ierr = VecDestroy(th->res);CHKERRQ(ierr); 80 ierr = PetscFree(th);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 /* 85 This defines the nonlinear equation that is to be solved with SNES 86 G(U) = F[t0+T*dt, U, (U-U0)*shift] = 0 87 */ 88 #undef __FUNCT__ 89 #define __FUNCT__ "TSThetaFunction" 90 static PetscErrorCode TSThetaFunction(SNES snes,Vec x,Vec y,void *ctx) 91 { 92 TS ts = (TS)ctx; 93 TS_Theta *th = (TS_Theta*)ts->data; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,th->Xold,x);CHKERRQ(ierr); 98 ierr = TSComputeIFunction(ts,th->stage_time,x,th->Xdot,y);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "TSThetaJacobian" 104 static PetscErrorCode TSThetaJacobian(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,void *ctx) 105 { 106 TS ts = (TS)ctx; 107 TS_Theta *th = (TS_Theta*)ts->data; 108 PetscErrorCode ierr; 109 110 PetscFunctionBegin; 111 /* th->Xdot will have already been computed in TSThetaFunction */ 112 ierr = TSComputeIJacobian(ts,th->stage_time,x,th->Xdot,th->shift,A,B,str);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "TSSetUp_Theta" 119 static PetscErrorCode TSSetUp_Theta(TS ts) 120 { 121 TS_Theta *th = (TS_Theta*)ts->data; 122 PetscErrorCode ierr; 123 124 PetscFunctionBegin; 125 ierr = VecDuplicate(ts->vec_sol,&th->Xold);CHKERRQ(ierr); 126 ierr = VecDuplicate(ts->vec_sol,&th->X);CHKERRQ(ierr); 127 ierr = VecDuplicate(ts->vec_sol,&th->Xdot);CHKERRQ(ierr); 128 ierr = VecDuplicate(ts->vec_sol,&th->res);CHKERRQ(ierr); 129 ierr = SNESSetFunction(ts->snes,th->res,TSThetaFunction,ts);CHKERRQ(ierr); 130 ierr = SNESSetJacobian(ts->snes,ts->A,ts->B,TSThetaJacobian,ts);CHKERRQ(ierr); 131 PetscFunctionReturn(0); 132 } 133 /*------------------------------------------------------------*/ 134 135 #undef __FUNCT__ 136 #define __FUNCT__ "TSSetFromOptions_Theta" 137 static PetscErrorCode TSSetFromOptions_Theta(TS ts) 138 { 139 TS_Theta *th = (TS_Theta*)ts->data; 140 PetscErrorCode ierr; 141 142 PetscFunctionBegin; 143 ierr = PetscOptionsHead("SUNDIALS ODE solver options");CHKERRQ(ierr); 144 { 145 ierr = PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,PETSC_NULL);CHKERRQ(ierr); 146 } 147 ierr = PetscOptionsTail();CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "TSView_Theta" 153 static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer) 154 { 155 TS_Theta *th = (TS_Theta*)ts->data; 156 PetscTruth iascii; 157 PetscErrorCode ierr; 158 159 PetscFunctionBegin; 160 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 161 if (iascii) { 162 ierr = PetscViewerASCIIPrintf(viewer," Theta=%G\n",th->Theta);CHKERRQ(ierr); 163 } else { 164 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for TS_Theta",((PetscObject)viewer)->type_name); 165 } 166 PetscFunctionReturn(0); 167 } 168 169 /* ------------------------------------------------------------ */ 170 /*MC 171 TS_Theta - DAE solver using the implicit Theta method 172 173 Level: beginner 174 175 .seealso: TSCreate(), TS, TSSetType() 176 177 M*/ 178 EXTERN_C_BEGIN 179 #undef __FUNCT__ 180 #define __FUNCT__ "TSCreate_Theta" 181 PetscErrorCode PETSCTS_DLLEXPORT TSCreate_Theta(TS ts) 182 { 183 TS_Theta *th; 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 ts->ops->destroy = TSDestroy_Theta; 188 ts->ops->view = TSView_Theta; 189 ts->ops->setup = TSSetUp_Theta; 190 ts->ops->step = TSStep_Theta; 191 ts->ops->setfromoptions = TSSetFromOptions_Theta; 192 193 ierr = SNESCreate(((PetscObject)ts)->comm,&ts->snes);CHKERRQ(ierr); 194 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 195 196 ierr = PetscNewLog(ts,TS_Theta,&th);CHKERRQ(ierr); 197 ts->data = (void*)th; 198 199 th->Theta = 0.5; 200 201 PetscFunctionReturn(0); 202 } 203 EXTERN_C_END 204