Commits

BarryFSmith  committed b3d3934

added -ts_monitor_envelope and TSMonitorEnvelope() to track extremes of solution vector
added support to extchem.c to print maximums of species
updated sample run for extchem.c for displaying only "important species"

  • Participants
  • Parent commits 80666b6

Comments (0)

Files changed (4)

File include/petsc-private/tsimpl.h

   void           *transformctx;
 };
 
+struct _n_TSMonitorEnvelopeCtx {
+  Vec max,min;
+};
+
 #endif

File include/petscts.h

 PETSC_EXTERN PetscErrorCode TSMonitorLGTimeStep(TS,PetscInt,PetscReal,Vec,void *);
 PETSC_EXTERN PetscErrorCode TSMonitorLGSolution(TS,PetscInt,PetscReal,Vec,void *);
 PETSC_EXTERN PetscErrorCode TSMonitorLGSetVariableNames(TS,const char * const*);
+PETSC_EXTERN PetscErrorCode TSMonitorLGGetVariableNames(TS,const char *const **);
 PETSC_EXTERN PetscErrorCode TSMonitorLGSetDisplayVariables(TS,const char * const*);
 PETSC_EXTERN PetscErrorCode TSMonitorLGSetTransform(TS,PetscErrorCode (*)(void*,Vec,Vec),void*);
 PETSC_EXTERN PetscErrorCode TSMonitorLGError(TS,PetscInt,PetscReal,Vec,void *);
 PETSC_EXTERN PetscErrorCode TSMonitorLGSNESIterations(TS,PetscInt,PetscReal,Vec,void *);
 PETSC_EXTERN PetscErrorCode TSMonitorLGKSPIterations(TS,PetscInt,PetscReal,Vec,void *);
 
+typedef struct _n_TSMonitorEnvelopeCtx*  TSMonitorEnvelopeCtx;
+PETSC_EXTERN PetscErrorCode TSMonitorEnvelopeCtxCreate(TS,TSMonitorEnvelopeCtx*);
+PETSC_EXTERN PetscErrorCode TSMonitorEnvelope(TS,PetscInt,PetscReal,Vec,void*);
+PETSC_EXTERN PetscErrorCode TSMonitorEnvelopeGetBounds(TS,Vec*,Vec*);
+PETSC_EXTERN PetscErrorCode TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx*);
+
 typedef struct _n_TSMonitorSPEigCtx*  TSMonitorSPEigCtx;
 PETSC_EXTERN PetscErrorCode TSMonitorSPEigCtxCreate(MPI_Comm,const char[],const char[],int,int,int,int,PetscInt,TSMonitorSPEigCtx *);
 PETSC_EXTERN PetscErrorCode TSMonitorSPEigCtxDestroy(TSMonitorSPEigCtx*);

File src/ts/examples/tutorials/extchem.c

        cp $PETSC_DIR/externallibaries/tchem/data/periodictable.dat .
 
     Run with
-   ./extchem -Tini 1500 -ts_arkimex_fully_implicit -ts_max_snes_failures -1 -ts_adapt_monitor -ts_adapt_dt_max 1e-4 -ts_arkimex_type 4 -ts_monitor_lg_solution -ts_final_time .005 -draw_pause -2 -lg_indicate_data_points false -ts_lg_monitor_solution_variables CH4,O2,N2,AR,CO,O
+   ./extchem -Tini 1500 -ts_arkimex_fully_implicit -ts_max_snes_failures -1 -ts_adapt_monitor -ts_adapt_dt_max 1e-4 -ts_arkimex_type 4 -ts_monitor_lg_solution -ts_final_time .005 -draw_pause -2 -lg_indicate_data_points false -ts_lg_monitor_solution_variables H2,O2,H2O,CH4,CO,CO2,C2H2,N2  -ts_monitor_envelope
 
     The solution for component i = 0 is the tempature.
 
 };
 
 
+static PetscErrorCode FormMoleFraction(User,Vec,Vec);
 static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
 static PetscErrorCode FormRHSJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
 static PetscErrorCode FormInitialSolution(TS,Vec,void*);
   ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps\n",TSConvergedReasons[reason],ftime,steps);CHKERRQ(ierr);
 
+  {
+    Vec                max;
+    const char * const *names;
+    PetscInt           i;
+    const PetscReal    *bmax;
+
+    ierr = TSMonitorEnvelopeGetBounds(ts,&max,NULL);CHKERRQ(ierr);
+    if (max) {
+      ierr = TSMonitorLGGetVariableNames(ts,&names);CHKERRQ(ierr);
+      ierr = VecGetArrayRead(max,&bmax);CHKERRQ(ierr);
+      ierr = PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n");CHKERRQ(ierr);
+      for (i=1; i<user.Nspec; i++) {
+        if (bmax[i] > .01) {ierr = PetscPrintf(PETSC_COMM_SELF,"%s %g\n",names[i],bmax[i]);CHKERRQ(ierr);}
+      }
+      ierr = VecRestoreArrayRead(max,&bmax);CHKERRQ(ierr);
+    }
+  }
+
   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      Free work space.
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   for (i=0; i<sizeof(initial)/sizeof(initial[0]); i++) {
     int ispec = TC_getSpos(initial[i].name, strlen(initial[i].name));
     if (ispec < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Could not find species %s",initial[i].name);
-    ierr = PetscPrintf(PETSC_COMM_SELF,"Species %d: %s\n",i,initial[i].name);CHKERRQ(ierr);
+    ierr = PetscPrintf(PETSC_COMM_SELF,"Species %d: %s %g\n",i,initial[i].name,initial[i].massfrac);CHKERRQ(ierr);
     x[1+ispec] = initial[i].massfrac;
   }
   ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);

File src/ts/interface/ts.c

 .  -ts_monitor_draw_solution_phase - Monitor solution graphically with phase diagram
 .  -ts_monitor_draw_error - Monitor error graphically
 .  -ts_monitor_solution_binary <filename> - Save each solution to a binary file
--  -ts_monitor_solution_vtk <filename.vts> - Save each time step to a binary file, use filename-%%03D.vts
+.  -ts_monitor_solution_vtk <filename.vts> - Save each time step to a binary file, use filename-%%03D.vts
+-  -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
 
    Developer Note: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
 
     ierr = TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy);CHKERRQ(ierr);
   }
 
+  ierr = PetscOptionsName("-ts_monitor_envelope","Monitor maximum and minimum value of each component of the solution","TSMonitorEnvelope",&opt);CHKERRQ(ierr);
+  if (opt) {
+    TSMonitorEnvelopeCtx ctx;
+
+    ierr = TSMonitorEnvelopeCtxCreate(ts,&ctx);CHKERRQ(ierr);
+    ierr = TSMonitorSet(ts,TSMonitorEnvelope,ctx,(PetscErrorCode (*)(void**))TSMonitorEnvelopeCtxDestroy);CHKERRQ(ierr);
+  }
+
   /*
      This code is all wrong. One is creating objects inside the TSSetFromOptions() so if run with the options gui
      will bleed memory. Also one is using a PetscOptionsBegin() inside a PetscOptionsBegin()
 }
 #endif
 
-
-
 #undef __FUNCT__
 #define __FUNCT__ "TSMonitorLGSolution"
 /*@C
 .  ptime - current time
 -  lg - a line graph object
 
+   Options Database:
+.   -ts_lg_monitor_solution_variables
+
    Level: intermediate
 
     Notes: each process in a parallel run displays its component solutions in a separate window
 
    Input Parameters:
 +  ts - the TS context
-.  names - the names of the components, final string must be NULL
+-  names - the names of the components, final string must be NULL
 
    Level: intermediate
 
       TSMonitorLGCtx  ctx = ts->monitorcontext[i];
       ierr = PetscStrArrayDestroy(&ctx->names);CHKERRQ(ierr);
       ierr = PetscStrArrayallocpy(names,&ctx->names);CHKERRQ(ierr);
+      break;
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSMonitorLGGetVariableNames"
+/*@C
+   TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot
+
+   Collective on TS
+
+   Input Parameter:
+.  ts - the TS context
+
+   Output Parameter:
+.  names - the names of the components, final string must be NULL
+
+   Level: intermediate
+
+.keywords: TS,  vector, monitor, view
+
+.seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
+@*/
+PetscErrorCode  TSMonitorLGGetVariableNames(TS ts,const char *const **names)
+{
+  PetscInt       i;
+
+  PetscFunctionBegin;
+  *names = NULL;
+  for (i=0; i<ts->numbermonitors; i++) {
+    if (ts->monitor[i] == TSMonitorLGSolution) {
+      TSMonitorLGCtx  ctx = ts->monitorcontext[i];
+      *names = (const char *const *)ctx->names;
+      break;
     }
   }
   PetscFunctionReturn(0);
         }
         j++;
       }
+      break;
     }
   }
   PetscFunctionReturn(0);
   ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
+
+/* ------------------------------------------------------------------------*/
+#undef __FUNCT__
+#define __FUNCT__ "TSMonitorEnvelopeCtxCreate"
+/*@C
+   TSMonitorEnvelopeCtxCreate - Creates a context for use with TSMonitorEnvelope()
+
+   Collective on TS
+
+   Input Parameters:
+.  ts  - the ODE solver object
+
+   Output Parameter:
+.  ctx - the context
+
+   Level: intermediate
+
+.keywords: TS, monitor, line graph, residual, seealso
+
+.seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()
+
+@*/
+PetscErrorCode  TSMonitorEnvelopeCtxCreate(TS ts,TSMonitorEnvelopeCtx *ctx)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = PetscNew(struct _n_TSMonitorEnvelopeCtx,ctx);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSMonitorEnvelope"
+/*@C
+   TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution
+
+   Collective on TS
+
+   Input Parameters:
++  ts - the TS context
+.  step - current time-step
+.  ptime - current time
+-  ctx - the envelope context
+
+   Options Database:
+.  -ts_monitor_envelope
+
+   Level: intermediate
+
+   Notes: after a solve you can use TSMonitorEnvelopeGetBounds() to access the envelope
+
+.keywords: TS,  vector, monitor, view
+
+.seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorEnvelopeGetBounds()
+@*/
+PetscErrorCode  TSMonitorEnvelope(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
+{
+  PetscErrorCode       ierr;
+  TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dummy;
+
+  PetscFunctionBegin;
+  if (!ctx->max) {
+    ierr = VecDuplicate(u,&ctx->max);CHKERRQ(ierr);
+    ierr = VecDuplicate(u,&ctx->min);CHKERRQ(ierr);
+    ierr = VecCopy(u,ctx->max);CHKERRQ(ierr);
+    ierr = VecCopy(u,ctx->min);CHKERRQ(ierr);
+  } else {
+    ierr = VecPointwiseMax(ctx->max,u,ctx->max);CHKERRQ(ierr);
+    ierr = VecPointwiseMin(ctx->min,u,ctx->min);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "TSMonitorEnvelopeGetBounds"
+/*@C
+   TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution
+
+   Collective on TS
+
+   Input Parameter:
+.  ts - the TS context
+
+   Output Parameter:
++  max - the maximum values
+-  min - the minimum values
+
+   Level: intermediate
+
+.keywords: TS,  vector, monitor, view
+
+.seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
+@*/
+PetscErrorCode  TSMonitorEnvelopeGetBounds(TS ts,Vec *max,Vec *min)
+{
+  PetscInt i;
+
+  PetscFunctionBegin;
+  if (max) *max = NULL;
+  if (min) *min = NULL;
+  for (i=0; i<ts->numbermonitors; i++) {
+    if (ts->monitor[i] == TSMonitorEnvelope) {
+      TSMonitorEnvelopeCtx  ctx = ts->monitorcontext[i];
+      if (max) *max = ctx->max;
+      if (min) *min = ctx->min;
+      break;
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSMonitorEnvelopeCtxDestroy"
+/*@C
+   TSMonitorEnvelopeCtxDestroy - Destroys a context that was created  with TSMonitorEnvelopeCtxCreate().
+
+   Collective on TSMonitorEnvelopeCtx
+
+   Input Parameter:
+.  ctx - the monitor context
+
+   Level: intermediate
+
+.keywords: TS, monitor, line graph, destroy
+
+.seealso: TSMonitorLGCtxCreate(),  TSMonitorSet(), TSMonitorLGTimeStep();
+@*/
+PetscErrorCode  TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = VecDestroy(&(*ctx)->min);CHKERRQ(ierr);
+  ierr = VecDestroy(&(*ctx)->max);CHKERRQ(ierr);
+  ierr = PetscFree(*ctx);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}