Actual source code: tsfwdsen.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/tsimpl.h>
3: PetscLogEvent TS_ForwardStep;
5: /*@
6: TSForwardSetUp - Sets up the internal data structures for the later use
7: of forward sensitivity analysis
9: Collective on TS
11: Input Parameter:
12: . ts - the TS context obtained from TSCreate()
14: Level: advanced
16: .keywords: TS, forward sensitivity, setup
18: .seealso: TSCreate(), TSDestroy(), TSSetUp()
19: @*/
20: PetscErrorCode TSForwardSetUp(TS ts)
21: {
25: if (ts->forwardsetupcalled) return(0);
27: if (ts->vec_costintegral && !ts->vecs_integral_sensi && !ts->vecs_integral_sensip ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSForwardSetIntegralGradients() before TSSetCostIntegrand()");
29: if (ts->vecs_integral_sensi || ts->vecs_integral_sensip) {
30: VecDuplicateVecs(ts->vec_sol,ts->numcost,&ts->vecs_drdy);
31: }
32: if (ts->vecs_integral_sensip) {
33: VecDuplicateVecs(ts->vecs_integral_sensip[0],ts->numcost,&ts->vecs_drdp);
34: }
35: if (ts->ops->forwardsetup) {
36: (*ts->ops->forwardsetup)(ts);
37: }
38: ts->forwardsetupcalled = PETSC_TRUE;
39: return(0);
40: }
42: /*@C
43: TSForwardSetRHSJacobianP - Sets the function that computes the Jacobian of G w.r.t. the parameters p where y_t = G(y,p,t), as well as the location to store the vector array.
45: Logically Collective on TS
47: Input Parameters:
48: + ts - The TS context obtained from TSCreate()
49: - func - The function
51: Calling sequence of func:
52: $ func (TS ts,PetscReal t,Vec y,Vec* a,void *ctx);
53: + t - current timestep
54: . y - input vector (current ODE solution)
55: . a - output vector array
56: - ctx - [optional] user-defined function context
58: Level: intermediate
60: Notes: the number of vectors in a is the same as the number of parameters and each vector is of the same size as the system dimension.
62: .keywords: TS, forward sensitivity
64: .seealso: TSForwardSetSensitivities(), TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
65: @*/
66: PetscErrorCode TSForwardSetRHSJacobianP(TS ts,Vec* a,PetscErrorCode (*func)(TS,PetscReal,Vec,Vec*,void*),void *ctx)
67: {
72: ts->vecsrhsjacobianp = func;
73: ts->vecsrhsjacobianpctx = ctx;
74: if(a) ts->vecs_jacp = a;
75: return(0);
76: }
78: /*@C
79: TSForwardComputeRHSJacobianP - Runs the user-defined JacobianP function.
81: Collective on TS
83: Input Parameters:
84: . ts - The TS context obtained from TSCreate()
86: Level: developer
88: .keywords: TS, forward sensitivity
90: .seealso: TSForwardSetRHSJacobianP()
91: @*/
92: PetscErrorCode TSForwardComputeRHSJacobianP(TS ts,PetscReal t,Vec X,Vec* A )
93: {
100: if (ts->vecsrhsjacobianp) {
101: PetscStackPush("TS user JacobianP function for sensitivity analysis");
102: (*ts->vecsrhsjacobianp)(ts,t,X,A,ts->vecsrhsjacobianpctx);
103: PetscStackPop;
104: }
105: return(0);
106: }
108: /*@
109: TSForwardSetIntegralGradients - Set the vectors holding forward sensitivities of the integral term.
111: Input Parameter:
112: . ts- the TS context obtained from TSCreate()
113: . numfwdint- number of integrals
114: . v = the vectors containing the gradients for each integral wrt initial values
115: . vp = the vectors containing the gradients for each integral wrt parameters
117: Level: intermediate
119: .keywords: TS, forward sensitivity
121: .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
122: @*/
123: PetscErrorCode TSForwardSetIntegralGradients(TS ts,PetscInt numfwdint,Vec *vp,Vec *v)
124: {
127: if (ts->numcost && ts->numcost!=numfwdint) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2rd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand()");
128: if (!ts->numcost) ts->numcost = numfwdint;
130: ts->vecs_integral_sensi = v;
131: ts->vecs_integral_sensip = vp;
132: return(0);
133: }
135: /*@
136: TSForwardGetIntegralGradients - Returns the forward sensitivities ofthe integral term.
138: Input Parameter:
139: . ts- the TS context obtained from TSCreate()
141: Output Parameter:
142: . v = the vectors containing the gradients for each integral wrt initial values
143: . vp = the vectors containing the gradients for each integral wrt parameters
145: Level: intermediate
147: .keywords: TS, forward sensitivity
149: .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
150: @*/
151: PetscErrorCode TSForwardGetIntegralGradients(TS ts,PetscInt *numfwdint,Vec **v,Vec **vp)
152: {
156: if (numfwdint) *numfwdint = ts->numcost;
157: if (v) *v = ts->vecs_integral_sensi;
158: if (vp) *vp = ts->vecs_integral_sensip;
159: return(0);
160: }
162: /*@
163: TSForwardStep - Compute the forward sensitivity for one time step.
165: Collective on TS
167: Input Arguments:
168: . ts - time stepping context
170: Level: advanced
172: Notes:
173: This function cannot be called until TSStep() has been completed.
175: .keywords: TS, forward sensitivity
177: .seealso: TSForwardSetSensitivities(), TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardSetUp()
178: @*/
179: PetscErrorCode TSForwardStep(TS ts)
180: {
183: if (!ts->ops->forwardstep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide forward sensitivity analysis",((PetscObject)ts)->type_name);
184: PetscLogEventBegin(TS_ForwardStep,ts,0,0,0);
185: (*ts->ops->forwardstep)(ts);
186: PetscLogEventEnd(TS_ForwardStep,ts,0,0,0);
187: return(0);
188: }
190: /*@
191: TSForwardSetSensitivities - Sets the initial value of the trajectory sensitivities of solution w.r.t. the problem parameters and initial values.
193: Logically Collective on TS and Vec
195: Input Parameters:
196: + ts - the TS context obtained from TSCreate()
197: . nump - number of parameters
198: . sp - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
199: . num - number of initial values
200: - s - sensitivities with respect to the (selected) initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
202: Level: beginner
204: Notes:
205: Forward sensitivity is also called 'trajectory sensitivity' in some fields such as power systems.
206: This function turns on a flag to trigger TSSolve() to compute forward sensitivities automatically.
207: You must call this function before TSSolve().
208: The entries in these vectors must be correctly initialized with the values s_i = dy/dp|startingtime.
209: The two user-provided sensitivity vector arrays will be packed into one big array to simplify implementation.
211: .keywords: TS, timestep, set, forward sensitivity, initial values
213: .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
214: @*/
215: PetscErrorCode TSForwardSetSensitivities(TS ts,PetscInt nump,Vec *sp,PetscInt num,Vec *s)
216: {
217: PetscInt i;
222: ts->forward_solve = PETSC_TRUE;
223: ts->num_parameters = sp ? nump:0;
224: ts->num_initialvalues = s ? num:0;
225: /* pack fwdsensi and fwdsensip into a big array */
226: if (!ts->vecs_fwdsensipacked) {
227: PetscMalloc1(num+nump,&ts->vecs_fwdsensipacked);
228: }
229: for (i=0; i<num; i++) ts->vecs_fwdsensipacked[i] = s[i];
230: for (i=0; i<nump; i++) ts->vecs_fwdsensipacked[i+num] = sp[i];
231: return(0);
232: }
234: /*@
235: TSForwardGetSensitivities - Returns the trajectory sensitivities
237: Not Collective, but Vec returned is parallel if TS is parallel
239: Output Parameter:
240: + ts - the TS context obtained from TSCreate()
241: . nump - number of parameters
242: . sp - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
243: . num - number of initial values
244: - s - sensitivities with respect to the (selected) initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
246: Level: intermediate
248: .keywords: TS, forward sensitivity
250: .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
251: @*/
252: PetscErrorCode TSForwardGetSensitivities(TS ts,PetscInt *nump,Vec **sp,PetscInt *num,Vec **s)
253: {
255: if (nump) *nump = ts->num_parameters;
256: if (num) *num = ts->num_initialvalues;
257: if (sp) *sp = &ts->vecs_fwdsensipacked[(*num)];
258: if (s) *s = ts->vecs_fwdsensipacked;
259: return(0);
260: }