Actual source code: ex7.c
petsc-3.6.4 2016-04-12
1: static char help[] = "Solves DAE with integrator only on non-algebraic terms \n";
3: #include <petscts.h>
5: /*
6: \dot{U} = f(U,V)
7: F(U,V) = 0
9: Same as ex6.c except the user provided functions take input values as a single vector instead of two vectors
10: */
15: /*
16: f(U,V) = U + V
18: */
19: PetscErrorCode f(PetscReal t,Vec UV,Vec F)
20: {
21: PetscErrorCode ierr;
22: const PetscScalar *u,*v;
23: PetscScalar *f;
24: PetscInt n,i;
27: VecGetLocalSize(UV,&n);
28: n = n/2;
29: VecGetArrayRead(UV,&u);
30: v = u + n;
31: VecGetArray(F,&f);
32: for (i=0; i<n; i++) f[i] = u[i] + v[i];
33: VecRestoreArrayRead(UV,&u);
34: VecRestoreArray(F,&f);
35: return(0);
36: }
40: /*
41: F(U,V) = U - V
43: */
44: PetscErrorCode F(PetscReal t,Vec UV,Vec F)
45: {
46: PetscErrorCode ierr;
47: const PetscScalar *u,*v;
48: PetscScalar *f;
49: PetscInt n,i;
52: VecGetLocalSize(UV,&n);
53: n = n/2;
54: VecGetArrayRead(UV,&u);
55: v = u + n;
56: VecGetArray(F,&f);
57: for (i=0; i<n; i++) f[i] = u[i] - v[i];
58: VecRestoreArrayRead(UV,&u);
59: VecRestoreArray(F,&f);
60: return(0);
61: }
63: typedef struct {
64: PetscReal t;
65: SNES snes;
66: Vec UV,V;
67: VecScatter scatterU,scatterV;
68: PetscErrorCode (*f)(PetscReal,Vec,Vec);
69: PetscErrorCode (*F)(PetscReal,Vec,Vec);
70: } AppCtx;
72: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
73: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);
77: int main(int argc,char **argv)
78: {
80: AppCtx ctx;
81: TS ts;
82: Vec tsrhs,U;
83: IS is;
84: PetscInt I;
85: PetscMPIInt rank;
87: PetscInitialize(&argc,&argv,(char*)0,help);
88: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
89: TSCreate(PETSC_COMM_WORLD,&ts);
90: TSSetProblemType(ts,TS_NONLINEAR);
91: TSSetType(ts,TSEULER);
92: TSSetFromOptions(ts);
93: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&tsrhs);
94: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
95: TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
96: ctx.f = f;
98: SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
99: SNESSetFromOptions(ctx.snes);
100: SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
101: SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
102: ctx.F = F;
103: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);
105: /* Create scatters to move between separate U and V representation and UV representation of solution */
106: VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&ctx.UV);
107: I = 2*rank;
108: ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
109: VecScatterCreate(U,NULL,ctx.UV,is,&ctx.scatterU);
110: ISDestroy(&is);
111: I = 2*rank + 1;
112: ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
113: VecScatterCreate(ctx.V,NULL,ctx.UV,is,&ctx.scatterV);
114: ISDestroy(&is);
116: VecSet(U,1.0);
117: TSSolve(ts,U);
119: VecDestroy(&ctx.V);
120: VecDestroy(&ctx.UV);
121: VecScatterDestroy(&ctx.scatterU);
122: VecScatterDestroy(&ctx.scatterV);
123: VecDestroy(&tsrhs);
124: VecDestroy(&U);
125: SNESDestroy(&ctx.snes);
126: TSDestroy(&ts);
127: PetscFinalize();
128: return 0;
129: }
133: /*
134: Defines the RHS function that is passed to the time-integrator.
136: Solves F(U,V) for V and then computes f(U,V)
138: */
139: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
140: {
141: AppCtx *ctx = (AppCtx*)actx;
145: ctx->t = t;
146: VecScatterBegin(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
147: VecScatterEnd(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
148: SNESSolve(ctx->snes,NULL,ctx->V);
149: VecScatterBegin(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
150: VecScatterEnd(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
151: (*ctx->f)(t,ctx->UV,F);
152: return(0);
153: }
157: /*
158: Defines the nonlinear function that is passed to the nonlinear solver
160: */
161: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
162: {
163: AppCtx *ctx = (AppCtx*)actx;
167: VecScatterBegin(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
168: VecScatterEnd(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
169: (*ctx->F)(ctx->t,ctx->UV,F);
170: return(0);
171: }