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, &GTSctx);                                                  */
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, &gts);
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), &gts->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), &gts->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: }