Actual source code: gts.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: gts.c,v 1.9 2000/01/10 03:54:19 knepley Exp $";
3: #endif
5: #include src/ts/tsimpl.h
6: #include src/gvec/gvecimpl.h
7: #include gsolver.h
9: /* Logging support */
10: int GTS_Reform, GTS_Reallocate;
12: /*
13: This file provides routines for grid TS objects. They support solution of
14: systems of time-dependent nonlinear equations generated from a mesh and
15: discretization specified by a Grid object.
16: */
18: #undef __FUNCT__
20: /*@C
21: GTSDestroy - Destroys a grid TS.
23: Collective on GTS
25: Input Parameters:
26: . ts - The nonlinear solver context
28: Level: beginner
30: .keywords: grid TS, destroy
31: .seealso: TSDestroy(), GTSDuplicate()
32: @*/
33: int GTSDestroy(GTS ts)
34: {
35: TSProblemType type;
36: int ierr;
40: if (--ts->refct > 0)
41: return(0);
42: TSGetProblemType(ts, &type);
43: if (type == TS_LINEAR) {
44: GMatDestroy(ts->A);
45: }
46: PetscFree(ts->isExplicit);
47: PetscFree(ts->Iindex);
48: if (ts->snes) {
49: GSNESDestroy(ts->snes);
50: ts->snes = PETSC_NULL;
51: }
52: TSDestroy(ts);
53: return(0);
54: }
56: #undef __FUNCT__
58: /*@
59: GTSView - Views a grid TS.
61: Collective on GTS
63: Input Parameters:
64: . ts - The grid TS
65: . viewer - An optional visualization context
67: Options Database Key:
68: $ -ts_view : calls TSView() at end of TSSolve()
70: Notes:
71: The available visualization contexts include
72: $ VIEWER_STDOUT_SELF - standard output (default)
73: $ VIEWER_STDOUT_WORLD - synchronized standard
74: $ output where only the first processor opens
75: $ the file. All other processors send their
76: $ data to the first processor to print.
78: The user can open alternative vistualization contexts with
79: $ PetscViewerFileOpenASCII() - output vector to a specified file
81: Level: beginner
83: .keywords: view, visualize, output, print, write, draw
84: .seealso: TSView()
85: @*/
86: int GTSView(GTS ts, PetscViewer viewer)
87: {
88: Grid grid;
89: FILE *fd;
90: int numActiveFields;
91: int field, func, op;
92: PetscTruth isascii;
93: int ierr;
97: if (!viewer) {
98: viewer = PETSC_VIEWER_STDOUT_SELF;
99: } else {
101: }
103: GTSGetGrid(ts, &grid);
104: PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &isascii);
105: if (isascii == PETSC_TRUE) {
106: GridView(grid, viewer);
107: PetscViewerFlush(viewer);
108: PetscViewerASCIIGetPointer(viewer, &fd);
109: PetscFPrintf(ts->comm, fd, "GTS Object:n");
110: GridGetNumActiveFields(grid, &numActiveFields);
111: if (numActiveFields == 1) {
112: PetscFPrintf(ts->comm, fd, " %d active field:n", numActiveFields);
113: } else {
114: PetscFPrintf(ts->comm, fd, " %d active fields:n", numActiveFields);
115: }
116: for(field = 0; field < grid->numFields; field++) {
117: if (grid->fields[field].isActive == PETSC_FALSE) continue;
118: if (grid->fields[field].name) {
119: PetscFPrintf(grid->comm, fd, " %s field with %d components:n", grid->fields[field].name, grid->fields[field].numComp);
120: } else {
121: PetscFPrintf(grid->comm, fd, " field %d with %d components:n", field, grid->fields[field].numComp);
122: }
123: if (ts->isExplicit[field]) PetscFPrintf(ts->comm, fd, " Explicitly time dependentn");
124: for(func = 0; func < grid->numRhsFuncs; func++) {
125: if (grid->rhsFuncs[func].field == field) {
126: PetscFPrintf(grid->comm, fd, " Rhs function with scalar multiplier %gn", PetscRealPart(grid->rhsFuncs[func].alpha));
127: }
128: }
129: for(op = 0; op < grid->numRhsOps; op++) {
130: if (grid->rhsOps[op].field == field) {
131: if (grid->rhsOps[op].nonlinearOp != PETSC_NULL) {
132: PetscFPrintf(grid->comm, fd, " Rhs nonlinear operator with scalar multiplier %gn",
133: PetscRealPart(grid->rhsOps[field].alpha));
134: } else {
135: PetscFPrintf(grid->comm, fd, " Rhs operator %d with scalar multiplier %gn",
136: grid->rhsOps[op].op, PetscRealPart(grid->rhsOps[op].alpha));
137: }
138: }
139: }
140: for(op = 0; op < grid->numMatOps; op++) {
141: if (grid->matOps[op].field == field) {
142: PetscFPrintf(grid->comm, fd, " Jacobian operator %d with scalar multiplier %gn",
143: grid->matOps[op].op, PetscRealPart(grid->matOps[op].alpha));
144: }
145: }
146: }
147: TSView(ts, viewer);
148: }
150: return(0);
151: }
153: #undef __FUNCT__
155: /*@
156: GTSSerialize - This function stores or recreates a grid TS using a viewer for a binary file.
158: Collective on grid
160: Input Parameters:
161: . comm - The communicator for the grid TS object
162: . viewer - The viewer context
163: . store - This flag is PETSC_TRUE is data is being written, otherwise it will be read
165: Output Parameter:
166: . ts - The grid TS
168: Level: beginner
170: .keywords: grid TS, serialize
171: .seealso: GridSerialize()
172: @*/
173: int GTSSerialize(Grid grid, GTS *ts, PetscViewer viewer, PetscTruth store)
174: {
175: MPI_Comm comm;
176: int ierr;
182: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
183: GSolverInitializePackage(PETSC_NULL);
184: #endif
187: PetscObjectGetComm((PetscObject) grid, &comm);
188: /* We need a way to communicate the grid to the TS */
189: if (store == PETSC_FALSE) {
190: *ts = (GTS) grid;
191: }
192: TSSerialize(comm, ts, viewer, store);
193: /* Newly created object should have grid as a parent */
194: if (store == PETSC_FALSE) {
195: PetscObjectCompose((PetscObject) *ts, "Grid", (PetscObject) grid);
197: /* Set boundary conditions */
198: TSSetRhsBC(*ts, GTSRhsBC);
199: TSSetSolutionBC(*ts, GTSSolutionBC);
200: /* Save space/time by reducing the system at the element level */
201: if ((*ts)->problem_type == TS_NONLINEAR) {
202: GridSetReduceElement(grid, PETSC_TRUE);
203: }
205: /* Set update functions */
206: TSSetPreStep(*ts, GTSPreStep);
207: TSSetUpdate(*ts, GTSUpdate);
208: TSSetPostStep(*ts, GTSPostStep);
209: }
211: return(0);
212: }
214: #undef __FUNCT__
216: /*@C
217: GTSDuplicate - Duplicates a grid TS.
219: Collective on GTS
221: Input Parameter:
222: . ts - The grid TS
224: Output Parameter:
225: . newts - The new grid TS
227: Level: beginner
229: .keywords: grid ts, duplicate
230: .seealso: TSDuplicate(), GTSDestroy()
231: @*/
232: int GTSDuplicate(GTS ts, GTS *newts)
233: {
236: SETERRQ(PETSC_ERR_SUP, " ");
237: #if 0
238: return(0);
239: #endif
240: }
242: #undef __FUNCT__
244: /*@C
245: GTSReform - Cleans up any auxilliary structures in the GTS. Should
246: be called after GridReform(), but before any user specific changes
247: to the grid. Usually followed by GTSReallocate().
249: Collective on GTS
251: Input Parameter:
252: . ts - The grid TS
254: Level: advanced
256: .keywords: grid ts, reform
257: .seealso: GTSReallocate(), GridReform()
258: @*/
259: int GTSReform(GTS ts)
260: {
265: if (ts->ops->reform) {
266: PetscLogEventBegin(GTS_Reform, ts, 0, 0, 0);
267: (*ts->ops->reform)(ts);
268: PetscLogEventEnd(GTS_Reform, ts, 0, 0, 0);
269: }
270: return(0);
271: }
273: #undef __FUNCT__
275: /*@C
276: GTSReallocate - Reallocates the storage for a grid TS. Should be
277: called after GTSReform() and any user updating of the grid.
279: Collective on GTS
281: Input Parameters:
282: . ts - The grid TS
283: . x - The new solution vector
285: Level: advanced
287: .keywords: grid ts, reallocate
288: .seealso: GTSReform(), GridReform()
289: @*/
290: int GTSReallocate(GTS ts, GVec x)
291: {
297: TSSetSolution(ts, x);
298: if (ts->ops->reallocate) {
299: PetscLogEventBegin(GTS_Reallocate, ts, 0, 0, 0);
300: (*ts->ops->reallocate)(ts);
301: PetscLogEventEnd(GTS_Reallocate, ts, 0, 0, 0);
302: }
303: return(0);
304: }
306: #undef __FUNCT__
308: /*@
309: GTSGetGrid - Gets the grid from a grid TS.
311: Not collective
313: Input Parameter:
314: . ts - The grid TS
316: Output Parameter:
317: . grid - The grid context
319: Level: intermediate
321: .keywords: grid ts, get
322: .seealso: GVecGetGrid(), GMatGetGrid(), GSNESGetGrid()
323: @*/
324: int GTSGetGrid(GTS ts, Grid *grid)
325: {
332: PetscObjectQuery((PetscObject) ts, "Grid", (PetscObject *) grid);
334: return(0);
335: }
337: #undef __FUNCT__
339: /*@
340: GTSGetInitialTimeStep - Returns the initial time step
342: Not collective
344: Input Parameter:
345: . ts - The grid TS
347: Output Parameter:
348: . dt - The initial time step
350: Level: intermediate
352: .keywords: grid ts, get, time step
353: .seealso: TSSetInitialTimeStep()
354: @*/
355: int GTSGetInitialTimeStep(GTS ts, double *dt)
356: {
360: *dt = ts->initial_time_step;
361: return(0);
362: }
364: #undef __FUNCT__
366: /*@
367: GTSGetTimeDependence - Determines whether or not a field has
368: explicit time dependence.
370: Not collective
372: Input Parameters:
373: . ts - The grid TS
374: . field - The field
376: Output Parameter:
377: . flag - Indicates whether or not the field has explicit time dependence
379: Level: intermediate
381: .keywords: grid ts, get
382: .seealso: GTSSetTimeDependence()
383: @*/
384: int GTSGetTimeDependence(GTS ts, int field, PetscTruth *flag)
385: {
386: Grid grid;
387: int numFields;
388: int ierr;
393: GTSGetGrid(ts, &grid);
395: GridGetNumFields(grid, &numFields);
396: if ((field < 0) || (field > numFields)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid field number");
397: *flag = ts->isExplicit[field];
398: return(0);
399: }
401: #undef __FUNCT__
403: /*@
404: GTSSetContext - Sets a user-defined context for use in calculation.
406: Collective on GTS
408: Input Parameters:
409: . ts - The grid TS
410: . ctx - The context
412: Level: intermediate
414: .keywords: grid ts, context
415: .seealso: GTSGetContext()
416: @*/
417: int GTSSetContext(GTS ts, void *ctx)
418: {
421: ts->funP = ctx;
422: ts->jacP = ctx;
423: return(0);
424: }
426: #undef __FUNCT__
428: /*@
429: GTSSetTimeDependence - Sets the explicit time dependence of a field.
431: Collective on GTS
433: Input Parameters:
434: . ts - The grid TS
435: . field - The field
436: . flag - Indicates whether or not the field has explicit time dependence
438: Level: intermediate
440: .keywords: grid ts, get
441: .seealso: GTSGetTimeDependence()
442: @*/
443: int GTSSetTimeDependence(GTS ts, int field, PetscTruth flag)
444: {
445: Grid grid;
446: int numFields;
447: int ierr;
451: GTSGetGrid(ts, &grid);
453: GridGetNumFields(grid, &numFields);
454: if ((field < 0) || (field > numFields)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid field number");
455: ts->isExplicit[field] = flag;
456: return(0);
457: }
459: #undef __FUNCT__
461: /*@
462: GTSEvaluateRhs - Constructs the Rhs vector for U_t = F(U, t).
464: Collective on GTS
466: Input Parameters:
467: . ts - The grid TS
468: . x - The current iterate
469: . ctx - The optional user context
471: Output Parameter:
472: . f - The function value
474: Level: advanced
476: .keywords: grid ts, rhs
477: .seealso: GSNESEvaluateRhs()
478: @*/
479: int GTSEvaluateRhs(GTS ts, double t, GVec x, GVec f, PetscObject ctx)
480: {
481: Grid grid;
482: PetscObject oldCtx;
483: PetscScalar zero = 0.0;
484: int ierr;
488: GTSGetGrid(ts, &grid);
490: /* Initialize vector */
491: VecSet(&zero, f);
492: /* Setup new context */
493: GTSCreateContext(ts, t, ctx, &oldCtx);
494: /* Form function */
495: GridEvaluateRhs(grid, x, f, ctx);
496: /* Cleanup */
497: GTSDestroyContext(ts, ctx, oldCtx);
498: return(0);
499: }
501: #undef __FUNCT__
503: /*@
504: GTSEvaluateJacobian - Constructs the Jacobian matrix for U_t = F(U,t)
506: Collective on GTS
508: Input Parameters:
509: . ts - The grid TS
510: . x - The current iterate
511: . flag - The matrix structure flag
512: . ctx - The optional user context
514: Output Parameters:
515: . J - The Jacobian matrix
516: . M - The preconditioner for the Jacobian matrix, usually the same as J
518: Level: advanced
520: .keywords: grid ts, evaluate, jacobian
521: .seealso: SNESComputeJacobian(), GridEvaluateSystemMatrix()
522: @*/
523: int GTSEvaluateJacobian(GTS ts, double t, GVec x, GMat *J, GMat *M, MatStructure *flag, PetscObject ctx)
524: {
525: Grid grid;
526: PetscObject oldCtx;
527: int ierr;
531: GTSGetGrid(ts, &grid);
533: /* Initialize matrix */
534: MatZeroEntries(*J);
535: /* Setup new context */
536: GTSCreateContext(ts, t, ctx, &oldCtx);
537: /* Form function */
538: GridEvaluateSystemMatrix(grid, x, J, M, flag, ctx);
539: /* Apply boundary conditions */
540: GMatSetBoundary(*J, 1.0, ctx);
541: /* Cleanup */
542: GTSDestroyContext(ts, ctx, oldCtx);
543: return(0);
544: }
546: #undef __FUNCT__
548: /*@
549: GTSEvaluateSystemMatrix - Constructs the system matrix for U_t = A(U,t) U
551: Collective on GTS
553: Input Parameters:
554: . ts - The grid TS
555: . flag - The matrix structure flag
556: . ctx - The optional user context
558: Output Parameters:
559: . A - The system matrix
560: . M - The preconditioner for the system matrix, usually the same as A
562: Level: advanced
564: .keywords: grid ts, evaluate, jacobian
565: .seealso: SNESComputeJacobian(), GridEvaluateSystemMatrix()
566: @*/
567: int GTSEvaluateSystemMatrix(GTS ts, double t, GMat *A, GMat *M, MatStructure *flag, PetscObject ctx)
568: {
569: Grid grid;
570: PetscObject oldCtx;
571: int ierr;
575: GTSGetGrid(ts, &grid);
577: /* Initialize matrix */
578: MatZeroEntries(*A);
579: /* Setup new context */
580: GTSCreateContext(ts, t, ctx, &oldCtx);
581: /* Form function */
582: GridEvaluateSystemMatrix(grid, PETSC_NULL, A, M, flag, ctx);
583: /* Apply boundary conditions */
584: GMatSetBoundary(*A, 1.0, ctx);
585: /* Cleanup */
586: GTSDestroyContext(ts, ctx, oldCtx);
587: return(0);
588: }
590: #undef __FUNCT__
592: /*@
593: GTSCalcBCValues - This function calculates the boundary values. It
594: is normally called once a timestep when using time dependent boundary
595: conditions.
597: Collective on GTS
599: Input Parameter:
600: . ts - The GTS
602: Level: advanced
604: .keywords: grid ts, reduction, boundary conditions
605: .seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
606: @*/
607: int GTSCalcBCValues(GTS ts)
608: {
609: Grid grid;
610: PetscObject oldCtx;
611: int ierr;
615: GTSGetGrid(ts, &grid);
616: /* Setup new context */
617: GTSCreateContext(ts, ts->ptime, (PetscObject) ts->funP, &oldCtx);
618: /* Form function */
619: GridCalcBCValues(grid, PETSC_TRUE, ts->funP);
620: /* Cleanup */
621: GTSDestroyContext(ts, (PetscObject) ts->funP, oldCtx);
622: return(0);
623: }
625: #undef __FUNCT__
627: /*@
628: GTSRhsBC - This sets the boundary conditions registered with the grid
629: on the Rhs vector, or zero boundary conditions for nonlinear systems.
631: Collective on GTS
633: Input Parameters:
634: . ts - The TS context obtained from TSCreate()
635: . ctx - The user-supplied context
637: Output Parameter:
638: . rhs - The Rhs
640: Level: advanced
642: .keywords: grid ts, timestep, set, boundary conditions
643: .seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
644: @*/
645: int GTSRhsBC(GTS ts, GVec rhs, void *ctx)
646: {
647: TSProblemType type;
648: PetscObject oldCtx;
649: int ierr;
654: TSGetProblemType(ts, &type);
655: /* Setup new context */
656: GTSCreateContext(ts, ts->ptime, (PetscObject) ctx, &oldCtx);
657: switch(type)
658: {
659: case TS_LINEAR:
660: GVecSetBoundary(rhs, ctx);
661: break;
662: case TS_NONLINEAR:
663: GVecSetBoundaryZero(rhs, ctx);
664: /* TSGetSolution(ts, &u); */
665: /* GVecSetBoundaryDifference(rhs, u, >Sctx); */
666: break;
667: default:
668: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid TS problem type %d", type);
669: }
670: /* Cleanup */
671: GTSDestroyContext(ts, (PetscObject) ctx, oldCtx);
672: return(0);
673: }
675: #undef __FUNCT__
677: /*@
678: GTSSolutionBC - This sets the boundary conditions registered with the grid
679: on the solution vector.
681: Collective on GTS
683: Input Parameters:
684: . ts - The TS context obtained from TSCreate()
685: . ctx - The user-supplied context
687: Output Parameter:
688: . sol - The solution
690: Level: advanced
692: .keywords: grid ts, timestep, set, boundary conditions
693: .seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
694: @*/
695: int GTSSolutionBC(GTS ts, GVec sol, void *ctx)
696: {
697: PetscObject oldCtx;
698: int ierr;
703: /* Setup new context */
704: GTSCreateContext(ts, ts->ptime, (PetscObject) ctx, &oldCtx);
705: /* Apply solution BC */
706: GVecSetBoundary(sol, ctx);
707: /* Cleanup */
708: GTSDestroyContext(ts, (PetscObject) ctx, oldCtx);
709: return(0);
710: }
712: #undef __FUNCT__
714: /*@
715: GTSSolutionBCforGSNES - This functions allows a GSNES
716: imbedded in a GTS to correctly set time dependent boundary
717: conditions on the initial solution at each time step.
719: Collective on GTS
721: Input Parameters:
722: . ts - The TS context obtained from TSCreate()
723: . ctx - The user-supplied context
725: Output Parameter:
726: . sol - The solution
728: Level: advanced
730: .keywords: grid ts, timestep, set, boundary conditions
731: .seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
732: @*/
733: int GTSSolutionBCforGSNES(GSNES snes, GVec sol, void *ctx)
734: {
735: GTS ts = (GTS) ctx;
741: GTSSolutionBC(ts, sol, ts->funP);
742: return(0);
743: }
745: #undef __FUNCT__
747: /*@
748: GTSPreStep - This is a general purpose function which
749: is called once at the beginning of time stepping. It
750: currently performs context wrapping and grid manipulation.
752: Collective on GTS
754: Input Parameters:
755: . ts - The GTS context
757: Level: advanced
759: .keywords: grid ts, timestep
760: .seealso: GTSUpdate(), GTSPostStep()
761: @*/
762: int GTSPreStep(GTS ts)
763: {
764: Grid grid;
765: PetscScalar mdt = 1.0/ts->time_step;
766: PetscTruth isTimeDependent;
767: TSProblemType type;
768: int numFields;
769: int field, f;
770: int ierr;
775: GTSGetGrid(ts, &grid);
776: TSGetProblemType(ts, &type);
777: if (type == TS_NONLINEAR) {
778: /* Scale Rhs so we get u_t - F = 0 */
779: GridScaleRhs(grid, -1.0);
781: /* Add the identity operator and scale to get J_{sys} = - J_{F} where J_{F} is the given Jacobian of F */
782: GridScaleSystemMatrix(grid, -1.0);
783: }
785: /* Add the identity operator and scale to get J = Delta t^{-1} I + J_{sys} where J_{sys} is the system matrix */
786: GridGetNumActiveFields(grid, &numFields);
787: for(f = 0; f < numFields; f++) {
788: GridGetActiveField(grid, f, &field);
789: GTSGetTimeDependence(ts, field, &isTimeDependent);
790: if (isTimeDependent) {
791: GridAddMatOperator(grid, field, field, IDENTITY, mdt, PETSC_FALSE, &ts->Iindex[field]);
792: }
793: }
795: return(0);
796: }
798: #undef __FUNCT__
800: /*@
801: GTSUpdate - This is a general purpose function which
802: is called at the beginning of every time step. It
803: currently calculates the mesh acceleration and moves
804: the mesh.
806: Collective on GTS
808: Input Parameters:
809: . ts - The GTS context
810: . t - The current time
812: Output Parameter:
813: . dt - The current time step
815: Level: advanced
817: .keywords: grid ts, timestep
818: .seealso: GTSPreStep(), GTSPostStep()
819: @*/
820: int GTSUpdate(GTS ts, PetscReal t, PetscReal *dt)
821: {
822: Grid grid;
823: Mesh mesh;
824: MeshMover mover;
825: PetscTruth isMoving;
826: int ierr;
831: GTSGetGrid(ts, &grid);
832: GridGetMesh(grid, &mesh);
833: MeshGetMovement(mesh, &isMoving);
834: if (isMoving == PETSC_TRUE) {
835: MeshGetMover(mesh, &mover);
836: /* Calculate mesh acceleration */
837: MeshMoverCalcNodeAccelerations(mover, PETSC_TRUE);
838: /* Move the mesh */
839: MeshMoveMesh(mesh, ts->time_step);
840: }
841: return(0);
842: }
844: #undef __FUNCT__
846: /*@
847: GTSPostStep - This is a general purpose function which
848: is called once at the end of time stepping. It currently
849: performs context unwrapping and grid manipulation.
851: Collective on GTS
853: Input Parameters:
854: . ts - The GTS context
856: Level: advanced
858: .keywords: grid ts, timestep
859: .seealso: GTSPreStep(), GTSUpdate()
860: @*/
861: int GTSPostStep(GTS ts)
862: {
863: Grid grid;
864: TSProblemType type;
865: PetscTruth isTimeDependent;
866: int numFields;
867: int field, f;
868: int ierr;
873: GTSGetGrid(ts, &grid);
874: TSGetProblemType(ts, &type);
875: if (type == TS_NONLINEAR) {
876: /* Scale Rhs back to the original state */
877: GridScaleRhs(grid, -1.0);
879: /* Remove the identity and scale back */
880: GridScaleSystemMatrix(grid, -1.0);
881: }
883: GridGetNumActiveFields(grid, &numFields);
884: for(f = 0; f < numFields; f++) {
885: GridGetActiveField(grid, f, &field);
886: GTSGetTimeDependence(ts, field, &isTimeDependent);
887: if (isTimeDependent) {
888: GridRemoveMatOperator(grid, ts->Iindex[field]);
889: }
890: }
892: return(0);
893: }
895: #undef __FUNCT__
897: /*@
898: GTSCreate - Creates a grid TS associated with a
899: particular discretization context.
901: Collective on Grid
903: Input Parameter:
904: . grid - The discretization context
905: . ctx - The optional user context
907: Output Parameter:
908: . ts - The integrator context
910: Level: beginner
912: Note:
913: The TS type is TS_NONLINEAR if any nonlinear operators have
914: been activated on the grid.
916: .keywords: grid ts, create
917: .seealso: TSCreate(), TSSetType()
918: @*/
919: int GTSCreate(Grid grid, void *ctx, GTS *ts)
920: {
921: MPI_Comm comm;
922: TSProblemType type;
923: GTS gts;
924: GMat A;
925: int numFields, numRhsOps;
926: int field, op;
927: int ierr;
932: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
933: GSolverInitializePackage(PETSC_NULL);
934: #endif
936: /* Get TS type */
937: type = TS_LINEAR;
938: GridGetNumRhsOperators(grid, &numRhsOps);
939: for(op = 0; op < grid->numRhsOps; op++) {
940: if (grid->rhsOps[op].nonlinearOp != PETSC_NULL) {
941: type = TS_NONLINEAR;
942: break;
943: }
944: }
946: /* Create the TS */
947: PetscObjectGetComm((PetscObject) grid, &comm);
948: TSCreate(comm, >s);
949: TSSetProblemType(gts, type);
950: gts->isGTS = PETSC_TRUE;
951: gts->bops->destroy = (int (*)(PetscObject)) GTSDestroy;
952: gts->bops->view = (int (*)(PetscObject, PetscViewer)) GTSView;
954: /* Set boundary conditions */
955: TSSetRhsBC(gts, GTSRhsBC);
956: TSSetSolutionBC(gts, GTSSolutionBC);
957: /* Save space/time by reducing the system at the element level */
958: if (type == TS_NONLINEAR) {
959: GridSetReduceElement(grid, PETSC_TRUE);
960: }
962: /* Set update functions */
963: TSSetPreStep(gts, GTSPreStep);
964: TSSetUpdate(gts, GTSUpdate);
965: TSSetPostStep(gts, GTSPostStep);
967: /* Create the array which indicates time dependent fields */
968: GridGetNumFields(grid, &numFields);
969: PetscMalloc(numFields * sizeof(PetscTruth), >s->isExplicit);
970: PetscLogObjectMemory(gts, numFields * sizeof(PetscTruth));
971: for(field = 0; field < numFields; field++) {
972: gts->isExplicit[field] = PETSC_TRUE;
973: }
975: /* Create the array which holds the index of the identity operator for each time dependent field */
976: PetscMalloc(numFields * sizeof(int), >s->Iindex);
977: PetscLogObjectMemory(gts, numFields * sizeof(int));
978: for(field = 0; field < numFields; field++)
979: gts->Iindex[field] = -1;
981: /* Setup TS data structures */
982: gts->funP = ctx;
983: gts->jacP = ctx;
984: if (type == TS_LINEAR) {
985: PetscTruth reduceSystem;
987: GMatCreate(grid, &A);
988: GridGetReduceSystem(grid, &reduceSystem);
989: if (reduceSystem == PETSC_FALSE) {
990: MatSetOption(A, MAT_YES_NEW_NONZERO_LOCATIONS);
991: }
992: gts->A = A;
993: gts->B = A;
994: PetscLogObjectParent(grid, A);
995: }
997: PetscObjectCompose((PetscObject) gts, "Grid", (PetscObject) grid);
998: *ts = gts;
999: return(0);
1000: }
1002: #undef __FUNCT__
1004: /*@
1005: GTSSolutionMonitor - Monitors solution at each GTS iteration.
1007: Collective on GTS
1009: Input Parameters:
1010: . ts - The integrator context
1011: . step - The current time step
1012: . ltime - The current time
1013: . sol - The current solution vector
1014: . ctx - The viewer
1016: Level: intermediate
1018: .keywords: grid ts, monitor, solution
1019: .seealso: TSDefaultMonitor(),TSSetMonitor(),GTSResidualMonitor(),
1020: GTSErrorMonitor()
1021: @*/
1022: int GTSSolutionMonitor(GTS ts, int step, PetscReal ltime, GVec sol, void *ctx)
1023: {
1024: PetscViewer viewer = (PetscViewer) ctx;
1025: int ierr;
1030: GVecView(sol, viewer);
1031: return(0);
1032: }
1034: #undef __FUNCT__
1036: /*@
1037: GTSErrorMonitor - Displays the error at each iteration.
1039: Collective on GTS
1041: Input Parameters:
1042: + ts - The integrator context
1043: . step - The current time step
1044: . ltime - The current time
1045: . sol - The current solution vector
1046: - monCtx - The GTSErrorMonitorCxt
1048: Notes:
1049: The final argument to TSSetMonitor() with this routine must be
1050: a pointer to a GTSErrorMonitorCtx.
1052: Level: intermediate
1054: .keywords: grid ts, monitor, error
1055: .seealso: TSDefaultMonitor(),TSSetMonitor(),GTSSolutionMonitor(), GTSResidualMonitor()
1056: @*/
1057: int GTSErrorMonitor(GTS ts, int step, PetscReal ltime, GVec sol, void *monCtx)
1058: {
1059: GTSErrorMonitorCtx *errorCtx = (GTSErrorMonitorCtx *) monCtx;
1060: PetscObject oldCtx;
1061: PetscScalar minusOne = -1.0;
1062: GVec e;
1063: MPI_Comm comm;
1064: FILE *file;
1065: PetscReal norm_2, norm_max, sol_norm_2;
1066: int ierr;
1072: GTSCreateContext(ts, ltime, (PetscObject) errorCtx->ctx, &oldCtx);
1073: (*errorCtx->solutionFunc)(errorCtx->solution, errorCtx->ctx);
1074: GVecGetWorkGVec(errorCtx->solution, &e);
1075: VecWAXPY(&minusOne, sol, errorCtx->solution, e);
1076: GVecView(e, errorCtx->error_viewer);
1077: GTSDestroyContext(ts, (PetscObject) errorCtx->ctx, oldCtx);
1079: /* Compute 2-norm and max-norm of error */
1080: if (errorCtx->norm_error_viewer) {
1081: PetscObjectGetComm((PetscObject) ts, &comm);
1082: PetscViewerASCIIGetPointer(errorCtx->norm_error_viewer, &file);
1083: VecNorm(e, NORM_2, &norm_2);
1084: VecNorm(e, NORM_MAX, &norm_max);
1085: VecNorm(errorCtx->solution, NORM_2, &sol_norm_2);
1086: PetscFPrintf(comm, file, "Timestep %d time %g error 2-norm %g error max-norm %g exact sol 2-norm %gn",
1087: step, ltime, norm_2, norm_max, sol_norm_2);
1088: }
1089: GVecRestoreWorkGVec(errorCtx->solution, &e);
1090: return(0);
1091: }
1093: EXTERN_C_BEGIN
1094: #undef __FUNCT__
1096: int GTSOptionsChecker_Private(GTS ts)
1097: {
1098: char *prefix;
1099: int loc[4], nmax;
1100: MPI_Comm comm;
1101: PetscViewer viewer;
1102: PetscDraw draw;
1103: PetscTruth opt;
1104: int ierr;
1107: TSGetOptionsPrefix(ts, &prefix);
1108: PetscObjectGetComm((PetscObject) ts, &comm);
1110: nmax = 4;
1111: loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
1112: PetscOptionsGetIntArray(prefix, "-gvec_ts_solutionmonitor", loc, &nmax, &opt);
1113: if (opt) {
1114: PetscViewerDrawOpen(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);
1115: PetscViewerDrawGetDraw(viewer, 0, &draw);
1116: PetscDrawSetTitle(draw, "Approx. Solution");
1117: PetscLogObjectParent(ts, (PetscObject) viewer);
1118: TSSetMonitor(ts, GTSSolutionMonitor, (void *) viewer, PETSC_NULL);
1119: }
1120: #if 0
1121: nmax = 4;
1122: loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
1123: PetscOptionsGetIntArray(prefix, "-gvec_ts_errormonitor", loc, &nmax, &opt);
1124: if (opt) {
1125: PetscViewerDrawOpenX(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);
1126: PetscViewerDrawGetDraw(viewer, &draw);
1127: PetscDrawSetTitle(draw, "Error");
1128: PetscLogObjectParent(ts, (PetscObject) viewer);
1129: TSSetMonitor(ts, GTSErrorMonitor, (void *) viewer);
1130: }
1131: #endif
1132: /*--------------------------------------------------------------------- */
1133: PetscOptionsHasName(PETSC_NULL, "-help", &opt);
1134: if (opt) {
1135: char *pprefix;
1136: int len = 2;
1137: if (prefix != PETSC_NULL) {
1138: PetscStrlen(prefix, &len);
1139: }
1140: PetscMalloc((len+2) * sizeof(char), &pprefix);
1141: PetscStrcpy(pprefix, "-");
1142: if (prefix != PETSC_NULL)
1143: PetscStrcat(pprefix, prefix);
1144: PetscPrintf(comm," Additional TS Monitor options for grid vectorsn");
1145: PetscPrintf(comm," %sgvec_ts_solutionmonitorn", pprefix);
1146: PetscPrintf(comm," %sgvec_ts_errormonitor -- not working yetn", pprefix);
1147: PetscFree(pprefix);
1148: }
1149: return(0);
1150: }
1151: EXTERN_C_END
1153: /*---------------------------------------------- Context Functions --------------------------------------------------*/
1154: EXTERN_C_BEGIN
1155: #undef __FUNCT__
1157: int GTSContextDestroy_Private(GTSContext ctx)
1158: {
1160: if (--ctx->refct > 0)
1161: return(0);
1162: PetscLogObjectDestroy(ctx);
1163: PetscHeaderDestroy(ctx);
1164: return(0);
1165: }
1166: EXTERN_C_END
1168: #undef __FUNCT__
1170: /*@
1171: GTSCreateContext - This functions creates a GTSContext attached to the given context.
1173: Collective on GTS
1175: Input Parameters:
1176: . ts - The GTS
1177: . ltime - The time for the context
1178: . ctx - The parent context
1180: Output Parameter:
1181: . oldCtx - The previous GTSContext
1183: Level: developer
1185: .keywords: grid ts, context, create
1186: .seealso: GTSDestroyContext(), GTSCreateConstraintContext()
1187: @*/
1188: int GTSCreateContext(GTS ts, double ltime, PetscObject ctx, PetscObject *oldCtx)
1189: {
1190: GTSContext GTSctx;
1191: int ierr;
1194: /* Preserve previous context */
1195: PetscObjectQuery(ctx, "GTSContext", oldCtx);
1196: if (*oldCtx != PETSC_NULL) {
1197: PetscObjectReference(*oldCtx);
1198: }
1199: /* Setup new context */
1200: PetscHeaderCreate(GTSctx, _GTSContext, int, TS_COOKIE, -1, "context", ts->comm, GTSContextDestroy_Private, 0);
1201: PetscLogObjectCreate(GTSctx);
1202: GTSctx->t = ltime;
1203: GTSctx->dt = ts->time_step;
1204: GTSctx->sol = ts->vec_sol;
1205: GTSctx->nwork = ts->nwork;
1206: GTSctx->work = ts->work;
1207: PetscObjectCompose(ctx, "GTSContext", (PetscObject) GTSctx);
1208: /* This makes sure the context is destroyed when the object is destroyed */
1209: PetscObjectDereference((PetscObject) GTSctx);
1210: return(0);
1211: }
1213: #undef __FUNCT__
1215: /*@
1216: GTSDestroyContext - This functions destroys a GTSContext attached to the given context and
1217: reattaches the old context.
1219: Collective on GTS
1221: Input Parameters:
1222: . ts - The GTS
1223: . ctx - The parent context
1224: . oldCtx - The previous GTSContext
1226: Level: developer
1228: .keywords: grid ts, context, destroy
1229: .seealso: GTSCreateContext(), GTSCreateConstraintContext()
1230: @*/
1231: int GTSDestroyContext(GTS ts, PetscObject ctx, PetscObject oldCtx)
1232: {
1236: if (oldCtx != PETSC_NULL) {
1237: /* Restore previous context -- New one is destroyed automatically */
1238: PetscObjectCompose(ctx, "GTSContext", oldCtx);
1239: PetscObjectDereference(oldCtx);
1240: }
1241: return(0);
1242: }
1244: #undef __FUNCT__
1246: /*@
1247: GTSCreateConstraintContext - This functions creates a GTSContext attached to the
1248: constraint context of the underlying grid.
1250: Collective on GTS
1252: Input Parameter:
1253: . ts - The GTS
1255: Level: developer
1257: .keywords: grid ts, context
1258: .seealso: GTSCreateContext(), GTSDestroyContext()
1259: @*/
1260: int GTSCreateConstraintContext(GTS ts)
1261: {
1262: Grid grid;
1263: PetscTruth isConstrained;
1264: PetscConstraintObject constCtx;
1265: PetscObject oldCtx;
1266: int ierr;
1270: if (ts->isGTS == PETSC_FALSE)
1271: return(0);
1272: GTSGetGrid(ts, &grid);
1273: GridIsConstrained(grid, &isConstrained);
1274: if (isConstrained == PETSC_TRUE) {
1275: GridGetConstraintContext(grid, &constCtx);
1276: GTSCreateContext(ts, ts->ptime, (PetscObject) constCtx, &oldCtx);
1277: /* Remove the extra reference from GTSCreateContext() */
1278: if (oldCtx != PETSC_NULL) {
1279: PetscObjectDereference(oldCtx);
1280: }
1281: }
1282: return(0);
1283: }