PetscPostTimeStep.cpp File Reference

Post-time-step function. More...

#include <stdio.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mpivars_cpp.h>
#include <simulation_object.h>
#include <petscinterface.h>

Go to the source code of this file.

Macros
#define	__FUNCT__   "PetscPostTimeStep"

Functions
int	OutputSolution (void *, int, double)
PetscErrorCode	PetscPostTimeStep (TS ts)

Detailed Description

Post-time-step function.

Author

Debojyoti Ghosh

Definition in file PetscPostTimeStep.cpp.

Macro Definition Documentation

__FUNCT__

#define __FUNCT__   "PetscPostTimeStep"

Definition at line 18 of file PetscPostTimeStep.cpp.

Function Documentation

OutputSolution()

int OutputSolution	(	void *	s,
		int	nsims,
		double	a_time
	)

Write solutions to file

Write out the solution to file

Parameters

s	Array of simulation objects of type SimulationObject
nsims	Number of simulation objects
a_time	Current simulation time

Definition at line 27 of file OutputSolution.cpp.

{
  SimulationObject* simobj = (SimulationObject*) s;
  int ns;
  _DECLARE_IERR_;

  for (ns = 0; ns < nsims; ns++) {

    HyPar*        solver = &(simobj[ns].solver);
    MPIVariables* mpi    = &(simobj[ns].mpi);

    if ((!solver->WriteOutput) && (strcmp(solver->plot_solution,"yes"))) continue;

    /* time integration module may have auxiliary arrays to write out, so get them */
    int NSolutions = 0;
    IERR TimeGetAuxSolutions(&NSolutions,NULL,solver,-1,ns); CHECKERR(ierr);
    if (NSolutions > 10) NSolutions = 10;

    int  nu;
    char fname_root[_MAX_STRING_SIZE_];
    char aux_fname_root[_MAX_STRING_SIZE_];
    strcpy(fname_root, solver->op_fname_root);
    strcpy(aux_fname_root, solver->aux_op_fname_root);

    if (nsims > 1) {
      char index[_MAX_STRING_SIZE_];
      GetStringFromInteger(ns, index, (int)log10(nsims)+1);
      strcat(fname_root, "_");
      strcat(fname_root, index);
      strcat(aux_fname_root, "_");
      strcat(aux_fname_root, index);
    }

    for (nu=0; nu<NSolutions; nu++) {

      double  *uaux = NULL;
      IERR TimeGetAuxSolutions(&NSolutions,&uaux,solver,nu,ns); CHECKERR(ierr);

      IERR WriteArray(  solver->ndims,
                        solver->nvars,
                        solver->dim_global,
                        solver->dim_local,
                        solver->ghosts,
                        solver->x,
                        uaux,
                        solver,
                        mpi,
                        aux_fname_root ); CHECKERR(ierr);

      aux_fname_root[2]++;
    }

#if defined(HAVE_CUDA)
    if (solver->use_gpu) {
      /* Copy values from GPU memory to CPU memory for writing. */
      gpuMemcpy(solver->x, solver->gpu_x, sizeof(double)*solver->size_x, gpuMemcpyDeviceToHost);

#ifdef CUDA_VAR_ORDERDING_AOS
      gpuMemcpy(  solver->u, 
                  solver->gpu_u, 
                  sizeof(double)*solver->ndof_cells_wghosts, 
                  gpuMemcpyDeviceToHost );
#else
      double *h_u = (double *) malloc(sizeof(double)*solver->ndof_cells_wghosts);
      gpuMemcpy(h_u, solver->gpu_u, sizeof(double)*solver->ndof_cells_wghosts, gpuMemcpyDeviceToHost);
      for (int i=0; i<solver->npoints_local_wghosts; i++) {
        for (int v=0; v<solver->nvars; v++) {
          solver->u[i*solver->nvars+v] = h_u[i+v*solver->npoints_local_wghosts];
        }
      }
      free(h_u);
#endif
    }
#endif

    WriteArray(  solver->ndims,
                 solver->nvars,
                 solver->dim_global,
                 solver->dim_local,
                 solver->ghosts,
                 solver->x,
                 solver->u,
                 solver,
                 mpi,
                 fname_root );

    if (!strcmp(solver->plot_solution, "yes")) {
      PlotArray(   solver->ndims,
                   solver->nvars,
                   solver->dim_global,
                   solver->dim_local,
                   solver->ghosts,
                   solver->x,
                   solver->u,
                   a_time,
                   solver,
                   mpi,
                   fname_root );
    }

    /* increment the index string, if required */
    if ((!strcmp(solver->output_mode,"serial")) && (!strcmp(solver->op_overwrite,"no"))) {
      IncrementFilenameIndex(solver->filename_index,solver->index_length);
    }

  }

  return(0);
}

PetscPostTimeStep()

PetscErrorCode PetscPostTimeStep

(

TS

ts

)

Function called after every time step

Parameters

ts	Time integrator object

Definition at line 21 of file PetscPostTimeStep.cpp.

{
  PETScContext* context(nullptr);

  PetscFunctionBegin;

  TSGetApplicationContext(ts,&context);
  if (!context) {
    fprintf(stderr,"Error in PetscPreTimeStep: Null context!\n");
    return(1);
  }
  SimulationObject* sim = (SimulationObject*) context->simobj;
  int nsims = context->nsims;

  Vec Y;
  TSGetSolution(ts,&Y);

  double waqt;
  TSGetTime(ts,&waqt);
  int iter;
  TSGetStepNumber(ts,&iter);
  double dt;
  TSGetTimeStep(ts,&dt);

  context->tic++;

  double max_cfl = 0.0;
  double total_norm = 0.0;
  for (int ns = 0; ns < nsims; ns++) {

    HyPar* solver = &(sim[ns].solver);
    MPIVariables* mpi = &(sim[ns].mpi);

    /* get solution */
    TransferVecFromPETSc(solver->u,Y,context,ns,context->offsets[ns]);

    /* Call any physics-specific post-step function */
    if (solver->PostStep) solver->PostStep(solver->u,solver,mpi,waqt,iter);

    /* Calculate CFL and diffusion number */
    double local_max_cfl  = -1.0, global_max_cfl  = -1.0;
    if (solver->ComputeCFL) local_max_cfl = solver->ComputeCFL(solver,mpi,dt,waqt);
    MPIMax_double(&global_max_cfl ,&local_max_cfl ,1,&mpi->world);
    if (global_max_cfl > max_cfl) max_cfl = global_max_cfl;

    /* Calculate norm of the change in the solution for this time step */
    _ArrayAXPY_(solver->u,-1.0,solver->u0,(solver->npoints_local_wghosts*solver->nvars));
    double sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
                                  solver->ghosts,solver->index,solver->u0);
    double global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
    double norm = sqrt((global_sum/(double)solver->npoints_global));
    total_norm += norm;
  
    if (!strcmp(solver->ConservationCheck,"yes")) {
      /* calculate volume integral of the solution at this time step */
      solver->VolumeIntegralFunction(solver->VolumeIntegral,solver->u,solver,mpi);
      /* calculate surface integral of the flux at this time step */
      solver->BoundaryIntegralFunction(solver,mpi);
      /* calculate the conservation error at this time step       */
      solver->CalculateConservationError(solver,mpi);
    }

  }

  gettimeofday(&context->iter_end_time,NULL);
  long long walltime;
  walltime = (  (context->iter_end_time.tv_sec * 1000000 + context->iter_end_time.tv_usec)
              - (context->iter_start_time.tv_sec * 1000000 + context->iter_start_time.tv_usec));
  context->iter_wctime = (double) walltime / 1000000.0;
  context->ti_runtime += context->iter_wctime;

  if ((!context->rank) && (iter%sim[0].solver.screen_op_iter == 0)) {

    if (nsims > 1) {
      printf("--\n");
      printf("iter=%7d,  t=%1.3e\n", iter, waqt);
      printf("  dt=%1.3E  ", dt);
      printf("  CFL=%1.3E, ", max_cfl);
      printf("  norm=%1.4E\n", total_norm);
      printf("  wctime=%1.1E (s)\n",context->iter_wctime);
    } else {
      printf("iter=%7d  ", iter);
      printf("dt=%1.3E  ", dt);
      printf("t=%1.3E  ", waqt);
      printf("CFL=%1.3E  ", max_cfl );
      printf("norm=%1.4E  ", total_norm);
      printf("wctime: %1.1E (s)  ",context->iter_wctime);
    }

    /* calculate and print conservation error */
    if (!strcmp(sim[0].solver.ConservationCheck,"yes")) {

      double error = 0;
      for  (int ns = 0; ns < nsims; ns++) {
        for (int v=0; v<sim[ns].solver.nvars; v++) {
          error += (sim[ns].solver.ConservationError[v] * sim[ns].solver.ConservationError[v]);
        }
      }
      error = sqrt(error);
      printf("  cons_err=%1.4E\n", error);

    } else {

      if (nsims == 1) printf("\n");

    }
  
    /* print physics-specific info, if available */
    for( int ns = 0; ns < nsims; ns++) {
      if (sim[ns].solver.PrintStep) {
        if (nsims > 1) printf("Physics-specific output for domain %d:\n", ns);
        sim[ns].solver.PrintStep( &(sim[ns].solver),
                                  &(sim[ns].mpi),
                                  waqt );
      }
    }

    if (nsims > 1) {
      printf("--\n");
      printf("\n");
    }

  }

  /* Write intermediate solution to file */
  if (iter%sim[0].solver.file_op_iter == 0) { 
    for (int ns = 0; ns < nsims; ns++) {
      HyPar* solver = &(sim[ns].solver);
      MPIVariables* mpi = &(sim[ns].mpi);
      if (solver->PhysicsOutput) solver->PhysicsOutput(solver,mpi,waqt);
    }
    OutputSolution(sim,nsims,waqt);
#ifdef with_librom
    context->op_times_arr.push_back(waqt);
#endif
    context->tic=0;
  }

  PetscFunctionReturn(0);
}