simulation.h File Reference

Base class for simulation and declarations for C functions. More...

#include <simulation_object.h>
#include <mpivars_cpp.h>
#include <petscinterface.h>

Go to the source code of this file.

Data Structures
class	Simulation
	Base class for a simulation. More...

Functions
int	ReadInputs (void *, int, int)
int	WriteInputs (void *, int, int)
int	Initialize (void *, int)
int	Initialize_GPU (void *, int)
int	InitialSolution (void *, int)
int	InitializeBoundaries (void *, int)
int	InitializeImmersedBoundaries (void *, int)
int	InitializePhysics (void *, int)
int	InitializePhysicsData (void *, int, int, int *)
int	InitializeSolvers (void *, int)
int	Cleanup (void *, int)
void	SimWriteErrors (void *, int, int, double, double)
int	SolvePETSc (void *, int, int, int)
	Integrate in time with PETSc. More...
int	Solve (void *, int, int, int)

Detailed Description

Base class for simulation and declarations for C functions.

Author

Debojyoti Ghosh

Definition in file simulation.h.

Function Documentation

int ReadInputs	(	void *	s,
		int	nsims,
		int	rank
	)

Read the input parameters

Read the simulation inputs from the file solver.inp. Rank 0 reads in the inputs and broadcasts them to all the processors.

The format of solver.inp is as follows:

begin
    <keyword>   <value>
    <keyword>   <value>
    <keyword>   <value>
    ...
    <keyword>   <value>
end

where the list of keywords and their type are:

Keyword name	Type	Variable	Default value
ndims	int	HyPar::ndims	1
nvars	int	HyPar::nvars	1
size	int[ndims]	HyPar::dim_global	must be specified
iproc	int[ndims]	MPIVariables::iproc	must be specified (see notes below)
ghost	int	HyPar::ghosts	1
n_iter	int	HyPar::n_iter	0
restart_iter	int	HyPar::restart_iter	0
time_scheme	char[]	HyPar::time_scheme	euler
time_scheme_type	char[]	HyPar::time_scheme_type	none
hyp_space_scheme	char[]	HyPar::spatial_scheme_hyp	1
hyp_flux_split	char[]	HyPar::SplitHyperbolicFlux	no
hyp_interp_type	char[]	HyPar::interp_type	characteristic
par_space_type	char[]	HyPar::spatial_type_par	nonconservative-1stage
par_space_scheme	char[]	HyPar::spatial_scheme_par	2
dt	double	HyPar::dt	0.0
conservation_check	char[]	HyPar::ConservationCheck	no
screen_op_iter	int	HyPar::screen_op_iter	1
file_op_iter	int	HyPar::file_op_iter	1000
op_file_format	char[]	HyPar::op_file_format	text
ip_file_type	char[]	HyPar::ip_file_type	ascii
input_mode	char[]	HyPar::input_mode	serial
output_mode	char[]	HyPar::output_mode	serial
op_overwrite	char[]	HyPar::op_overwrite	no
plot_solution	char[]	HyPar::plot_solution	no
model	char[]	HyPar::model	must be specified
immersed_body	char[]	HyPar::ib_filename	"none"
size_exact	int[ndims]	HyPar::dim_global_ex	HyPar::dim_global
use_gpu	char[]	HyPar::use_gpu	no
gpu_device_no	int	HyPar::gpu_device_no	-1

Notes:

• "ndims" must be specified before "size".
• the input "iproc" is ignored when running a sparse grids simulation.

• if "input_mode" or "output_mode" are set to "parallel" or "mpi-io", the number of I/O ranks must be specified right after as an integer. For example:

  begin
      ...
      input_mode  parallel 4
      ...
  end
  

  This means that 4 MPI ranks will participate in file I/O (assuming total MPI ranks is more than 4) (see ReadArrayParallel(), WriteArrayParallel(), ReadArrayMPI_IO() ).

  • The number of I/O ranks specified for "input_mode" and "output_mode" must be same. Otherwise, the value for the one specified last will be used.
  • The number of I/O ranks must be such that the total number of MPI ranks is an integer multiple. Otherwise, the code will use only 1 I/O rank.
• If any of the keywords are not present, the default value is used, except the ones whose default values say "must be specified". Thus, keywords that are not required for a particular simulation may be left out of the solver.inp input file. For example,
  • a Euler1D simulation does not need "par_space_type" or "par_space_scheme" because it does not have a parabolic term.
  • unless a conservation check is required, "conservation_check" can be left out and the code will not check for conservation.
  • "immersed_body" need not be specified if there are no immersed bodies present. NOTE: However, if it is specified, and a file of that filename does not exist, it will result in an error.

Parameters

s	Array of simulation objects of type SimulationObject of size nsims
nsims	Number of simulation objects
rank	MPI rank of this process

Definition at line 93 of file ReadInputs.c.

{
  SimulationObject *sim = (SimulationObject*) s;
  int n, ferr    = 0;

  if (sim == NULL) {
    printf("Error: simulation object array is NULL!\n");
    printf("Please consider killing this run.\n");
    return(1);
  }

  if (!rank) {

    /* set some default values for optional inputs */
    for (n = 0; n < nsims; n++) {
      sim[n].solver.ndims           = 1;
      sim[n].solver.nvars           = 1;
      sim[n].solver.ghosts          = 1;
      sim[n].solver.dim_global      = NULL;
      sim[n].solver.dim_local       = NULL;
      sim[n].solver.dim_global_ex   = NULL;
      sim[n].mpi.iproc              = NULL;
      sim[n].mpi.N_IORanks          = 1;
      sim[n].solver.dt              = 0.0;
      sim[n].solver.n_iter          = 0;
      sim[n].solver.restart_iter    = 0;
      sim[n].solver.screen_op_iter  = 1;
      sim[n].solver.file_op_iter    = 1000;
      sim[n].solver.write_residual  = 0;
      sim[n].solver.flag_ib         = 0;
#if defined(HAVE_CUDA)
      sim[n].solver.use_gpu         = 0;
      sim[n].solver.gpu_device_no   = -1;
#endif
      strcpy(sim[n].solver.time_scheme        ,"euler"         );
      strcpy(sim[n].solver.time_scheme_type   ," "             );
      strcpy(sim[n].solver.spatial_scheme_hyp ,"1"             );
      strcpy(sim[n].solver.spatial_type_par   ,_NC_1STAGE_     );
      strcpy(sim[n].solver.spatial_scheme_par ,"2"             );
      strcpy(sim[n].solver.interp_type        ,"characteristic");
      strcpy(sim[n].solver.ip_file_type       ,"ascii"         );
      strcpy(sim[n].solver.input_mode         ,"serial"        );
      strcpy(sim[n].solver.output_mode        ,"serial"        );
      strcpy(sim[n].solver.op_file_format     ,"text"          );
      strcpy(sim[n].solver.op_overwrite       ,"no"            );
      strcpy(sim[n].solver.plot_solution      ,"no"            );
      strcpy(sim[n].solver.model              ,"none"          );
      strcpy(sim[n].solver.ConservationCheck  ,"no"            );
      strcpy(sim[n].solver.SplitHyperbolicFlux,"no"            );
      strcpy(sim[n].solver.ib_filename        ,"none"          );
    }

    /* open the file */
    FILE *in;
    printf("Reading solver inputs from file \"solver.inp\".\n");
    in = fopen("solver.inp","r");
    if (!in) {
      fprintf(stderr,"Error: File \"solver.inp\" not found.\n");
      fprintf(stderr,"Please consider killing this run.\n");
      return(1);
    }

    /* reading solver inputs */
    char word[_MAX_STRING_SIZE_];
    ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);

    if (!strcmp(word, "begin")){

      while (strcmp(word, "end")) {

          ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);

        if (!strcmp(word, "ndims")) {

          ferr = fscanf(in,"%d",&(sim[0].solver.ndims)); if (ferr != 1) return(1);
          sim[0].solver.dim_global    = (int*) calloc (sim[0].solver.ndims,sizeof(int));
          sim[0].mpi.iproc            = (int*) calloc (sim[0].solver.ndims,sizeof(int));
          sim[0].solver.dim_global_ex = (int*) calloc (sim[0].solver.ndims,sizeof(int));

          int n;
          for (n = 1; n < nsims; n++) {
            sim[n].solver.ndims = sim[0].solver.ndims;
            sim[n].solver.dim_global    = (int*) calloc (sim[n].solver.ndims,sizeof(int));
            sim[n].mpi.iproc            = (int*) calloc (sim[n].solver.ndims,sizeof(int));
            sim[n].solver.dim_global_ex = (int*) calloc (sim[n].solver.ndims,sizeof(int));
          }

        } else if (!strcmp(word, "nvars")) {

          ferr = fscanf(in,"%d",&(sim[0].solver.nvars));
          for (int n = 1; n < nsims; n++) sim[n].solver.nvars = sim[0].solver.nvars;

        } else if   (!strcmp(word, "size")) {

          for (int n = 0; n < nsims; n++) {
            if (!sim[n].solver.dim_global) {
              fprintf(stderr,"Error in ReadInputs(): dim_global not allocated for n=%d.\n", n);
              fprintf(stderr,"Please specify ndims before dimensions.\n"         );
              return(1);
            } else {
              for (int i=0; i<sim[n].solver.ndims; i++) {
                ferr = fscanf(in,"%d",&(sim[n].solver.dim_global[i]));
                if (ferr != 1) {
                  fprintf(stderr,"Error in ReadInputs() while reading grid sizes for domain %d.\n", n);
                  return(1);
                }
                sim[n].solver.dim_global_ex[i] = sim[n].solver.dim_global[i];
              }
            }
          }

        } else if   (!strcmp(word, "size_exact")) {

          for (int n = 0; n < nsims; n++) {
            if (!sim[n].solver.dim_global_ex) {
              fprintf(stderr,"Error in ReadInputs(): dim_global_ex not allocated for n=%d.\n", n);
              fprintf(stderr,"Please specify ndims before dimensions.\n"         );
              return(1);
            } else {
              for (int i=0; i<sim[n].solver.ndims; i++) {
                ferr = fscanf(in,"%d",&(sim[n].solver.dim_global_ex[i]));
                if (ferr != 1) {
                  fprintf(stderr,"Error in ReadInputs() while reading exact solution grid sizes for domain %d.\n", n);
                  return(1);
                }
              }
            }
          }

        } else if (!strcmp(word, "iproc")) {

          int n;
          for (n = 0; n < nsims; n++) {
            if (!sim[n].mpi.iproc) {
              fprintf(stderr,"Error in ReadInputs(): iproc not allocated for n=%d.\n", n);
              fprintf(stderr,"Please specify ndims before iproc.\n"         );
              return(1);
            } else {
              int i;
              for (i=0; i<sim[n].solver.ndims; i++) {
                ferr = fscanf(in,"%d",&(sim[n].mpi.iproc[i]));
                if (ferr != 1) {
                  fprintf(stderr,"Error in ReadInputs() while reading iproc for domain %d.\n", n);
                  return(1);
                }
              }
            }
          }

              } else if (!strcmp(word, "ghost")) {

          ferr = fscanf(in,"%d",&(sim[0].solver.ghosts));

          int n;
          for (n = 1; n < nsims; n++) sim[n].solver.ghosts = sim[0].solver.ghosts;

        } else if (!strcmp(word, "n_iter")) {

          ferr = fscanf(in,"%d",&(sim[0].solver.n_iter));

          int n;
          for (n = 1; n < nsims; n++) sim[n].solver.n_iter = sim[0].solver.n_iter;

        } else if (!strcmp(word, "restart_iter")) {

          ferr = fscanf(in,"%d",&(sim[0].solver.restart_iter));

          int n;
          for (n = 1; n < nsims; n++) sim[n].solver.restart_iter = sim[0].solver.restart_iter;

        } else if (!strcmp(word, "time_scheme")) {

          ferr = fscanf(in,"%s",sim[0].solver.time_scheme);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.time_scheme, sim[0].solver.time_scheme);

        }   else if (!strcmp(word, "time_scheme_type" )) {

          ferr = fscanf(in,"%s",sim[0].solver.time_scheme_type);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.time_scheme_type, sim[0].solver.time_scheme_type);

        }   else if (!strcmp(word, "hyp_space_scheme")) {

          ferr = fscanf(in,"%s",sim[0].solver.spatial_scheme_hyp);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.spatial_scheme_hyp, sim[0].solver.spatial_scheme_hyp);

        }   else if (!strcmp(word, "hyp_flux_split")) {

          ferr = fscanf(in,"%s",sim[0].solver.SplitHyperbolicFlux);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.SplitHyperbolicFlux, sim[0].solver.SplitHyperbolicFlux);

        }   else if (!strcmp(word, "hyp_interp_type")) {

          ferr = fscanf(in,"%s",sim[0].solver.interp_type);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.interp_type, sim[0].solver.interp_type);

        }   else if (!strcmp(word, "par_space_type")) {

          ferr = fscanf(in,"%s",sim[0].solver.spatial_type_par);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.spatial_type_par, sim[0].solver.spatial_type_par);

        }   else if (!strcmp(word, "par_space_scheme")) {

          ferr = fscanf(in,"%s",sim[0].solver.spatial_scheme_par);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.spatial_scheme_par, sim[0].solver.spatial_scheme_par);

        }   else if (!strcmp(word, "dt")) {

          ferr = fscanf(in,"%lf",&(sim[0].solver.dt));

          int n;
          for (n = 1; n < nsims; n++) sim[n].solver.dt = sim[0].solver.dt;

        }   else if (!strcmp(word, "conservation_check" )) {

          ferr = fscanf(in,"%s",sim[0].solver.ConservationCheck);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.ConservationCheck, sim[0].solver.ConservationCheck);

        }   else if (!strcmp(word, "screen_op_iter")) {

          ferr = fscanf(in,"%d",&(sim[0].solver.screen_op_iter));

          int n;
          for (n = 1; n < nsims; n++) sim[n].solver.screen_op_iter = sim[0].solver.screen_op_iter;

        }   else if (!strcmp(word, "file_op_iter")) {

          ferr = fscanf(in,"%d",&(sim[0].solver.file_op_iter));

          int n;
          for (n = 1; n < nsims; n++) sim[n].solver.file_op_iter = sim[0].solver.file_op_iter;

        }   else if (!strcmp(word, "op_file_format")) {

          ferr = fscanf(in,"%s",sim[0].solver.op_file_format);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.op_file_format, sim[0].solver.op_file_format);

        }   else if (!strcmp(word, "ip_file_type")) {

          ferr = fscanf(in,"%s",sim[0].solver.ip_file_type);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.ip_file_type, sim[0].solver.ip_file_type);

        }   else if (!strcmp(word, "input_mode")) {

          ferr = fscanf(in,"%s",sim[0].solver.input_mode);
          if (strcmp(sim[0].solver.input_mode,"serial")) ferr = fscanf(in,"%d",&(sim[0].mpi.N_IORanks));

          int n;
          for (n = 1; n < nsims; n++) {
            strcpy(sim[n].solver.input_mode, sim[0].solver.input_mode);
            if (strcmp(sim[n].solver.input_mode,"serial")) sim[n].mpi.N_IORanks = sim[0].mpi.N_IORanks;
          }

            } else if (!strcmp(word, "output_mode"))  {

          ferr = fscanf(in,"%s",sim[0].solver.output_mode);
          if (strcmp(sim[0].solver.output_mode,"serial")) ferr = fscanf(in,"%d",&(sim[0].mpi.N_IORanks));

          int n;
          for (n = 1; n < nsims; n++) {
            strcpy(sim[n].solver.output_mode, sim[0].solver.output_mode);
            if (strcmp(sim[n].solver.output_mode,"serial")) sim[n].mpi.N_IORanks = sim[0].mpi.N_IORanks;
          }

        } else if   (!strcmp(word, "op_overwrite")) {

          ferr = fscanf(in,"%s",sim[0].solver.op_overwrite);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.op_overwrite, sim[0].solver.op_overwrite);

        } else if   (!strcmp(word, "plot_solution")) {

          ferr = fscanf(in,"%s",sim[0].solver.plot_solution);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.plot_solution, sim[0].solver.plot_solution);

        }   else if (!strcmp(word, "model")) {

          ferr = fscanf(in,"%s",sim[0].solver.model);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.model, sim[0].solver.model);

        }   else if (!strcmp(word, "immersed_body")) {

          ferr = fscanf(in,"%s",sim[0].solver.ib_filename);

          int n;
          for (n = 1; n < nsims; n++) strcpy(sim[n].solver.ib_filename, sim[0].solver.ib_filename);

        }
#if defined(HAVE_CUDA)
        else if (!strcmp(word, "use_gpu")) {
          ferr = fscanf(in,"%s",word);
          if (!strcmp(word, "yes") || !strcmp(word, "true")) sim[0].solver.use_gpu = 1;

          int n;
          for (n = 1; n < nsims; n++) sim[n].solver.use_gpu = sim[0].solver.use_gpu;
        } else if (!strcmp(word, "gpu_device_no")) {
          ferr = fscanf(in,"%d", &sim[0].solver.gpu_device_no);

          int n;
          for (n = 1; n < nsims; n++) sim[n].solver.gpu_device_no = sim[0].solver.gpu_device_no;
        }
#endif
        else if (strcmp(word, "end")) {

          char useless[_MAX_STRING_SIZE_];
          ferr = fscanf(in,"%s",useless);
          printf("Warning: keyword %s in file \"solver.inp\" with value %s not recognized or extraneous. Ignoring.\n",
                  word,useless);

        }
        if (ferr != 1) return(1);

      }

    } else {

          fprintf(stderr,"Error: Illegal format in file \"solver.inp\".\n");
      return(1);

    }

    /* close the file */
    fclose(in);

    /* some checks */
    for (n = 0; n < nsims; n++) {

      if (sim[n].solver.screen_op_iter <= 0)  sim[n].solver.screen_op_iter = 1;
      if (sim[n].solver.file_op_iter <= 0)    sim[n].solver.file_op_iter   = sim[n].solver.n_iter;

      if ((sim[n].solver.ndims != 3) && (strcmp(sim[n].solver.ib_filename,"none"))) {
        printf("Warning: immersed boundaries not implemented for ndims = %d. ",sim[n].solver.ndims);
        printf("Ignoring input for \"immersed_body\" (%s).\n",sim[n].solver.ib_filename);
        strcpy(sim[n].solver.ib_filename,"none");
      }
      sim[n].solver.flag_ib = strcmp(sim[n].solver.ib_filename,"none");

      /* restart only supported for binary output files */
      if ((sim[n].solver.restart_iter != 0) && strcmp(sim[n].solver.op_file_format,"binary")) {
        if (!sim[n].mpi.rank) fprintf(stderr,"Error in ReadInputs(): Restart is supported only for binary output files.\n");
        return(1);
      }
    }
  }

#ifndef serial
  for (n = 0; n < nsims; n++) {

    /* Broadcast the input parameters */
    MPIBroadcast_integer(&(sim[n].solver.ndims),1,0,&(sim[n].mpi.world));
    if (sim[n].mpi.rank) {
      sim[n].solver.dim_global    = (int*) calloc (sim[n].solver.ndims,sizeof(int));
      sim[n].mpi.iproc            = (int*) calloc (sim[n].solver.ndims,sizeof(int));
      sim[n].solver.dim_global_ex = (int*) calloc (sim[n].solver.ndims,sizeof(int));
    }
    MPIBroadcast_integer(&(sim[n].solver.nvars)         ,1                  ,0,&(sim[n].mpi.world));
    MPIBroadcast_integer( sim[n].solver.dim_global      ,sim[n].solver.ndims,0,&(sim[n].mpi.world));
    MPIBroadcast_integer( sim[n].solver.dim_global_ex   ,sim[n].solver.ndims,0,&(sim[n].mpi.world));
    MPIBroadcast_integer( sim[n].mpi.iproc              ,sim[n].solver.ndims,0,&(sim[n].mpi.world));
    MPIBroadcast_integer(&(sim[n].mpi.N_IORanks)        ,1                  ,0,&(sim[n].mpi.world));
    MPIBroadcast_integer(&(sim[n].solver.ghosts)        ,1                  ,0,&(sim[n].mpi.world));
    MPIBroadcast_integer(&(sim[n].solver.n_iter)        ,1                  ,0,&(sim[n].mpi.world));
    MPIBroadcast_integer(&(sim[n].solver.restart_iter)  ,1                  ,0,&(sim[n].mpi.world));
    MPIBroadcast_integer(&(sim[n].solver.screen_op_iter),1                  ,0,&(sim[n].mpi.world));
    MPIBroadcast_integer(&(sim[n].solver.file_op_iter)  ,1                  ,0,&(sim[n].mpi.world));
    MPIBroadcast_integer(&(sim[n].solver.flag_ib)       ,1                  ,0,&(sim[n].mpi.world));
#if defined(HAVE_CUDA)
    MPIBroadcast_integer(&(sim[n].solver.use_gpu)       ,1                  ,0,&(sim[n].mpi.world));
    MPIBroadcast_integer(&(sim[n].solver.gpu_device_no) ,1                  ,0,&(sim[n].mpi.world));
#endif
    MPIBroadcast_character(sim[n].solver.time_scheme        ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.time_scheme_type   ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.spatial_scheme_hyp ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.interp_type        ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.spatial_type_par   ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.spatial_scheme_par ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.ConservationCheck  ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.SplitHyperbolicFlux,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.op_file_format     ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.ip_file_type       ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.input_mode         ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.output_mode        ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.op_overwrite       ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.plot_solution      ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.model              ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));
    MPIBroadcast_character(sim[n].solver.ib_filename        ,_MAX_STRING_SIZE_,0,&(sim[n].mpi.world));

    MPIBroadcast_double(&(sim[n].solver.dt),1,0,&(sim[n].mpi.world));
  }
#endif

  return 0;
}

int WriteInputs	(	void *	s,
		int	nsims,
		int	rank
	)

Write the input parameters

Write the simulation inputs read from the file solver.inp.

Parameters

s	Array of simulation objects of type SimulationObject of size nsims
nsims	Number of simulation objects
rank	MPI rank of this process

Definition at line 15 of file WriteInputs.c.

{
  SimulationObject *sim = (SimulationObject*) s;
  int n;

  if (sim == NULL)  return 0;

  if (!rank) {

    printf("  No. of dimensions                          : %d\n",sim[0].solver.ndims);
    printf("  No. of variables                           : %d\n",sim[0].solver.nvars);
    if (nsims > 1) {
      printf("  Domain sizes:\n");
      for (int n = 0; n < nsims; n++) {
        printf("    domain %3d - ", n);
        for (int i=0; i<sim[n].solver.ndims; i++) printf ("%d ",sim[n].solver.dim_global[i]);
        printf("\n");
      }
#ifndef serial
        printf("  Processes along each dimension:\n");
      for (int n = 0; n < nsims; n++) {
        printf("    domain %3d - ", n);
        for (int i=0; i<sim[n].solver.ndims; i++) printf ("%d ",sim[n].mpi.iproc[i]);
        printf("\n");
      }
#endif
      printf("  Exact solution domain sizes:\n");
      for (int n = 0; n < nsims; n++) {
        printf("    domain %3d - ", n);
        for (int i=0; i<sim[n].solver.ndims; i++) printf ("%d ",sim[n].solver.dim_global_ex[i]);
        printf("\n");
      }
    } else {
        printf("  Domain size                                : ");
      for (int i=0; i<sim[0].solver.ndims; i++) printf ("%d ",sim[0].solver.dim_global[i]);
      printf("\n");
#ifndef serial
        printf("  Processes along each dimension             : ");
      for (int i=0; i<sim[0].solver.ndims; i++) printf ("%d ",sim[0].mpi.iproc[i]);
      printf("\n");
#endif
        printf("  Exact solution domain size                 : ");
      for (int i=0; i<sim[0].solver.ndims; i++) printf ("%d ",sim[0].solver.dim_global_ex[i]);
      printf("\n");
    }
      printf("  No. of ghosts pts                          : %d\n"     ,sim[0].solver.ghosts              );
      printf("  No. of iter.                               : %d\n"     ,sim[0].solver.n_iter              );
      printf("  Restart iteration                          : %d\n"     ,sim[0].solver.restart_iter        );
#ifdef with_petsc
    if (sim[0].solver.use_petscTS)
      printf("  Time integration scheme                    : PETSc \n"                            );
    else {
      printf("  Time integration scheme                    : %s ",sim[0].solver.time_scheme             );
      if (strcmp(sim[0].solver.time_scheme,_FORWARD_EULER_)) {
        printf("(%s)",sim[0].solver.time_scheme_type                                                    );
      }
      printf("\n");
    }
#else
    printf("  Time integration scheme                    : %s ",sim[0].solver.time_scheme               );
    if (strcmp(sim[0].solver.time_scheme,_FORWARD_EULER_)) {
      printf("(%s)",sim[0].solver.time_scheme_type                                                      );
    }
    printf("\n");
#endif
    printf("  Spatial discretization scheme (hyperbolic) : %s\n"     ,sim[0].solver.spatial_scheme_hyp  );
    printf("  Split hyperbolic flux term?                : %s\n"     ,sim[0].solver.SplitHyperbolicFlux );
    printf("  Interpolation type for hyperbolic term     : %s\n"     ,sim[0].solver.interp_type         );
    printf("  Spatial discretization type   (parabolic ) : %s\n"     ,sim[0].solver.spatial_type_par    );
    printf("  Spatial discretization scheme (parabolic ) : %s\n"     ,sim[0].solver.spatial_scheme_par  );
    printf("  Time Step                                  : %E\n"     ,sim[0].solver.dt                  );
    printf("  Check for conservation                     : %s\n"     ,sim[0].solver.ConservationCheck   );
    printf("  Screen output iterations                   : %d\n"     ,sim[0].solver.screen_op_iter      );
    printf("  File output iterations                     : %d\n"     ,sim[0].solver.file_op_iter        );
    printf("  Initial solution file type                 : %s\n"     ,sim[0].solver.ip_file_type        );
    printf("  Initial solution read mode                 : %s"       ,sim[0].solver.input_mode          );
    if (strcmp(sim[0].solver.input_mode,"serial"))    printf("  [%d file IO rank(s)]\n",sim[0].mpi.N_IORanks  );
    else                                        printf("\n");
    printf("  Solution file write mode                   : %s"       ,sim[0].solver.output_mode         );
    if (strcmp(sim[0].solver.output_mode,"serial"))   printf("  [%d file IO rank(s)]\n",sim[0].mpi.N_IORanks  );
    else                                        printf("\n");
    printf("  Solution file format                       : %s\n"     ,sim[0].solver.op_file_format      );
    printf("  Overwrite solution file                    : %s\n"     ,sim[0].solver.op_overwrite        );
#if defined(HAVE_CUDA)
    printf("  Use GPU                                    : %s\n"     ,(sim[0].solver.use_gpu == 1)? "yes" : "no");
    printf("  GPU device no                              : %d\n"     ,(sim[0].solver.gpu_device_no));
#endif
    printf("  Physical model                             : %s\n"     ,sim[0].solver.model               );
    if (sim[0].solver.flag_ib) {
      printf("  Immersed Body                              : %s\n"     ,sim[0].solver.ib_filename         );
    }
  }

  return 0;
}

int Initialize	(	void *	s,
		int	nsims
	)

Initialize the solver

Initialization function called at the beginning of a simulation. This function does the following:

• allocates memory for MPI related arrays
• initializes the values for MPI variables
• creates sub-communicators and communication groups
• allocates memory for arrays to store solution, right-hand-side, flux, and other working vectors.
• initializes function counters to zero

Parameters

s	Array of simulation objects of type SimulationObject
nsims	Number of simulation objects

Definition at line 26 of file Initialize.c.

{
  SimulationObject* simobj = (SimulationObject*) s;
  int i,d,n;

  if (nsims == 0) {
    return 1;
  }

#if defined(HAVE_CUDA)
  if (simobj[0].solver.use_gpu && (simobj[0].solver.gpu_device_no >= 0)) {
      gpuSetDevice(simobj[0].solver.gpu_device_no);
  }
#endif

  if (!simobj[0].mpi.rank)  printf("Partitioning domain and allocating data arrays.\n");

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

    /* this is a full initialization, not a barebones one */
    simobj[n].is_barebones = 0;

    /* allocations */
    simobj[n].mpi.ip           = (int*) calloc (simobj[n].solver.ndims,sizeof(int));
    simobj[n].mpi.is           = (int*) calloc (simobj[n].solver.ndims,sizeof(int));
    simobj[n].mpi.ie           = (int*) calloc (simobj[n].solver.ndims,sizeof(int));
    simobj[n].mpi.bcperiodic   = (int*) calloc (simobj[n].solver.ndims,sizeof(int));
    simobj[n].solver.dim_local = (int*) calloc (simobj[n].solver.ndims,sizeof(int));
    simobj[n].solver.isPeriodic= (int*) calloc (simobj[n].solver.ndims,sizeof(int));

#if defined(HAVE_CUDA)
    simobj[n].mpi.wctime = 0;
    simobj[n].mpi.wctime_total = 0;
#endif

#ifndef serial
    _DECLARE_IERR_;

    /* Domain partitioning */
    int total_proc = 1;
    for (i=0; i<simobj[n].solver.ndims; i++) total_proc *= simobj[n].mpi.iproc[i];
    if (simobj[n].mpi.nproc != total_proc) {
      fprintf(stderr,"Error on rank %d: total number of processes is not consistent ", simobj[n].mpi.rank);
      fprintf(stderr,"with number of processes along each dimension.\n");
      if (nsims > 1) fprintf(stderr,"for domain %d.\n", n);
      fprintf(stderr,"mpiexec was called with %d processes, ",simobj[n].mpi.nproc);
      fprintf(stderr,"total number of processes from \"solver.inp\" is %d.\n", total_proc);
      return(1);
    }

    /* calculate ndims-D rank of each process (ip[]) from rank in MPI_COMM_WORLD */
    IERR MPIRanknD( simobj[n].solver.ndims,
                    simobj[n].mpi.rank,
                    simobj[n].mpi.iproc,
                    simobj[n].mpi.ip); CHECKERR(ierr);

    /* calculate local domain sizes along each dimension */
    for (i=0; i<simobj[n].solver.ndims; i++) {
      simobj[n].solver.dim_local[i] = MPIPartition1D( simobj[n].solver.dim_global[i],
                                                      simobj[n].mpi.iproc[i],
                                                      simobj[n].mpi.ip[i] );
    }

    /* calculate local domain limits in terms of global domain */
    IERR MPILocalDomainLimits(  simobj[n].solver.ndims,
                                simobj[n].mpi.rank,
                                &(simobj[n].mpi),
                                simobj[n].solver.dim_global,
                                simobj[n].mpi.is,
                                simobj[n].mpi.ie  );
    CHECKERR(ierr);

    /* create sub-communicators for parallel computations along grid lines in each dimension */
    IERR MPICreateCommunicators(simobj[n].solver.ndims,&(simobj[n].mpi)); CHECKERR(ierr);

    /* initialize periodic BC flags to zero */
    for (i=0; i<simobj[n].solver.ndims; i++) simobj[n].mpi.bcperiodic[i] = 0;

    /* create communication groups */
    IERR MPICreateIOGroups(&(simobj[n].mpi)); CHECKERR(ierr);

#else

    for (i=0; i<simobj[n].solver.ndims; i++) {
      simobj[n].mpi.ip[i]            = 0;
      simobj[n].solver.dim_local[i]  = simobj[n].solver.dim_global[i];
      simobj[n].mpi.iproc[i]         = 1;
      simobj[n].mpi.is[i]            = 0;
      simobj[n].mpi.ie[i]            = simobj[n].solver.dim_local[i];
      simobj[n].mpi.bcperiodic[i]    = 0;
    }

#endif

    simobj[n].solver.npoints_global
      = simobj[n].solver.npoints_local
      = simobj[n].solver.npoints_local_wghosts
      = 1;
    for (i=0; i<simobj[n].solver.ndims; i++) {
      simobj[n].solver.npoints_global *= simobj[n].solver.dim_global[i];
      simobj[n].solver.npoints_local *= simobj[n].solver.dim_local [i];
      simobj[n].solver.npoints_local_wghosts *= (simobj[n].solver.dim_local[i]+2*simobj[n].solver.ghosts);
    }

    /* Allocations */
    simobj[n].solver.index = (int*) calloc ((short)simobj[n].solver.ndims,sizeof(int));
    simobj[n].solver.stride_with_ghosts = (int*) calloc ((short)simobj[n].solver.ndims,sizeof(int));
    simobj[n].solver.stride_without_ghosts = (int*) calloc ((short)simobj[n].solver.ndims,sizeof(int));
    int accu1 = 1, accu2 = 1;
    for (i=0; i<simobj[n].solver.ndims; i++) {
      simobj[n].solver.stride_with_ghosts[i]    = accu1;
      simobj[n].solver.stride_without_ghosts[i] = accu2;
      accu1 *= (simobj[n].solver.dim_local[i]+2*simobj[n].solver.ghosts);
      accu2 *=  simobj[n].solver.dim_local[i];
    }

#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      gpuMalloc((void**)&simobj[n].solver.gpu_dim_local, simobj[n].solver.ndims*sizeof(int));
      gpuMemcpy(  simobj[n].solver.gpu_dim_local, 
                  simobj[n].solver.dim_local, 
                  simobj[n].solver.ndims*sizeof(int), 
                  gpuMemcpyHostToDevice );
    }
#endif

    /* state variables */
    int size = 1;
    for (i=0; i<simobj[n].solver.ndims; i++) {
      size *= (simobj[n].solver.dim_local[i]+2*simobj[n].solver.ghosts);
    }
    simobj[n].solver.ndof_cells_wghosts = simobj[n].solver.nvars*size;
    simobj[n].solver.u = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      gpuMalloc((void**)&simobj[n].solver.gpu_u, simobj[n].solver.nvars*size*sizeof(double));
      gpuMemset(simobj[n].solver.gpu_u, 0, simobj[n].solver.nvars*size*sizeof(double));
    }
#endif
#ifdef with_petsc
    if (simobj[n].solver.use_petscTS) {
      simobj[n].solver.u0      = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
      simobj[n].solver.uref    = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
      simobj[n].solver.rhsref  = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
      simobj[n].solver.rhs     = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
    } else simobj[n].solver.u0 = simobj[n].solver.uref = simobj[n].solver.rhsref = simobj[n].solver.rhs = NULL;
#endif
#ifdef with_librom
    simobj[n].solver.u_rom_predicted = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
#endif
    simobj[n].solver.hyp     = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
    simobj[n].solver.par     = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
    simobj[n].solver.source  = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
    simobj[n].solver.iblank  = (double*) calloc (size              ,sizeof(double));

#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      gpuMalloc((void**)&simobj[n].solver.hyp, simobj[n].solver.nvars*size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.par, simobj[n].solver.nvars*size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.source, simobj[n].solver.nvars*size*sizeof(double));
      gpuMemset(simobj[n].solver.hyp, 0, simobj[n].solver.nvars*size*sizeof(double));
      gpuMemset(simobj[n].solver.par, 0, simobj[n].solver.nvars*size*sizeof(double));
      gpuMemset(simobj[n].solver.source, 0, simobj[n].solver.nvars*size*sizeof(double));
    } else {
#endif
      simobj[n].solver.hyp     = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
      simobj[n].solver.par     = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
      simobj[n].solver.source  = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
#if defined(HAVE_CUDA)
    }
#endif

    simobj[n].solver.iblank  = (double*) calloc (size,sizeof(double));
#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      gpuMalloc((void**)&simobj[n].solver.gpu_iblank, size*sizeof(double));
      gpuMemset(simobj[n].solver.gpu_iblank, 0, size*sizeof(double));
    }
#endif
    
    /* grid */
    size = 0;
    for (i=0; i<simobj[n].solver.ndims; i++) {
      size += (simobj[n].solver.dim_local[i]+2*simobj[n].solver.ghosts);
    }
    simobj[n].solver.x     = (double*) calloc (size,sizeof(double));
    simobj[n].solver.dxinv = (double*) calloc (size,sizeof(double));
    simobj[n].solver.size_x = size;
#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      gpuMalloc((void**)&simobj[n].solver.gpu_x, size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.gpu_dxinv, size*sizeof(double));
      gpuMemset(simobj[n].solver.gpu_x, 0, size*sizeof(double));
      gpuMemset(simobj[n].solver.gpu_dxinv, 0, size*sizeof(double));
    }
#endif
    
    /* cell-centered arrays needed to compute fluxes */
    size = 1;  
    for (i=0; i<simobj[n].solver.ndims; i++) {
      size *= (simobj[n].solver.dim_local[i]+2*simobj[n].solver.ghosts);
    }
#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      gpuMalloc((void**)&simobj[n].solver.fluxC, simobj[n].solver.nvars*size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.uC, simobj[n].solver.nvars*size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.Deriv1, simobj[n].solver.nvars*size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.Deriv2, simobj[n].solver.nvars*size*sizeof(double));
      gpuMemset(simobj[n].solver.fluxC, 0, simobj[n].solver.nvars*size*sizeof(double));
      gpuMemset(simobj[n].solver.uC, 0, simobj[n].solver.nvars*size*sizeof(double));
      gpuMemset(simobj[n].solver.Deriv1, 0, simobj[n].solver.nvars*size*sizeof(double));
      gpuMemset(simobj[n].solver.Deriv2, 0, simobj[n].solver.nvars*size*sizeof(double));
    } else {
#endif
      simobj[n].solver.uC     = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
      simobj[n].solver.fluxC  = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
      simobj[n].solver.Deriv1 = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
      simobj[n].solver.Deriv2 = (double*) calloc (simobj[n].solver.nvars*size,sizeof(double));
#if defined(HAVE_CUDA)
    }
#endif

    /* node-centered arrays needed to compute fluxes */
    size = 1;  for (i=0; i<simobj[n].solver.ndims; i++) size *= (simobj[n].solver.dim_local[i]+1);
    size *= simobj[n].solver.nvars;
    simobj[n].solver.ndof_nodes = size;
#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      gpuMalloc((void**)&simobj[n].solver.fluxI, size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.uL, size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.uR, size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.fL, size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.fR, size*sizeof(double));
      gpuMemset(simobj[n].solver.fluxI, 0, size*sizeof(double));
      gpuMemset(simobj[n].solver.uL, 0, size*sizeof(double));
      gpuMemset(simobj[n].solver.uR, 0, size*sizeof(double));
      gpuMemset(simobj[n].solver.fL, 0, size*sizeof(double));
      gpuMemset(simobj[n].solver.fR, 0, size*sizeof(double));
    } else {
#endif
      simobj[n].solver.fluxI = (double*) calloc (size,sizeof(double));
      simobj[n].solver.uL    = (double*) calloc (size,sizeof(double));
      simobj[n].solver.uR    = (double*) calloc (size,sizeof(double));
      simobj[n].solver.fL    = (double*) calloc (size,sizeof(double));
      simobj[n].solver.fR    = (double*) calloc (size,sizeof(double));
#if defined(HAVE_CUDA)
    }
#endif
    
    /* allocate MPI send/receive buffer arrays */
    int bufdim[simobj[n].solver.ndims], maxbuf = 0;
    for (d = 0; d < simobj[n].solver.ndims; d++) {
      bufdim[d] = 1;
      for (i = 0; i < simobj[n].solver.ndims; i++) {
        if (i == d) bufdim[d] *= simobj[n].solver.ghosts;
        else        bufdim[d] *= simobj[n].solver.dim_local[i];
      }
      if (bufdim[d] > maxbuf) maxbuf = bufdim[d];
    }
    maxbuf *= (simobj[n].solver.nvars*simobj[n].solver.ndims);
    simobj[n].mpi.maxbuf  = maxbuf;
    simobj[n].mpi.sendbuf = (double*) calloc (2*simobj[n].solver.ndims*maxbuf,sizeof(double));
    simobj[n].mpi.recvbuf = (double*) calloc (2*simobj[n].solver.ndims*maxbuf,sizeof(double));
#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      simobj[n].mpi.cpu_dim = (int *) calloc(simobj[n].solver.ndims, sizeof(int));
      _ArrayCopy1D_(simobj[n].solver.dim_local, simobj[n].mpi.cpu_dim, simobj[n].solver.ndims);
      gpuMalloc((void**)&simobj[n].mpi.gpu_sendbuf, 2*simobj[n].solver.ndims*simobj[n].mpi.maxbuf*sizeof(double));
      gpuMalloc((void**)&simobj[n].mpi.gpu_recvbuf, 2*simobj[n].solver.ndims*simobj[n].mpi.maxbuf*sizeof(double));
      gpuMemset(simobj[n].mpi.gpu_sendbuf, 0, 2*simobj[n].solver.ndims*simobj[n].mpi.maxbuf*sizeof(double));
      gpuMemset(simobj[n].mpi.gpu_recvbuf, 0, 2*simobj[n].solver.ndims*simobj[n].mpi.maxbuf*sizeof(double));
    }
#endif
    
    /* allocate the volume and boundary integral arrays */
    simobj[n].solver.VolumeIntegral        = (double*) calloc (simobj[n].solver.nvars  ,sizeof(double));
    simobj[n].solver.VolumeIntegralInitial = (double*) calloc (simobj[n].solver.nvars  ,sizeof(double));
    simobj[n].solver.TotalBoundaryIntegral = (double*) calloc (simobj[n].solver.nvars,sizeof(double));
    simobj[n].solver.ConservationError     = (double*) calloc (simobj[n].solver.nvars,sizeof(double));
    for (i=0; i<simobj[n].solver.nvars; i++) simobj[n].solver.ConservationError[i] = -1;
#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      int total_offset = 0;
      for (d=0; d<simobj[n].solver.ndims; d++) {
          simobj[n].solver.gpu_npoints_boundary_offset[d] = total_offset;
          simobj[n].solver.gpu_npoints_boundary[d] = 1;

          for (i=0; i<simobj[n].solver.ndims; i++) {
              if (i != d) simobj[n].solver.gpu_npoints_boundary[d] *= simobj[n].solver.dim_local[i];
          }
          total_offset += 2*simobj[n].solver.gpu_npoints_boundary[d];
      }
      simobj[n].solver.StageBoundaryBuffer_size = (total_offset*simobj[n].solver.nvars);
      gpuMalloc((void**)&simobj[n].solver.StageBoundaryBuffer, simobj[n].solver.StageBoundaryBuffer_size*sizeof(double));
      gpuMemset(simobj[n].solver.StageBoundaryBuffer, 0, simobj[n].solver.StageBoundaryBuffer_size*sizeof(double));

      size = 2*simobj[n].solver.ndims*simobj[n].solver.nvars;
      gpuMalloc((void**)&simobj[n].solver.StageBoundaryIntegral, size*sizeof(double));
      gpuMalloc((void**)&simobj[n].solver.StepBoundaryIntegral, size*sizeof(double));
      gpuMemset(simobj[n].solver.StageBoundaryIntegral, 0, size*sizeof(double));
      gpuMemset(simobj[n].solver.StepBoundaryIntegral, 0, size*sizeof(double));
    } else {
#endif
      simobj[n].solver.StageBoundaryIntegral = (double*) calloc (2*simobj[n].solver.ndims*simobj[n].solver.nvars,sizeof(double));
      simobj[n].solver.StepBoundaryIntegral  = (double*) calloc (2*simobj[n].solver.ndims*simobj[n].solver.nvars,sizeof(double));
#if defined(HAVE_CUDA)
    }
#endif

    /* initialize function call counts to zero */
    simobj[n].solver.count_hyp
      = simobj[n].solver.count_par
      = simobj[n].solver.count_sou
      = 0;
#ifdef with_petsc
    simobj[n].solver.count_RHSFunction
      = simobj[n].solver.count_IFunction
      = simobj[n].solver.count_IJacobian
      = simobj[n].solver.count_IJacFunction
      = 0;
#endif

    /* Initialize iblank to 1*/
    _ArraySetValue_(simobj[n].solver.iblank,simobj[n].solver.npoints_local_wghosts,1);
#if defined(HAVE_CUDA)
    if (simobj[n].solver.use_gpu) {
      gpuArraySetValue(simobj[n].solver.gpu_iblank, simobj[n].solver.npoints_local_wghosts, 1.0);
    }
#endif

  }

  return 0;
}

int Initialize_GPU	(	void *	,
		int
	)

Initialize the solver

int InitialSolution	(	void *	s,
		int	nsims
	)

Read the initial solution

Read in initial solution from file, and compute grid spacing and volume integral of the initial solution

Parameters

s	Array of simulation objects of type SimulationObject
nsims	Number of simulation objects

Definition at line 25 of file InitialSolution.c.

{
  SimulationObject* simobj = (SimulationObject*) s;
  int n, flag, d, i, offset, ierr;

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

    int ghosts = simobj[n].solver.ghosts;

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

    ierr = ReadArray( simobj[n].solver.ndims,
                      simobj[n].solver.nvars,
                      simobj[n].solver.dim_global,
                      simobj[n].solver.dim_local,
                      simobj[n].solver.ghosts,
                      &(simobj[n].solver),
                      &(simobj[n].mpi),
                      simobj[n].solver.x,
                      simobj[n].solver.u,
                      fname_root,
                      &flag );
    if (ierr) {
      fprintf(stderr, "Error in InitialSolution() on rank %d.\n",
              simobj[n].mpi.rank);
      return ierr;
    }
    if (!flag) {
      fprintf(stderr,"Error: initial solution file not found.\n");
      return(1);
    }
    CHECKERR(ierr);

    /* exchange MPI-boundary values of u between processors */
    MPIExchangeBoundariesnD(  simobj[n].solver.ndims,
                              simobj[n].solver.nvars,
                              simobj[n].solver.dim_local,
                              simobj[n].solver.ghosts,
                              &(simobj[n].mpi),
                              simobj[n].solver.u  );

    /* calculate dxinv */
    offset = 0;
    for (d = 0; d < simobj[n].solver.ndims; d++) {
      for (i = 0; i < simobj[n].solver.dim_local[d]; i++) {
        simobj[n].solver.dxinv[i+offset+ghosts]
          = 2.0 / (simobj[n].solver.x[i+1+offset+ghosts]-simobj[n].solver.x[i-1+offset+ghosts]);
      }
      offset += (simobj[n].solver.dim_local[d] + 2*ghosts);
    }

    /* exchange MPI-boundary values of dxinv between processors */
    offset = 0;
    for (d = 0; d < simobj[n].solver.ndims; d++) {
      ierr = MPIExchangeBoundaries1D( &(simobj[n].mpi),
                                      &(simobj[n].solver.dxinv[offset]),
                                      simobj[n].solver.dim_local[d],
                                      ghosts,
                                      d,
                                      simobj[n].solver.ndims ); CHECKERR(ierr);
      if (ierr) {
        fprintf(stderr, "Error in InitialSolution() on rank %d.\n",
                simobj[n].mpi.rank);
        return ierr;
      }
      offset += (simobj[n].solver.dim_local[d] + 2*ghosts);
    }

    /* fill in ghost values of dxinv at physical boundaries by extrapolation */
    offset = 0;
    for (d = 0; d < simobj[n].solver.ndims; d++) {
      double *dxinv = &(simobj[n].solver.dxinv[offset]);
      int    *dim = simobj[n].solver.dim_local;
      if (simobj[n].mpi.ip[d] == 0) {
        /* fill left boundary along this dimension */
        for (i = 0; i < ghosts; i++) dxinv[i] = dxinv[ghosts];
      }
      if (simobj[n].mpi.ip[d] == simobj[n].mpi.iproc[d]-1) {
        /* fill right boundary along this dimension */
        for (i = dim[d]+ghosts; i < dim[d]+2*ghosts; i++) dxinv[i] = dxinv[dim[d]+ghosts-1];
      }
      offset  += (dim[d] + 2*ghosts);
    }

    /* calculate volume integral of the initial solution */
    ierr = VolumeIntegral(  simobj[n].solver.VolumeIntegralInitial,
                            simobj[n].solver.u,
                            &(simobj[n].solver),
                            &(simobj[n].mpi) ); CHECKERR(ierr);
    if (ierr) {
      fprintf(stderr, "Error in InitialSolution() on rank %d.\n",
              simobj[n].mpi.rank);
      return ierr;
    }
    if (!simobj[n].mpi.rank) {
      if (nsims > 1) printf("Volume integral of the initial solution on domain %d:\n", n);
      else           printf("Volume integral of the initial solution:\n");
      for (d=0; d<simobj[n].solver.nvars; d++) {
        printf("%2d:  %1.16E\n",d,simobj[n].solver.VolumeIntegralInitial[d]);
      }
    }
    /* Set initial total boundary flux integral to zero */
    _ArraySetValue_(simobj[n].solver.TotalBoundaryIntegral,simobj[n].solver.nvars,0);

  }

#if defined(HAVE_CUDA)
  if (simobj[0].solver.use_gpu) {
    for (int n = 0; n < nsims; n++) {
      int npoints_local_wghosts = simobj[n].solver.npoints_local_wghosts;
      int nvars                 = simobj[n].solver.nvars;
      int size_x                = simobj[n].solver.size_x;

      gpuMemcpy(simobj[n].solver.gpu_x,
                simobj[n].solver.x, simobj[n].solver.size_x*sizeof(double),
                gpuMemcpyHostToDevice);
      gpuMemcpy(simobj[n].solver.gpu_dxinv, simobj[n].solver.dxinv,
                simobj[n].solver.size_x*sizeof(double),
                gpuMemcpyHostToDevice);

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

  return 0;
}

int InitializeBoundaries	(	void *	s,
		int	nsims
	)

Initialize the boundary conditions

This function initializes the variables and functions related to implementing the boundary conditions.

• Rank 0 reads in the boundary conditions file and broadcasts the information to all processors.
• Depending on the type of boundary, additional information is read in. For example, for Dirichlet boundary, the Dirichlet value is read in.
• Allocate and initialize arrays and variables related to implementing the boundary conditions.
• Each rank finds out if the subdomain it owns abuts any of the boundaries specified.

Note that boundary conditions are implemented as boundary objects of the type DomainBoundary.

Parameters

s	Array of simulation objects of type SimulationObject
nsims	number of simulation objects

Definition at line 36 of file InitializeBoundaries.c.

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

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

    DomainBoundary  *boundary = NULL;
    HyPar           *solver   = &(sim[ns].solver);
    MPIVariables    *mpi      = &(sim[ns].mpi);
    int nb, ferr;

    /* root process reads boundary condition file */
    if (!mpi->rank) {

      char filename[_MAX_STRING_SIZE_] = "boundary";
      char filename_backup[_MAX_STRING_SIZE_] = "boundary";
      if (nsims > 1) {
        char index[_MAX_STRING_SIZE_];
        GetStringFromInteger(ns, index, (int)log10(nsims)+1);
        strcat(filename, "_");
        strcat(filename, index);
      }
      strcat(filename, ".inp");
      strcat(filename_backup, ".inp");

      FILE *in;
      in = fopen(filename,"r");
      if (!in) {
        in = fopen(filename_backup, "r");
        if (!in) {
          fprintf(stderr,"Error: boundary condition file %s or %s not found.\n",
                  filename, filename_backup );
          return(1);
        } else {
          if (nsims > 1) printf("Domain %d: ", ns);
          printf("Reading boundary conditions from %s.\n", filename_backup);
        }
      } else {
        if (nsims > 1) printf("Domain %d: ", ns);
        printf("Reading boundary conditions from %s.\n", filename);
      }

      /* read number of boundary conditions and allocate */
      ferr = fscanf(in,"%d",&solver->nBoundaryZones); if (ferr != 1) return(1);
      boundary = (DomainBoundary*) calloc (solver->nBoundaryZones,sizeof(DomainBoundary));
      for (nb = 0; nb < solver->nBoundaryZones; nb++) {
        boundary[nb].DirichletValue = boundary[nb].SpongeValue
                                   = boundary[nb].FlowVelocity
                                   = boundary[nb].UnsteadyDirichletData
                                   = NULL;
        boundary[nb].UnsteadyDirichletSize = NULL;
      }

      /* read each boundary condition */
      for (nb = 0; nb < solver->nBoundaryZones; nb++) {
        int d, v;
        boundary[nb].xmin = (double*) calloc (solver->ndims,sizeof(double)); /* deallocated in BCCleanup.c */
        boundary[nb].xmax = (double*) calloc (solver->ndims,sizeof(double)); /* deallocated in BCCleanup.c */

        ferr = fscanf(in,"%s",boundary[nb].bctype); if (ferr != 1) return(1);
        ferr = fscanf(in,"%d",&boundary[nb].dim  ); if (ferr != 1) return(1);
        ferr = fscanf(in,"%d",&boundary[nb].face ); if (ferr != 1) return(1);
        for (d=0; d < solver->ndims; d++) {
          ferr = fscanf(in,"%lf %lf", &boundary[nb].xmin[d], &boundary[nb].xmax[d]);
          if (ferr != 2) return(1);
        }

        /* read in boundary type-specific additional data if required */

        if (!strcmp(boundary[nb].bctype,_DIRICHLET_)) {
          boundary[nb].DirichletValue = (double*) calloc (solver->nvars,sizeof(double));
                                       /* deallocated in BCCleanup.c */
          /* read the Dirichlet value for each variable on this boundary */
          for (v = 0; v < solver->nvars; v++) ferr = fscanf(in,"%lf",&boundary[nb].DirichletValue[v]);
        }

        if (!strcmp(boundary[nb].bctype,_SPONGE_)) {
          boundary[nb].SpongeValue = (double*) calloc (solver->nvars,sizeof(double));
                                       /* deallocated in BCCleanup.c */
          /* read the sponge value for each variable on this boundary */
          for (v = 0; v < solver->nvars; v++) ferr = fscanf(in,"%lf",&boundary[nb].SpongeValue[v]);
        }

        if (    (!strcmp(boundary[nb].bctype,_SLIP_WALL_))
            ||  (!strcmp(boundary[nb].bctype,_NOSLIP_WALL_)) ) {
          boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
                                       /* deallocated in BCCleanup.c */
          /* read the wall velocity */
          for (v = 0; v < solver->ndims; v++) ferr = fscanf(in,"%lf",&boundary[nb].FlowVelocity[v]);
        }

        if (    (!strcmp(boundary[nb].bctype,_SW_SLIP_WALL_))
            ||  (!strcmp(boundary[nb].bctype,_SW_NOSLIP_WALL_)) ) {
          boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
                                       /* deallocated in BCCleanup.c */
          /* read the wall velocity */
          for (v = 0; v < solver->ndims; v++) ferr = fscanf(in,"%lf",&boundary[nb].FlowVelocity[v]);
        }

        if (!strcmp(boundary[nb].bctype,_SUBSONIC_INFLOW_)) {
          boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
                                       /* deallocated in BCCleanup.c */
          /* read in the inflow density and velocity */
          ferr = fscanf(in,"%lf",&boundary[nb].FlowDensity);
          for (v = 0; v < solver->ndims; v++) ferr = fscanf(in,"%lf",&boundary[nb].FlowVelocity[v]);
        }

        if (!strcmp(boundary[nb].bctype,_SUBSONIC_OUTFLOW_)) {
          /* read in the outflow pressure */
          ferr = fscanf(in,"%lf",&boundary[nb].FlowPressure);
        }

        if (!strcmp(boundary[nb].bctype,_SUBSONIC_AMBIVALENT_)) {
          boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
                                       /* deallocated in BCCleanup.c */
          /* read in the inflow density, velocity, and pressure */
          ferr = fscanf(in,"%lf",&boundary[nb].FlowDensity);
          for (v = 0; v < solver->ndims; v++) ferr = fscanf(in,"%lf",&boundary[nb].FlowVelocity[v]);
          ferr = fscanf(in,"%lf",&boundary[nb].FlowPressure);
        }

        if (!strcmp(boundary[nb].bctype,_SUPERSONIC_INFLOW_)) {
          boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
                                       /* deallocated in BCCleanup.c */
          /* read in the inflow density, velocity and pressure */
          ferr = fscanf(in,"%lf",&boundary[nb].FlowDensity);
          for (v = 0; v < solver->ndims; v++) ferr = fscanf(in,"%lf",&boundary[nb].FlowVelocity[v]);
          ferr = fscanf(in,"%lf",&boundary[nb].FlowPressure);
        }

        if (!strcmp(boundary[nb].bctype,_TURBULENT_SUPERSONIC_INFLOW_)) {
          boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
                                       /* deallocated in BCCleanup.c */
          /* read in the inflow density, velocity and pressure */
          ferr = fscanf(in,"%lf",&boundary[nb].FlowDensity);
          for (v = 0; v < solver->ndims; v++) ferr = fscanf(in,"%lf",&boundary[nb].FlowVelocity[v]);
          ferr = fscanf(in,"%lf",&boundary[nb].FlowPressure);
          ferr = fscanf(in,"%s" , boundary[nb].UnsteadyDirichletFilename);
        }

        if (    (!strcmp(boundary[nb].bctype,_THERMAL_SLIP_WALL_))
            ||  (!strcmp(boundary[nb].bctype,_THERMAL_NOSLIP_WALL_)) ){
          boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
                                       /* deallocated in BCCleanup.c */
          /* read the wall velocity */
          for (v = 0; v < solver->ndims; v++) ferr = fscanf(in,"%lf",&boundary[nb].FlowVelocity[v]);
          /* read in the filename where temperature data is available */
          ferr = fscanf(in,"%s" , boundary[nb].UnsteadyTemperatureFilename);
        }

        /* if boundary is periodic, let the MPI and HyPar know */
        if (!strcmp(boundary[nb].bctype,_PERIODIC_)) {
          solver->isPeriodic[boundary[nb].dim] = 1;
        }
        /*
          The MPI function to exchange internal (MPI) boundary information will handle
          periodic boundaries ONLY IF number of process along that dimension is more
          than 1.
        */
        if ((!strcmp(boundary[nb].bctype,_PERIODIC_)) && (mpi->iproc[boundary[nb].dim] > 1)) {
          mpi->bcperiodic[boundary[nb].dim] = 1;
        }

        /* some checks */
        if (boundary[nb].dim >= solver->ndims) {
          fprintf(stderr,"Error in reading boundary condition %d: dim %d is invalid (ndims = %d).\n",
                  nb,boundary[nb].dim,solver->ndims);
          return(1);
        }
        printf("  Boundary %30s:  Along dimension %2d and face %+1d\n",
                  boundary[nb].bctype,boundary[nb].dim,boundary[nb].face);
      }

      fclose(in);
      printf("%d boundary condition(s) read.\n",solver->nBoundaryZones);
    }

    /* tell other processes how many BCs are there and let them allocate */
    IERR MPIBroadcast_integer(&solver->nBoundaryZones,1,0,&mpi->world); CHECKERR(ierr);
    if (mpi->rank) {
      boundary = (DomainBoundary*) calloc (solver->nBoundaryZones,sizeof(DomainBoundary));
      for (nb = 0; nb < solver->nBoundaryZones; nb++) {
        boundary[nb].xmin = (double*) calloc (solver->ndims,sizeof(double)); /* deallocated in BCCleanup.c */
        boundary[nb].xmax = (double*) calloc (solver->ndims,sizeof(double)); /* deallocated in BCCleanup.c */
        boundary[nb].DirichletValue = boundary[nb].SpongeValue
                                   = boundary[nb].FlowVelocity
                                   = boundary[nb].UnsteadyDirichletData
                                   = NULL;
        boundary[nb].UnsteadyDirichletSize = NULL;
      }
    }

    /* communicate BC data to other processes */
    for (nb = 0; nb < solver->nBoundaryZones; nb++) {
      IERR MPIBroadcast_character(boundary[nb].bctype,_MAX_STRING_SIZE_,0,&mpi->world); CHECKERR(ierr);
      IERR MPIBroadcast_integer  (&boundary[nb].dim  ,1                ,0,&mpi->world); CHECKERR(ierr);
      IERR MPIBroadcast_integer  (&boundary[nb].face ,1                ,0,&mpi->world); CHECKERR(ierr);
      IERR MPIBroadcast_double   (boundary[nb].xmin  ,solver->ndims    ,0,&mpi->world); CHECKERR(ierr);
      IERR MPIBroadcast_double   (boundary[nb].xmax  ,solver->ndims    ,0,&mpi->world); CHECKERR(ierr);
    }
    IERR MPIBroadcast_integer(solver->isPeriodic,solver->ndims,0,&mpi->world);CHECKERR(ierr);

    /* broadcast periodic boundary info for MPI to all processes */
    IERR MPIBroadcast_integer(mpi->bcperiodic,solver->ndims,0,&mpi->world);CHECKERR(ierr);

    /* On other processes, if necessary, allocate and receive boundary-type-specific data */
    for (nb = 0; nb < solver->nBoundaryZones; nb++) {
      if (!strcmp(boundary[nb].bctype,_DIRICHLET_)) {
        if (mpi->rank)  boundary[nb].DirichletValue = (double*) calloc (solver->nvars,sizeof(double));
        IERR MPIBroadcast_double(boundary[nb].DirichletValue,solver->nvars,0,&mpi->world); CHECKERR(ierr);
      }

      if (!strcmp(boundary[nb].bctype,_SPONGE_)) {
        if (mpi->rank)  boundary[nb].SpongeValue = (double*) calloc (solver->nvars,sizeof(double));
        IERR MPIBroadcast_double(boundary[nb].SpongeValue,solver->nvars,0,&mpi->world); CHECKERR(ierr);
      }

      if (    (!strcmp(boundary[nb].bctype,_SLIP_WALL_))
          ||  (!strcmp(boundary[nb].bctype,_NOSLIP_WALL_)) ) {
        if (mpi->rank) boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
        IERR MPIBroadcast_double(boundary[nb].FlowVelocity,solver->ndims,0,&mpi->world); CHECKERR(ierr);
      }

      if (    (!strcmp(boundary[nb].bctype,_SW_SLIP_WALL_))
          ||  (!strcmp(boundary[nb].bctype,_SW_NOSLIP_WALL_)) ) {
        if (mpi->rank) boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
        IERR MPIBroadcast_double(boundary[nb].FlowVelocity,solver->ndims,0,&mpi->world); CHECKERR(ierr);
      }

      if (!strcmp(boundary[nb].bctype,_SUBSONIC_INFLOW_)) {
        if (mpi->rank) boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
        IERR MPIBroadcast_double(&boundary[nb].FlowDensity,1            ,0,&mpi->world); CHECKERR(ierr);
        IERR MPIBroadcast_double(boundary[nb].FlowVelocity,solver->ndims,0,&mpi->world); CHECKERR(ierr);
      }

      if (!strcmp(boundary[nb].bctype,_SUBSONIC_OUTFLOW_)) {
        IERR MPIBroadcast_double(&boundary[nb].FlowPressure,1,0,&mpi->world); CHECKERR(ierr);
      }

      if (!strcmp(boundary[nb].bctype,_SUBSONIC_AMBIVALENT_)) {
        if (mpi->rank) boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
        IERR MPIBroadcast_double(&boundary[nb].FlowDensity,1            ,0,&mpi->world); CHECKERR(ierr);
        IERR MPIBroadcast_double(boundary[nb].FlowVelocity,solver->ndims,0,&mpi->world); CHECKERR(ierr);
        IERR MPIBroadcast_double(&boundary[nb].FlowPressure,1,0,&mpi->world); CHECKERR(ierr);
      }

      if (!strcmp(boundary[nb].bctype,_SUPERSONIC_INFLOW_)) {
        if (mpi->rank) boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
        IERR MPIBroadcast_double(&boundary[nb].FlowDensity ,1            ,0,&mpi->world); CHECKERR(ierr);
        IERR MPIBroadcast_double(boundary[nb].FlowVelocity ,solver->ndims,0,&mpi->world); CHECKERR(ierr);
        IERR MPIBroadcast_double(&boundary[nb].FlowPressure,1            ,0,&mpi->world); CHECKERR(ierr);
      }

      if (!strcmp(boundary[nb].bctype,_TURBULENT_SUPERSONIC_INFLOW_)) {
        if (mpi->rank) boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
        IERR MPIBroadcast_double(&boundary[nb].FlowDensity ,1            ,0,&mpi->world); CHECKERR(ierr);
        IERR MPIBroadcast_double(boundary[nb].FlowVelocity ,solver->ndims,0,&mpi->world); CHECKERR(ierr);
        IERR MPIBroadcast_double(&boundary[nb].FlowPressure,1            ,0,&mpi->world); CHECKERR(ierr);
        /* allocate arrays and read in unsteady boundary data */
        IERR BCReadTurbulentInflowData(&boundary[nb],mpi,solver->ndims,solver->nvars,solver->dim_local); CHECKERR(ierr);
      }

      if (    (!strcmp(boundary[nb].bctype,_THERMAL_SLIP_WALL_))
          ||  (!strcmp(boundary[nb].bctype,_THERMAL_NOSLIP_WALL_)) ) {
        if (mpi->rank) boundary[nb].FlowVelocity = (double*) calloc (solver->ndims,sizeof(double));
        IERR MPIBroadcast_double(boundary[nb].FlowVelocity,solver->ndims,0,&mpi->world); CHECKERR(ierr);
        /* allocate arrays and read in boundary temperature data */
        IERR BCReadTemperatureData(&boundary[nb],mpi,solver->ndims,solver->nvars,solver->dim_local); CHECKERR(ierr);
      }

    }

    solver->boundary = boundary;

    /* each process calculates its own part of these boundaries */
    IERR CalculateLocalExtent(solver,mpi); CHECKERR(ierr);

#if defined(HAVE_CUDA)
    int bounds[GPU_MAX_NDIMS];
    if (sim[0].solver.use_gpu) {
      for (nb = 0; nb < solver->nBoundaryZones; nb++) {
        _ArraySubtract1D_(bounds,boundary[nb].ie,boundary[nb].is,solver->ndims);
  
        _ArrayProduct1D_(bounds,solver->ndims,boundary[nb].gpu_npoints_bounds);
        boundary[nb].gpu_npoints_local_wghosts = solver->npoints_local_wghosts;
  
        _ArrayProduct1D_(bounds,solver->ndims,boundary[nb].gpu_npoints_bounds);
        boundary[nb].gpu_npoints_local_wghosts = solver->npoints_local_wghosts;
  
        gpuMalloc((void**)&boundary[nb].gpu_is, solver->ndims*sizeof(int));
        gpuMalloc((void**)&boundary[nb].gpu_ie, solver->ndims*sizeof(int));
        gpuMalloc((void**)&boundary[nb].gpu_bounds, solver->ndims*sizeof(int));
        gpuMemcpy(boundary[nb].gpu_is, boundary[nb].is, solver->ndims*sizeof(int), gpuMemcpyHostToDevice);
        gpuMemcpy(boundary[nb].gpu_ie, boundary[nb].ie, solver->ndims*sizeof(int), gpuMemcpyHostToDevice);
        gpuMemcpy(boundary[nb].gpu_bounds, bounds, solver->ndims*sizeof(int), gpuMemcpyHostToDevice);
        if (   (!strcmp(boundary[nb].bctype,_SLIP_WALL_))
            || (!strcmp(boundary[nb].bctype,_NOSLIP_WALL_)) ) {
            gpuMalloc((void**)&boundary[nb].gpu_FlowVelocity, solver->ndims*sizeof(double));
            gpuMemcpy(  boundary[nb].gpu_FlowVelocity, 
                        boundary[nb].FlowVelocity, 
                        solver->ndims*sizeof(double), 
                        gpuMemcpyHostToDevice);
        }
      }
    }
#endif

    /* initialize function pointers for each boundary condition */
    for (nb = 0; nb < solver->nBoundaryZones; nb++) {
#if defined(HAVE_CUDA)
      BCInitialize(&boundary[nb], solver->use_gpu);
#else
      BCInitialize(&boundary[nb], 0);
#endif
    }

  }

  return 0;
}

int InitializeImmersedBoundaries	(	void *	s,
		int	nsims
	)

Initialize the immersed boundary conditions

Initialize the immersed boundaries, if present.

• Read in immersed body from STL file.
• Allocate and set up ImmersedBoundary object.
• Identify blanked-out grid points based on immersed body geometry.
• Identify and make a list of immersed boundary points on each rank.
• For each immersed boundary point, find the "nearest" facet.

Parameters

s	Array of simulation objects of type SimulationObject
nsims	Number of simulation objects

Definition at line 22 of file InitializeImmersedBoundaries.c.

{
  SimulationObject* simobj = (SimulationObject*) s;
  int n;

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

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

    ImmersedBoundary *ib       = NULL;
    Body3D           *body     = NULL;

    int stat, d, ndims = solver->ndims;

    if ((!solver->flag_ib) || (ndims != _IB_NDIMS_)) {
      solver->ib = NULL;
      continue;
    }

    /* Read in immersed body from file */
    IBReadBodySTL(&body,solver->ib_filename,mpi,&stat);
    if (stat) {
      if (!mpi->rank) {
        fprintf(stderr,"Error in InitializeImmersedBoundaries(): Unable to ");
        fprintf(stderr,"read immersed body from file %s.\n",solver->ib_filename);
      }
      solver->flag_ib = 0;
      solver->ib = NULL;
      return(1);
    }
    IBComputeBoundingBox(body);

    /* allocate immersed boundary object and set it up */
    ib = (ImmersedBoundary*) calloc (1, sizeof(ImmersedBoundary));
    ib->tolerance = 1e-12;
    ib->delta     = 1e-6;
    ib->itr_max   = 500;
    ib->body      = body;
    solver->ib    = ib;

    int     offset_global, offset_local,
            *dim_local  = solver->dim_local,
            *dim_global = solver->dim_global,
            ghosts      = solver->ghosts,
            size        = dim_global[0] + dim_global[1] + dim_global[2],
            count       = 0;
    double  *Xg         = (double*) calloc(size,sizeof(double));

    /* assemble the global grid on rank 0 */
    offset_global = offset_local = 0;
    for (d=0; d<ndims; d++) {
      IERR MPIGatherArray1D(mpi,(mpi->rank?NULL:&Xg[offset_global]),
                            &solver->x[offset_local+ghosts],
                            mpi->is[d],mpi->ie[d],dim_local[d],0); CHECKERR(ierr);
      offset_global += dim_global[d];
      offset_local  += dim_local [d] + 2*ghosts;
    }
    /* send the global grid to other ranks */
    MPIBroadcast_double(Xg,size,0,&mpi->world);
  
    /* identify whether this is a 3D or "pseudo-2D" simulation */
    IBIdentifyMode(Xg,dim_global,solver->ib);
  
    /* identify grid points inside the immersed body */
    int count_inside_body = 0;
    count = IBIdentifyBody(solver->ib,dim_global,dim_local,ghosts,mpi,Xg,solver->iblank);
    MPISum_integer(&count_inside_body,&count,1,&mpi->world);
    free(Xg);

    /* At ghost points corresponding to the physical boundary, extrapolate from the interior 
       (this should also work for bodies that are adjacent to physical boundaries). At interior
       (MPI) boundaries, exchange iblank across MPI ranks.
    */
    int indexb[ndims], indexi[ndims], bounds[ndims], offset[ndims];
    for (d = 0; d < ndims; d++) {
      /* left boundary */
      if (!mpi->ip[d]) {
        _ArrayCopy1D_(dim_local,bounds,ndims); bounds[d] = ghosts;
        _ArraySetValue_(offset,ndims,0); offset[d] = -ghosts;
        int done = 0; _ArraySetValue_(indexb,ndims,0);
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,ndims); indexi[d] = ghosts-1-indexb[d];
          int p1; _ArrayIndex1DWO_(ndims,dim_local,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (ndims,dim_local,indexi,ghosts,p2);
          solver->iblank[p1] = solver->iblank[p2];
          _ArrayIncrementIndex_(ndims,bounds,indexb,done);
        }
      }
      /* right boundary */
      if (mpi->ip[d] == mpi->iproc[d]-1) {
        _ArrayCopy1D_(dim_local,bounds,ndims); bounds[d] = ghosts;
        _ArraySetValue_(offset,ndims,0); offset[d] = dim_local[d];
        int done = 0; _ArraySetValue_(indexb,ndims,0);
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,ndims); indexi[d] = dim_local[d]-1-indexb[d];
          int p1; _ArrayIndex1DWO_(ndims,dim_local,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (ndims,dim_local,indexi,ghosts,p2);
          solver->iblank[p1] = solver->iblank[p2];
          _ArrayIncrementIndex_(ndims,bounds,indexb,done);
        }
      }
    }
    MPIExchangeBoundariesnD(ndims,1,dim_local,ghosts,mpi,solver->iblank);
  
    /* identify and create a list of immersed boundary points on each rank */
    int count_boundary_points = 0;
    count = IBIdentifyBoundary(solver->ib,mpi,dim_local,ghosts,solver->iblank);
    MPISum_integer(&count_boundary_points,&count,1,&mpi->world);
  
    /* find the nearest facet for each immersed boundary point */
    double ld = 0, xmin, xmax, ymin, ymax, zmin, zmax;
    _GetCoordinate_(0,0             ,dim_local,ghosts,solver->x,xmin);
    _GetCoordinate_(0,dim_local[0]-1,dim_local,ghosts,solver->x,xmax);
    _GetCoordinate_(1,0             ,dim_local,ghosts,solver->x,ymin);
    _GetCoordinate_(1,dim_local[1]-1,dim_local,ghosts,solver->x,ymax);
    _GetCoordinate_(2,0             ,dim_local,ghosts,solver->x,zmin);
    _GetCoordinate_(2,dim_local[2]-1,dim_local,ghosts,solver->x,zmax);
    double xlen = xmax - xmin;
    double ylen = ymax - ymin;
    double zlen = zmax - zmin;
    ld = max3(xlen,ylen,zlen);
    count = IBNearestFacetNormal(solver->ib,mpi,solver->x,ld,dim_local,ghosts);
    if (count) {
      fprintf(stderr, "Error in InitializeImmersedBoundaries():\n");
      fprintf(stderr, "  IBNearestFacetNormal() returned with error code %d on rank %d.\n",
              count, mpi->rank);
      return(count);
    }

    /* For the immersed boundary points, find the interior points for extrapolation,
       and compute their interpolation coefficients */
    count = IBInterpCoeffs(solver->ib,mpi,solver->x,dim_local,ghosts,solver->iblank);
    if (count) {
      fprintf(stderr, "Error in InitializeImmersedBoundaries():\n");
      fprintf(stderr, "  IBInterpCoeffs() returned with error code %d on rank %d.\n",
              count, mpi->rank);
      return(count);
    }
  
    /* Create facet mapping */;
    count = IBCreateFacetMapping(ib,mpi,solver->x,dim_local,ghosts);
    if (count) {
      fprintf(stderr, "Error in InitializeImmersedBoundaries():\n");
      fprintf(stderr, "  IBCreateFacetMapping() returned with error code %d on rank %d.\n",
              count, mpi->rank);
      return(count);
    }
    
    /* Done */
    if (!mpi->rank) {
      double percentage;
      printf("Immersed body read from %s:\n",solver->ib_filename);
      if (nsims > 1) printf("For domain %d,\n", n);
      printf("    Number of facets: %d\n    Bounding box: [%3.1f,%3.1lf] X [%3.1f,%3.1lf] X [%3.1f,%3.1lf]\n",
             body->nfacets,body->xmin,body->xmax,body->ymin,body->ymax,body->zmin,body->zmax);
      percentage = ((double)count_inside_body)/((double)solver->npoints_global)*100.0;
      printf("    Number of grid points inside immersed body: %d (%4.1f%%).\n",count_inside_body,percentage);
      percentage = ((double)count_boundary_points)/((double)solver->npoints_global)*100.0;
      printf("    Number of immersed boundary points        : %d (%4.1f%%).\n",count_boundary_points,percentage);
      printf("    Immersed body simulation mode             : %s.\n", ib->mode);
    }

  }

  return(0);
}

int InitializePhysics	(	void *	s,
		int	nsims
	)

Initialize the physics

Initialize the physical model for a simulation: Depending on the physical model specified, this function calls the initialization function for that physical model. The latter is responsible for setting all the physics-specific functions that are required by the model.

Parameters

s	Array of simulation objects of type SimulationObject
nsims	number of simulation objects

Definition at line 38 of file InitializePhysics.c.

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

  if (nsims == 0) return 0;

  if (!sim[0].mpi.rank) {
    printf("Initializing physics. Model = \"%s\"\n",sim[0].solver.model);
  }

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

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

    /* Initialize physics-specific functions to NULL */
    solver->ComputeCFL            = NULL;
    solver->ComputeDiffNumber     = NULL;
    solver->FFunction             = NULL;
    solver->dFFunction            = NULL;
    solver->FdFFunction           = NULL;
    solver->GFunction             = NULL;
    solver->HFunction             = NULL;
    solver->SFunction             = NULL;
    solver->UFunction             = NULL;
    solver->JFunction             = NULL;
    solver->KFunction             = NULL;
    solver->Upwind                = NULL;
    solver->UpwinddF              = NULL;
    solver->UpwindFdF             = NULL;
    solver->PreStage              = NULL;
    solver->PostStage             = NULL;
    solver->PreStep               = NULL;
    solver->PostStep              = NULL;
    solver->PrintStep             = NULL;
    solver->PhysicsOutput         = NULL;
    solver->PhysicsInput          = NULL;
    solver->AveragingFunction     = NULL;
    solver->GetLeftEigenvectors   = NULL;
    solver->GetRightEigenvectors  = NULL;
    solver->IBFunction            = NULL;

    if (!strcmp(solver->model,_LINEAR_ADVECTION_DIFFUSION_REACTION_)) {
  
      solver->physics = (LinearADR*) calloc (1,sizeof(LinearADR));
      IERR LinearADRInitialize(solver,mpi); CHECKERR(ierr);

    } else if (!strcmp(solver->model,_BURGERS_)) {

    solver->physics = (Burgers*) calloc (1,sizeof(Burgers));
    IERR BurgersInitialize(solver,mpi); CHECKERR(ierr);

    } else if (!strcmp(solver->model,_FP_DOUBLE_WELL_)) {
  
      solver->physics = (FPDoubleWell*) calloc (1,sizeof(FPDoubleWell));
      IERR FPDoubleWellInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_FP_POWER_SYSTEM_)) {
  
      solver->physics = (FPPowerSystem*) calloc (1,sizeof(FPPowerSystem));
      IERR FPPowerSystemInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_FP_POWER_SYSTEM_1BUS_)) {
  
      solver->physics = (FPPowerSystem1Bus*) calloc (1,sizeof(FPPowerSystem1Bus));
      IERR FPPowerSystem1BusInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_FP_POWER_SYSTEM_3BUS_)) {
  
      solver->physics = (FPPowerSystem3Bus*) calloc (1,sizeof(FPPowerSystem3Bus));
      IERR FPPowerSystem3BusInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_EULER_1D_)) {
  
      solver->physics = (Euler1D*) calloc (1,sizeof(Euler1D));
      IERR Euler1DInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_EULER_2D_)) {
  
      solver->physics = (Euler2D*) calloc (1,sizeof(Euler2D));
      IERR Euler2DInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_NAVIER_STOKES_2D_)) {
  
      solver->physics = (NavierStokes2D*) calloc (1,sizeof(NavierStokes2D));
      IERR NavierStokes2DInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_NAVIER_STOKES_3D_)) {
  
      solver->physics = (NavierStokes3D*) calloc (1,sizeof(NavierStokes3D));
      IERR NavierStokes3DInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_NUMA2D_)) {
  
      solver->physics = (Numa2D*) calloc (1,sizeof(Numa2D));
      IERR Numa2DInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_NUMA3D_)) {
  
      solver->physics = (Numa3D*) calloc (1,sizeof(Numa3D));
      IERR Numa3DInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_SHALLOW_WATER_1D_)) {
  
      solver->physics = (ShallowWater1D*) calloc (1,sizeof(ShallowWater1D));
      IERR ShallowWater1DInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_SHALLOW_WATER_2D_)) {
  
      solver->physics = (ShallowWater2D*) calloc (1,sizeof(ShallowWater2D));
      IERR ShallowWater2DInitialize(solver,mpi); CHECKERR(ierr);
  
    } else if (!strcmp(solver->model,_VLASOV_)) {
  
      solver->physics = (Vlasov*) calloc (1,sizeof(Vlasov));
      IERR VlasovInitialize(solver,mpi); CHECKERR(ierr);
  
    }else {
  
      fprintf(stderr,"Error (domain %d): %s is not a supported physical model.\n",
              ns, solver->model);
      return(1);
  
    }
  
    /* some checks */
    if ( ( (solver->GetLeftEigenvectors == NULL) || (solver->GetRightEigenvectors == NULL) )
        && (!strcmp(solver->interp_type,_CHARACTERISTIC_)) && (solver->nvars > 1) ) {
      if (!mpi->rank) {
        fprintf(stderr,"Error (domain %d): Interpolation type is defined as characteristic ", ns);
        fprintf(stderr,"but physics initializations returned NULL pointers for ");
        fprintf(stderr,"Get(Left,Right)Eigenvectors needed for characteristic-based ");
        fprintf(stderr,"reconstruction.\n");
      }
      return(1);
    }
  
    if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
      if ((!solver->dFFunction) || (!solver->UpwinddF)) {
        if (!mpi->rank) {
          fprintf(stderr,"Error (domain %d): Splitting of hyperbolic flux requires a dFFunction ", ns);
          fprintf(stderr,"and its upwinding function UpwinddF.\n");
          fprintf(stderr,"Error: f(u) = [f(u) - df(u)] + df(u).\n");
          fprintf(stderr,"Error: dFFunction or UpwinddF (or both) is (are) NULL.\n");
        }
        return(1);
      }
      if (solver->FdFFunction && solver->UpwindFdF) solver->flag_fdf_specified = 1;
      else                                          solver->flag_fdf_specified = 0;
    }
  
    if ((solver->IBFunction == NULL) && (solver->flag_ib)) {
      if (!mpi->rank) {
        fprintf(stderr,"Error in InitializePhysics() (domain %d): Physical model %s does not yet have an immersed boundary treatment.\n",
                ns, solver->model);
      }
      return(1);
    }

  }

  return(0);
}

int InitializePhysicsData	(	void *	s,
		int	idx,
		int	nsims,
		int *	dim_data
	)

Initialize the physics data

For each simulation object, call the physics-specific function to read in any physics data that is not a part of the solution vector.

Parameters

s	Simulation object of type SimulationObject
idx	Index of this simulation object
nsims	Total number of simuations
dim_data	Dimenions of physics-specific data

Definition at line 12 of file InitializePhysicsData.c.

{
  SimulationObject *sim     = (SimulationObject*) s;
  HyPar            *solver  = &(sim->solver);
  MPIVariables     *mpi     = &(sim->mpi);

  if (solver->PhysicsInput) {
    int ierr = solver->PhysicsInput(solver, mpi, idx, nsims, dim_data);
    if (ierr) {
      fprintf(stderr, "Error in InitializePhysicsData():\n");
      fprintf(stderr, "  solver->PhysicsInput() returned error %d on rank %d\n", 
              ierr, mpi->rank);
      return ierr;
    }
  }

  return 0;
}

int InitializeSolvers	(	void *	s,
		int	nsims
	)

Initialize the solvers

This function initializes all solvers-specific function pointers depending on user input. The specific functions used for spatial discretization, time integration, and solution output are set here.

Parameters

s	Array of simulation objects of type SimulationObject
nsims	number of simulation objects

Definition at line 52 of file InitializeSolvers.c.

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

  if (nsims == 0) return 0;

  if (!sim[0].mpi.rank) {
    printf("Initializing solvers.\n");
  }

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

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

    solver->ApplyBoundaryConditions = ApplyBoundaryConditions;
    solver->ApplyIBConditions = ApplyIBConditions;
    solver->SourceFunction = SourceFunction;
#if defined(HAVE_CUDA)
    if (solver->use_gpu) {
      solver->HyperbolicFunction = gpuHyperbolicFunction;
    } else {
#endif
      solver->HyperbolicFunction = HyperbolicFunction;
#if defined(HAVE_CUDA)
    }
#endif
    solver->VolumeIntegralFunction      = VolumeIntegral;
    solver->BoundaryIntegralFunction    = BoundaryIntegral;
    solver->CalculateConservationError  = CalculateConservationError;
    solver->NonlinearInterp             = NonLinearInterpolation;

    /* choose the type of parabolic discretization */
    solver->ParabolicFunction         = NULL;
    solver->SecondDerivativePar       = NULL;
    solver->FirstDerivativePar        = NULL;
    solver->InterpolateInterfacesPar  = NULL;

#if defined(HAVE_CUDA)
    if (solver->use_gpu) {

      if (!strcmp(solver->spatial_type_par,_NC_2STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionNC2Stage;

        if (!strcmp(solver->spatial_scheme_par,_FOURTH_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar = gpuFirstDerivativeFourthOrderCentral;
        } else {
          fprintf(stderr,"ERROR (domain %d): scheme %s is not supported on GPU!",
                  ns, solver->spatial_scheme_par);
          return 1;
        }
  
      }

    } else {
#endif

      if (!strcmp(solver->spatial_type_par,_NC_1STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionNC1Stage;
        if (!strcmp(solver->spatial_scheme_par,_SECOND_ORDER_CENTRAL_)) {
          solver->SecondDerivativePar      = SecondDerivativeSecondOrderCentral;
        } else if (!strcmp(solver->spatial_scheme_par,_FOURTH_ORDER_CENTRAL_)) {
          solver->SecondDerivativePar      = SecondDerivativeFourthOrderCentral;
        } else {
          fprintf(stderr,"Error (domain %d): %s is not a supported ",
                  ns, solver->spatial_scheme_par);
          fprintf(stderr,"spatial scheme of type %s for the parabolic terms.\n",
                  solver->spatial_type_par);
        }
  
      } else if (!strcmp(solver->spatial_type_par,_NC_2STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionNC2Stage;
        if (!strcmp(solver->spatial_scheme_par,_SECOND_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar       = FirstDerivativeFirstOrder;
          /* why first order? see ParabolicFunctionNC2Stage.c. 2nd order central
             approximation to the 2nd derivative can be expressed as a conjugation
             of 1st order approximations to the 1st derivative (one forward and
             one backward) -- this prevents odd-even decoupling */
        } else if (!strcmp(solver->spatial_scheme_par,_FOURTH_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar       = FirstDerivativeFourthOrderCentral;
          /* why 4th order? I could not derive the decomposition of the
             4th order central approximation to the 2nd derivative! Some problems
             may show odd-even decoupling */
        } else {
          fprintf(stderr,"Error (domain %d): %s is not a supported ",
                  ns, solver->spatial_scheme_par);
          fprintf(stderr,"spatial scheme of type %s for the parabolic terms.\n",
                solver->spatial_type_par);
        }
  
      } else if (!strcmp(solver->spatial_type_par,_NC_1_5STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionNC1_5Stage;
        if (!strcmp(solver->spatial_scheme_par,_SECOND_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar       = FirstDerivativeSecondOrderCentral;
          solver->SecondDerivativePar      = SecondDerivativeSecondOrderCentral;
        } else if (!strcmp(solver->spatial_scheme_par,_FOURTH_ORDER_CENTRAL_)) {
          solver->FirstDerivativePar       = FirstDerivativeFourthOrderCentral;
          solver->SecondDerivativePar      = SecondDerivativeFourthOrderCentral;
        } else {
          fprintf(stderr,"Error (domain %d): %s is not a supported ",
                  ns, solver->spatial_scheme_par);
          fprintf(stderr,"spatial scheme of type %s for the parabolic terms.\n",
                solver->spatial_type_par);
        }
  
      } else if (!strcmp(solver->spatial_type_par,_CONS_1STAGE_)) {
  
        solver->ParabolicFunction = ParabolicFunctionCons1Stage;
        if (!strcmp(solver->spatial_scheme_par,_SECOND_ORDER_CENTRAL_)) {
          solver->InterpolateInterfacesPar = Interp2PrimSecondOrder;
        } else {
          fprintf(stderr,"Error (domain %d): %s is not a supported ",
                  ns, solver->spatial_scheme_par);
          fprintf(stderr,"spatial scheme of type %s for the parabolic terms.\n",
                solver->spatial_type_par);
        }
  
      } else {
  
        fprintf(stderr,"Error (domain %d): %s is not a supported ",
                ns, solver->spatial_type_par);
        fprintf(stderr,"spatial discretization type for the parabolic terms.\n");
        return(1);
  
      }

#if defined(HAVE_CUDA)
    }
#endif

    /* Spatial interpolation for hyperbolic term */
    solver->interp                = NULL;
    solver->compact               = NULL;
    solver->lusolver              = NULL;
    solver->SetInterpLimiterVar   = NULL;
    solver->flag_nonlinearinterp  = 1;
    if (strcmp(solver->interp_type,_CHARACTERISTIC_) && strcmp(solver->interp_type,_COMPONENTS_)) {
      fprintf(stderr,"Error in InitializeSolvers() (domain %d): %s is not a ",
              ns, solver->interp_type);
      fprintf(stderr,"supported interpolation type.\n");
      return(1);
    }

#if defined(HAVE_CUDA)
    if (solver->use_gpu) {

      if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_WENO_)) {
  
        /* Fifth order WENO scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          fprintf(stderr,
                  "Error (domain %d): characteristic-based WENO5 is not yet implemented on GPUs.\n",
                  ns );
          return 1;
        } else {
          solver->InterpolateInterfacesHyp = gpuInterp1PrimFifthOrderWENO;
        }
        solver->interp = (WENOParameters*) calloc(1,sizeof(WENOParameters));
        IERR WENOInitialize(solver,mpi,solver->spatial_scheme_hyp,solver->interp_type); CHECKERR(ierr);
        solver->flag_nonlinearinterp = !(((WENOParameters*)solver->interp)->no_limiting);
  
      } else {

        fprintf(stderr,
                "Error (domain %d): %s is a not a supported spatial interpolation scheme on GPUs.\n",
                ns, solver->spatial_scheme_hyp);
        return 1;
      }

    } else {
#endif

      if (!strcmp(solver->spatial_scheme_hyp,_FIRST_ORDER_UPWIND_)) {
  
        /* First order upwind scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFirstOrderUpwindChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFirstOrderUpwind;
        }
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_SECOND_ORDER_CENTRAL_)) {
  
        /* Second order central scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimSecondOrderCentralChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimSecondOrderCentral;
        }
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_SECOND_ORDER_MUSCL_)) {
  
        /* Second order MUSCL scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimSecondOrderMUSCLChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimSecondOrderMUSCL;
        }
        solver->interp = (MUSCLParameters*) calloc(1,sizeof(MUSCLParameters));
        IERR MUSCLInitialize(solver,mpi); CHECKERR(ierr);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_THIRD_ORDER_MUSCL_)) {
  
        /* Third order MUSCL scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimThirdOrderMUSCLChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimThirdOrderMUSCL;
        }
        solver->interp = (MUSCLParameters*) calloc(1,sizeof(MUSCLParameters));
        IERR MUSCLInitialize(solver,mpi); CHECKERR(ierr);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FOURTH_ORDER_CENTRAL_)) {
  
        /* Fourth order central scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFourthOrderCentralChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFourthOrderCentral;
        }
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_UPWIND_)) {
  
        /* Fifth order upwind scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderUpwindChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderUpwind;
        }
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_COMPACT_UPWIND_)) {
  
        /* Fifth order compact upwind scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderCompactUpwindChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderCompactUpwind;
        }
        solver->compact = (CompactScheme*) calloc(1,sizeof(CompactScheme));
        IERR CompactSchemeInitialize(solver,mpi,solver->interp_type);
        solver->lusolver = (TridiagLU*) calloc (1,sizeof(TridiagLU));
        IERR tridiagLUInit(solver->lusolver,&mpi->world);CHECKERR(ierr);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_WENO_)) {
  
        /* Fifth order WENO scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderWENOChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderWENO;
        }
        solver->interp = (WENOParameters*) calloc(1,sizeof(WENOParameters));
        IERR WENOInitialize(solver,mpi,solver->spatial_scheme_hyp,solver->interp_type); CHECKERR(ierr);
        solver->flag_nonlinearinterp = !(((WENOParameters*)solver->interp)->no_limiting);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
  
        /* Fifth order CRWENO scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderCRWENOChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderCRWENO;
        }
        solver->interp = (WENOParameters*) calloc(1,sizeof(WENOParameters));
        IERR WENOInitialize(solver,mpi,solver->spatial_scheme_hyp,solver->interp_type); CHECKERR(ierr);
        solver->flag_nonlinearinterp = !(((WENOParameters*)solver->interp)->no_limiting);
        solver->compact = (CompactScheme*) calloc(1,sizeof(CompactScheme));
        IERR CompactSchemeInitialize(solver,mpi,solver->interp_type);
        solver->lusolver = (TridiagLU*) calloc (1,sizeof(TridiagLU));
        IERR tridiagLUInit(solver->lusolver,&mpi->world);CHECKERR(ierr);
  
      } else if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_HCWENO_)) {
  
        /* Fifth order HCWENO scheme */
        if ((solver->nvars > 1) && (!strcmp(solver->interp_type,_CHARACTERISTIC_))) {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderHCWENOChar;
        } else {
          solver->InterpolateInterfacesHyp = Interp1PrimFifthOrderHCWENO;
        }
        solver->interp = (WENOParameters*) calloc(1,sizeof(WENOParameters));
        IERR WENOInitialize(solver,mpi,solver->spatial_scheme_hyp,solver->interp_type); CHECKERR(ierr);
        solver->flag_nonlinearinterp = !(((WENOParameters*)solver->interp)->no_limiting);
        solver->compact = (CompactScheme*) calloc(1,sizeof(CompactScheme));
        IERR CompactSchemeInitialize(solver,mpi,solver->interp_type);
        solver->lusolver = (TridiagLU*) calloc (1,sizeof(TridiagLU));
        IERR tridiagLUInit(solver->lusolver,&mpi->world);CHECKERR(ierr);
  
      } else {

        fprintf(stderr,"Error (domain %d): %s is a not a supported spatial interpolation scheme.\n",
                ns, solver->spatial_scheme_hyp);
        return(1);
      }

#if defined(HAVE_CUDA)
    }
#endif

    /* Time integration */
    solver->time_integrator = NULL;
#ifdef with_petsc
    if (solver->use_petscTS) {
      /* dummy -- not used */
      solver->TimeIntegrate = TimeForwardEuler;
      solver->msti = NULL;
    } else {
      if (!strcmp(solver->time_scheme,_FORWARD_EULER_)) {
        solver->TimeIntegrate = TimeForwardEuler;
        solver->msti = NULL;
      } else if (!strcmp(solver->time_scheme,_RK_)) {
        solver->TimeIntegrate = TimeRK;
        solver->msti = (ExplicitRKParameters*) calloc (1,sizeof(ExplicitRKParameters));
        IERR TimeExplicitRKInitialize(solver->time_scheme,solver->time_scheme_type,
                                      solver->msti,mpi); CHECKERR(ierr);
      } else if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
        solver->TimeIntegrate = TimeGLMGEE;
        solver->msti = (GLMGEEParameters*) calloc (1,sizeof(GLMGEEParameters));
        IERR TimeGLMGEEInitialize(solver->time_scheme,solver->time_scheme_type,
                                  solver->msti,mpi); CHECKERR(ierr);
      } else {
        fprintf(stderr,"Error (domain %d): %s is a not a supported time-integration scheme.\n",
                ns, solver->time_scheme);
        return(1);
      }
    }
#else
    if (!strcmp(solver->time_scheme,_FORWARD_EULER_)) {
      solver->TimeIntegrate = TimeForwardEuler;
      solver->msti = NULL;
    } else if (!strcmp(solver->time_scheme,_RK_)) {
      solver->TimeIntegrate = TimeRK;
      solver->msti = (ExplicitRKParameters*) calloc (1,sizeof(ExplicitRKParameters));
      IERR TimeExplicitRKInitialize(solver->time_scheme,solver->time_scheme_type,
                                    solver->msti,mpi); CHECKERR(ierr);
    } else if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
      solver->TimeIntegrate = TimeGLMGEE;
      solver->msti = (GLMGEEParameters*) calloc (1,sizeof(GLMGEEParameters));
      IERR TimeGLMGEEInitialize(solver->time_scheme,solver->time_scheme_type,
                                solver->msti,mpi); CHECKERR(ierr);
    } else {
      fprintf(stderr,"Error (domain %d): %s is a not a supported time-integration scheme.\n",
              ns, solver->time_scheme);
      return(1);
    }
#endif

    /* Solution output function */
    solver->WriteOutput    = NULL; /* default - no output */
    solver->filename_index = NULL;
    strcpy(solver->op_fname_root, "op");
#ifdef with_librom
    strcpy(solver->op_rom_fname_root, "op_rom");
#endif
    strcpy(solver->aux_op_fname_root, "ts0");
    if (!strcmp(solver->output_mode,"serial")) {
      solver->index_length = 5;
      solver->filename_index = (char*) calloc (solver->index_length+1,sizeof(char));
      int i; for (i=0; i<solver->index_length; i++) solver->filename_index[i] = '0';
      solver->filename_index[solver->index_length] = (char) 0;
      if (!strcmp(solver->op_file_format,"text")) {
        solver->WriteOutput = WriteText;
        strcpy(solver->solnfilename_extn,".dat");
      } else if (!strcmp(solver->op_file_format,"tecplot2d")) {
        solver->WriteOutput = WriteTecplot2D;
        strcpy(solver->solnfilename_extn,".dat");
      } else if (!strcmp(solver->op_file_format,"tecplot3d")) {
        solver->WriteOutput = WriteTecplot3D;
        strcpy(solver->solnfilename_extn,".dat");
      } else if ((!strcmp(solver->op_file_format,"binary")) || (!strcmp(solver->op_file_format,"bin"))) {
        solver->WriteOutput = WriteBinary;
        strcpy(solver->solnfilename_extn,".bin");
      } else if (!strcmp(solver->op_file_format,"none")) {
        solver->WriteOutput = NULL;
      } else {
        fprintf(stderr,"Error (domain %d): %s is not a supported file format.\n",
                ns, solver->op_file_format);
        return(1);
      }
      if ((!strcmp(solver->op_overwrite,"no")) && solver->restart_iter) {
        /* if it's a restart run, fast-forward the filename */
        int t;
        for (t=0; t<solver->restart_iter; t++)
          if ((t+1)%solver->file_op_iter == 0) IncrementFilenameIndex(solver->filename_index,solver->index_length);
      }
    } else if (!strcmp(solver->output_mode,"parallel")) {
      if (!strcmp(solver->op_file_format,"none")) solver->WriteOutput = NULL;
      else {
        /* only binary file writing supported in parallel mode */
        /* use post-processing scripts to convert              */
        solver->WriteOutput = WriteBinary;
        strcpy(solver->solnfilename_extn,".bin");
      }
    } else {
      fprintf(stderr,"Error (domain %d): %s is not a supported output mode.\n",
              ns, solver->output_mode);
      fprintf(stderr,"Should be \"serial\" or \"parallel\".    \n");
      return(1);
    }

    /* Solution plotting function */
    strcpy(solver->plotfilename_extn,".png");
#ifdef with_python
    solver->py_plt_func = NULL;
    solver->py_plt_func_args = NULL;
    {
      char python_plotting_fname[_MAX_STRING_SIZE_] = "plotSolution";
      PyObject* py_plot_name = PyUnicode_DecodeFSDefault(python_plotting_fname);
      PyObject* py_plot_module = PyImport_Import(py_plot_name);
      Py_DECREF(py_plot_name);
      if (py_plot_module) {
        solver->py_plt_func = PyObject_GetAttrString(py_plot_module, "plotSolution");
        if (!solver->py_plt_func) {
          if (!mpi->rank) {
            printf("Unable to load plotSolution function from Python module.\n");
          }
        } else {
          if (!mpi->rank) {
            printf("Loaded Python module for plotting.\n");
            printf("Loaded plotSolution function from Python module.\n");
          }
        }
      } else {
        if (!mpi->rank) {
          printf("Unable to load Python module for plotting.\n");
        }
      }
    }
#endif

  }

  return(0);
}

int Cleanup	(	void *	s,
		int	nsims
	)

Clean up: deallocate all arrays and objects

Cleans up and frees the memory after the completion of the simulation.

Parameters

s	Array of simulation objects of type SimulationObject
nsims	number of simulation objects

Definition at line 39 of file Cleanup.c.

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

  if (nsims == 0) return 0;

  if (!sim[0].mpi.rank) {
    printf("Deallocating arrays.\n");
  }

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

    if (sim[ns].is_barebones == 1) {
      fprintf(stderr, "Error in Cleanup(): object is barebones type.\n");
      return 1;
    }

    HyPar* solver = &(sim[ns].solver);
    MPIVariables* mpi = &(sim[ns].mpi);
    DomainBoundary* boundary = (DomainBoundary*) solver->boundary;
    int i;

    /* Clean up boundary zones */
    for (i = 0; i < solver->nBoundaryZones; i++) {
#if defined(HAVE_CUDA)
      BCCleanup(&boundary[i], solver->use_gpu);
#else
      BCCleanup(&boundary[i], 0);
#endif
    }
    free(solver->boundary);

    /* Clean up immersed boundaries */
    if (solver->flag_ib) {
      IERR IBCleanup(solver->ib);
      free(solver->ib);
    }

    /* Clean up any allocations in physical model */
    if (!strcmp(solver->model,_LINEAR_ADVECTION_DIFFUSION_REACTION_)) {
      IERR LinearADRCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_FP_DOUBLE_WELL_)) {
      IERR FPDoubleWellCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_FP_POWER_SYSTEM_)) {
      IERR FPPowerSystemCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_FP_POWER_SYSTEM_1BUS_)) {
      IERR FPPowerSystem1BusCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_FP_POWER_SYSTEM_3BUS_)) {
      IERR FPPowerSystem3BusCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_EULER_1D_)) {
      IERR Euler1DCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_EULER_2D_)) {
      IERR Euler2DCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_NAVIER_STOKES_2D_)) {
      IERR NavierStokes2DCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_NAVIER_STOKES_3D_)) {
      IERR NavierStokes3DCleanup(solver->physics); CHECKERR(ierr);
#if defined(HAVE_CUDA)
      if (solver->use_gpu) {
        IERR gpuNavierStokes3DCleanup(solver->physics); CHECKERR(ierr);
      }
#endif
    } else if (!strcmp(solver->model,_NUMA2D_)) {
      IERR Numa2DCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_NUMA3D_)) {
      IERR Numa3DCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_SHALLOW_WATER_1D_)) {
      IERR ShallowWater1DCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_SHALLOW_WATER_2D_)) {
      IERR ShallowWater2DCleanup(solver->physics); CHECKERR(ierr);
    } else if (!strcmp(solver->model,_VLASOV_)) {
      IERR VlasovCleanup(solver->physics); CHECKERR(ierr);
    }
    free(solver->physics);

    /* Clean up any allocations from time-integration */
#ifdef with_petsc
    if (!solver->use_petscTS) {
      if (!strcmp(solver->time_scheme,_RK_)) {
        IERR TimeExplicitRKCleanup(solver->msti); CHECKERR(ierr);
        free(solver->msti);
      } else if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
        IERR TimeGLMGEECleanup(solver->msti); CHECKERR(ierr);
        free(solver->msti);
      }
    }
#else
    if (!strcmp(solver->time_scheme,_RK_)) {
      IERR TimeExplicitRKCleanup(solver->msti); CHECKERR(ierr);
      free(solver->msti);
    } else if (!strcmp(solver->time_scheme,_GLM_GEE_)) {
      IERR TimeGLMGEECleanup(solver->msti); CHECKERR(ierr);
      free(solver->msti);
    }
#endif

    /* Clean up any spatial reconstruction related allocations */
    if (   (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_WENO_  ))
        || (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_))
        || (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_HCWENO_)) ) {
#if defined(HAVE_CUDA)
      IERR WENOCleanup(solver->interp, solver->use_gpu); CHECKERR(ierr);
#else
      IERR WENOCleanup(solver->interp, 0); CHECKERR(ierr);
#endif
    }
    if (solver->interp)   free(solver->interp);
    if (   (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_COMPACT_UPWIND_ ))
        || (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_         ))
        || (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_HCWENO_         )) ) {
      IERR CompactSchemeCleanup(solver->compact); CHECKERR(ierr);
    }
    if (solver->compact)  free(solver->compact);
    if (solver->lusolver) free(solver->lusolver);

    /* Free the communicators created */
    IERR MPIFreeCommunicators(solver->ndims,mpi); CHECKERR(ierr);

    /* These variables are allocated in Initialize.c */
    free(solver->dim_global);
    free(solver->dim_global_ex);
    free(solver->dim_local);
    free(solver->index);
    free(solver->u);
#ifdef with_petsc
    if (solver->u0)     free(solver->u0);
    if (solver->uref)   free(solver->uref);
    if (solver->rhsref) free(solver->rhsref);
    if (solver->rhs)    free(solver->rhs);
#endif
#ifdef with_librom
    free(solver->u_rom_predicted);
#endif
    free(solver->iblank);
    free(solver->x);
    free(solver->dxinv);
    free(solver->isPeriodic);
    free(mpi->iproc);
    free(mpi->ip);
    free(mpi->is);
    free(mpi->ie);
    free(mpi->bcperiodic);
    free(mpi->sendbuf);
    free(mpi->recvbuf);
    free(solver->VolumeIntegral);
    free(solver->VolumeIntegralInitial);
    free(solver->TotalBoundaryIntegral);
    free(solver->ConservationError);
    free(solver->stride_with_ghosts);
    free(solver->stride_without_ghosts);

#if defined(HAVE_CUDA)
    if (solver->use_gpu) {
      gpuFree(solver->hyp);
      gpuFree(solver->par);
      gpuFree(solver->source);
      gpuFree(solver->uC);
      gpuFree(solver->fluxC);
      gpuFree(solver->Deriv1);
      gpuFree(solver->Deriv2);
      gpuFree(solver->fluxI);
      gpuFree(solver->uL);
      gpuFree(solver->uR);
      gpuFree(solver->fL);
      gpuFree(solver->fR);
      gpuFree(solver->StageBoundaryBuffer);
      gpuFree(solver->StageBoundaryIntegral);
      gpuFree(solver->StepBoundaryIntegral);

      gpuFree(solver->gpu_dim_local);
      gpuFree(solver->gpu_iblank);
      gpuFree(solver->gpu_x);
      gpuFree(solver->gpu_dxinv);
      gpuFree(solver->gpu_u);
    } else {
#endif
      free(solver->hyp);
      free(solver->par);
      free(solver->source);
      free(solver->uC);
      free(solver->fluxC);
      free(solver->Deriv1);
      free(solver->Deriv2);
      free(solver->fluxI);
      free(solver->uL);
      free(solver->uR);
      free(solver->fL);
      free(solver->fR);
      free(solver->StageBoundaryIntegral);
      free(solver->StepBoundaryIntegral);
#if defined(HAVE_CUDA)
    }
#endif

    if (solver->filename_index) free(solver->filename_index);

  }

  return(0);
}

void SimWriteErrors	(	void *	s,
		int	nsims,
		int	rank,
		double	solver_runtime,
		double	main_runtime
	)

Write errors for each simulation

Writes out the errors and other data for each simulation.

Parameters

s	Array of simulations of type SimulationObject
nsims	Number of simulations
rank	MPI rank of this process
solver_runtime	Measured runtime of solver
main_runtime	Measured total runtime

Definition at line 18 of file SimulationWriteErrors.c.

{
  SimulationObject* sim = (SimulationObject*) s;
  int n;

  if (!rank) {

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

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

      char err_fname[_MAX_STRING_SIZE_],
           cons_fname[_MAX_STRING_SIZE_],
           fc_fname[_MAX_STRING_SIZE_];
      strcpy(err_fname,"errors");
      strcpy(cons_fname,"conservation");
      strcpy(fc_fname,"function_counts");
#ifdef with_librom
      char rom_diff_fname[_MAX_STRING_SIZE_]; 
      strcpy(rom_diff_fname,"pde_rom_diff");
#endif


      if (nsims > 1) {

        strcat(err_fname,"_");
        strcat(cons_fname,"_");
        strcat(fc_fname,"_");
#ifdef with_librom
        strcat(rom_diff_fname,"_");
#endif

        char index[_MAX_STRING_SIZE_];
        GetStringFromInteger(n, index, (int)log10(nsims)+1);

        strcat(err_fname,index);
        strcat(cons_fname,index);
        strcat(fc_fname,index);
#ifdef with_librom
        strcat(rom_diff_fname,index);
#endif
      }

      strcat(err_fname,".dat");
      strcat(cons_fname,".dat");
      strcat(fc_fname,".dat");
#ifdef with_librom
      strcat(rom_diff_fname,".dat");
#endif

      FILE *out; 
      /* write out solution errors and wall times to file */
      out = fopen(err_fname,"w");
      for (int d=0; d<sim[n].solver.ndims; d++) fprintf(out,"%4d ",sim[n].solver.dim_global[d]);
      for (int d=0; d<sim[n].solver.ndims; d++) fprintf(out,"%4d ",sim[n].mpi.iproc[d]);
      fprintf(out,"%1.16E  ",sim[n].solver.dt);
      fprintf(out,"%1.16E %1.16E %1.16E   ",sim[n].solver.error[0],sim[n].solver.error[1],sim[n].solver.error[2]);
      fprintf(out,"%1.16E %1.16E\n",solver_runtime,main_runtime);
      fclose(out);
      /* write out conservation errors to file */
      out = fopen(cons_fname,"w");
      for (int d=0; d<sim[n].solver.ndims; d++) fprintf(out,"%4d ",sim[n].solver.dim_global[d]);
      for (int d=0; d<sim[n].solver.ndims; d++) fprintf(out,"%4d ",sim[n].mpi.iproc[d]);
      fprintf(out,"%1.16E  ",sim[n].solver.dt);
      for (int d=0; d<sim[n].solver.nvars; d++) fprintf(out,"%1.16E ",sim[n].solver.ConservationError[d]);
      fprintf(out,"\n");
      fclose(out);
      /* write out function call counts to file */
      out = fopen(fc_fname,"w");
      fprintf(out,"%d\n",sim[n].solver.n_iter);
      fprintf(out,"%d\n",sim[n].solver.count_hyp);
      fprintf(out,"%d\n",sim[n].solver.count_par);
      fprintf(out,"%d\n",sim[n].solver.count_sou);
#ifdef with_petsc
      fprintf(out,"%d\n",sim[n].solver.count_RHSFunction);
      fprintf(out,"%d\n",sim[n].solver.count_IFunction);
      fprintf(out,"%d\n",sim[n].solver.count_IJacobian);
      fprintf(out,"%d\n",sim[n].solver.count_IJacFunction);
#endif
      fclose(out);
#ifdef with_librom
      /* write out solution errors and wall times to file */
      if (sim[n].solver.rom_diff_norms[0] >= 0) {
        out = fopen(rom_diff_fname,"w");
        for (int d=0; d<sim[n].solver.ndims; d++) fprintf(out,"%4d ",sim[n].solver.dim_global[d]);
        for (int d=0; d<sim[n].solver.ndims; d++) fprintf(out,"%4d ",sim[n].mpi.iproc[d]);
        fprintf(out,"%1.16E  ",sim[n].solver.dt);
        fprintf(out,"%1.16E %1.16E %1.16E   ",sim[n].solver.rom_diff_norms[0],sim[n].solver.rom_diff_norms[1],sim[n].solver.rom_diff_norms[2]);
        fprintf(out,"%1.16E %1.16E\n",solver_runtime,main_runtime);
        fclose(out);
      }
#endif

      /* print solution errors, conservation errors, and wall times to screen */
      if (sim[n].solver.error[0] >= 0) {
        printf("Computed errors for domain %d:\n", n);
        printf("  L1         Error           : %1.16E\n",sim[n].solver.error[0]);
        printf("  L2         Error           : %1.16E\n",sim[n].solver.error[1]);
        printf("  Linfinity  Error           : %1.16E\n",sim[n].solver.error[2]);
      }
      if (!strcmp(sim[n].solver.ConservationCheck,"yes")) {
        printf("Conservation Errors:\n");
        for (int d=0; d<sim[n].solver.nvars; d++) printf("\t%1.16E\n",sim[n].solver.ConservationError[d]);
      }
#ifdef with_librom
      if (sim[n].solver.rom_diff_norms[0] >= 0) {
        printf("Norms of the diff between ROM and PDE solutions for domain %d:\n", n);
        printf("  L1         Norm            : %1.16E\n",sim[n].solver.rom_diff_norms[0]);
        printf("  L2         Norm            : %1.16E\n",sim[n].solver.rom_diff_norms[1]);
        printf("  Linfinity  Norm            : %1.16E\n",sim[n].solver.rom_diff_norms[2]);
      }
#endif

    }

    printf("Solver runtime (in seconds): %1.16E\n",solver_runtime);
    printf("Total  runtime (in seconds): %1.16E\n",main_runtime);
    if (nsims > 1) printf("\n");

  }

  return;
}

int SolvePETSc	(	void *	s,
		int	nsims,
		int	rank,
		int	nproc
	)

Integrate in time with PETSc.

Solve the PDE using PETSc TS

This function integrates the semi-discrete ODE (obtained from discretizing the PDE in space) using the time-integration module of PETSc (https://petsc.org/release/docs/manualpages/TS/index.html). The time-integration context is set up using the parameters specified in the input file. However, they can also be specified using PETSc's command line inputs.

See PETSc's documentation and examples for more details on how to use its TS module. All functions and data types whose names start with Vec, Mat, PC, KSP, SNES, and TS are PETSc functions - refer to the PETSc documentation (usually googling with the function name shows the man page for that function on PETSc's website).

Parameters

s	Array of simulation objects of type SimulationObject
nsims	number of simulation objects
rank	MPI rank of this process
nproc	Number of MPI processes

Definition at line 50 of file SolvePETSc.cpp.

{
  SimulationObject* sim = (SimulationObject*) s;

  DM              dm; /* data management object */
  TS              ts; /* time integration object */
  Vec             Y,Z; /* PETSc solution vectors */
  Mat             A, B; /* Jacobian and preconditioning matrices */
  MatFDColoring   fdcoloring; /* coloring for sparse Jacobian computation */
  TSType          time_scheme; /* time integration method */
  TSProblemType   ptype; /* problem type - nonlinear or linear */

  int flag_mat_a = 0, 
      flag_mat_b = 0, 
      flag_fdcoloring = 0,
      iAuxSize = 0, i;

  PetscFunctionBegin;

  /* Register custom time-integration methods, if specified */
  PetscRegisterTIMethods(rank);
  if(!rank) printf("Setting up PETSc time integration... \n");

  /* create and set a PETSc context */
  PETScContext context;

  context.rank = rank;
  context.nproc = nproc;

  context.simobj = sim;
  context.nsims = nsims;

  /* default: everything explicit */
  context.flag_hyperbolic     = _EXPLICIT_; 
  context.flag_hyperbolic_f   = _EXPLICIT_; 
  context.flag_hyperbolic_df  = _EXPLICIT_; 
  context.flag_parabolic      = _EXPLICIT_; 
  context.flag_source         = _EXPLICIT_; 

  context.tic = 0;
  context.flag_is_linear = 0;
  context.globalDOF.clear();
  context.points.clear();
  context.ti_runtime = 0.0;
  context.waqt = 0.0;
  context.dt = sim[0].solver.dt;
  context.stage_times.clear();
  context.stage_index = 0;

#ifdef with_librom
  if (!rank) printf("Setting up libROM interface.\n");
  context.rom_interface = new libROMInterface( sim, 
                                               nsims, 
                                               rank, 
                                               nproc, 
                                               sim[0].solver.dt );
  context.rom_mode = ((libROMInterface*)context.rom_interface)->mode();
  context.op_times_arr.clear();
#endif

#ifdef with_librom
  if (      (context.rom_mode == _ROM_MODE_TRAIN_) 
        ||  (context.rom_mode == _ROM_MODE_INITIAL_GUESS_ )
        ||  (context.rom_mode == _ROM_MODE_NONE_ ) ) {

    if (context.rom_mode == _ROM_MODE_INITIAL_GUESS_) {
      ((libROMInterface*)context.rom_interface)->loadROM();
      ((libROMInterface*)context.rom_interface)->projectInitialSolution(sim);
    }
#endif
    PetscCreatePointList(&context);
  
    /* create and initialize PETSc solution vector and other parameters */
    /* PETSc solution vector does not have ghost points */
    VecCreate(MPI_COMM_WORLD,&Y);
    VecSetSizes(Y,context.ndofs,PETSC_DECIDE);
    VecSetUp(Y);
  
    /* copy initial solution to PETSc's vector */
    for (int ns = 0; ns < nsims; ns++) {
      TransferVecToPETSc( sim[ns].solver.u,
                          Y,
                          &context,
                          ns,
                          context.offsets[ns] );
    }
  
    /* Create the global DOF mapping for all the grid points */
    PetscGlobalDOF(&context);
  
    /* Define and initialize the time-integration object */
    TSCreate(MPI_COMM_WORLD,&ts);
    TSSetMaxSteps(ts,sim[0].solver.n_iter);
    TSSetMaxTime(ts,sim[0].solver.dt*sim[0].solver.n_iter);
    TSSetTimeStep(ts,sim[0].solver.dt);
    TSSetTime(ts,context.waqt);
    TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
    TSSetType(ts,TSBEULER);
  
    /* set default time step adaptivity to none */
    TSAdapt adapt;
    TSAdaptType adapt_type = TSADAPTNONE;
    TSGetAdapt(ts,&adapt);
    TSAdaptSetType(adapt,adapt_type);
  
    /* set options from input */
    TSSetFromOptions(ts);

    /* create DM */
    DMShellCreate(MPI_COMM_WORLD, &dm);
    DMShellSetGlobalVector(dm, Y);
    TSSetDM(ts, dm);
  
#ifdef with_librom
    TSAdaptGetType(adapt,&adapt_type);
    if (strcmp(adapt_type, TSADAPTNONE)) {
      if (!rank) printf("Warning: libROM interface not yet implemented for adaptive timestepping.\n");
    }
#endif
  
    /* Define the right and left -hand side functions for each time-integration scheme */
    TSGetType(ts,&time_scheme);
    TSGetProblemType(ts,&ptype);

    if (!strcmp(time_scheme,TSARKIMEX)) {
  
      /* implicit - explicit time integration */
  
      TSSetRHSFunction(ts,nullptr,PetscRHSFunctionIMEX,&context);
      TSSetIFunction  (ts,nullptr,PetscIFunctionIMEX,  &context);
  
      SNES     snes;
      KSP      ksp;
      PC       pc;
      SNESType snestype;
      TSGetSNES(ts,&snes);
      SNESGetType(snes,&snestype);

#ifdef with_librom
      if (context.rom_mode == _ROM_MODE_INITIAL_GUESS_) {
        SNESSetComputeInitialGuess(snes, PetscSetInitialGuessROM, &context);
      }
#endif
  
      context.flag_use_precon = 0;
      PetscOptionsGetBool(  nullptr,nullptr,
                            "-with_pc",
                            (PetscBool*)(&context.flag_use_precon),
                            nullptr );
  
      char precon_mat_type_c_st[_MAX_STRING_SIZE_] = "default";
      PetscOptionsGetString(  nullptr,
                              nullptr,
                              "-pc_matrix_type",
                              precon_mat_type_c_st,
                              _MAX_STRING_SIZE_,
                              nullptr );
      context.precon_matrix_type = std::string(precon_mat_type_c_st);
  
      if (context.flag_use_precon) {

        if (context.precon_matrix_type == "default") {

          /* Matrix-free representation of the Jacobian */
          flag_mat_a = 1;
          MatCreateShell( MPI_COMM_WORLD,
                          context.ndofs,
                          context.ndofs,
                          PETSC_DETERMINE,
                          PETSC_DETERMINE,
                          &context,
                          &A);
          if ((!strcmp(snestype,SNESKSPONLY)) || (ptype == TS_LINEAR)) {
            /* linear problem */
            context.flag_is_linear = 1;
            MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunctionIMEX_Linear);
            SNESSetType(snes,SNESKSPONLY);
          } else {
            /* nonlinear problem */
            context.flag_is_linear = 0;
            context.jfnk_eps = 1e-7;
            PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps,NULL);
            MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunctionIMEX_JFNK);
          }
          MatSetUp(A);
          /* check if Jacobian of the physical model is defined */
          for (int ns = 0; ns < nsims; ns++) {
            if ((!sim[ns].solver.JFunction) && (!sim[ns].solver.KFunction)) {
              if (!rank) {
                fprintf(stderr,"Error in SolvePETSc(): solver->JFunction  or solver->KFunction ");
                fprintf(stderr,"(point-wise Jacobians for hyperbolic or parabolic terms) must ");
                fprintf(stderr,"be defined for preconditioning.\n");
              }
              PetscFunctionReturn(1);
            }
          }
          /* Set up preconditioner matrix */
          flag_mat_b = 1;
          MatCreateAIJ( MPI_COMM_WORLD,
                        context.ndofs, 
                        context.ndofs,
                        PETSC_DETERMINE, 
                        PETSC_DETERMINE,
                        (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
                        2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
                        &B );
          MatSetBlockSize(B,sim[0].solver.nvars);
          /* Set the IJacobian function for TS */
          TSSetIJacobian(ts,A,B,PetscIJacobianIMEX,&context);

        } else if (context.precon_matrix_type == "fd") {

          flag_mat_a = 1;
          MatCreateSNESMF(snes,&A);
          flag_mat_b = 1;
          MatCreateAIJ( MPI_COMM_WORLD,
                        context.ndofs,
                        context.ndofs,
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
                        2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
                        &B);
          MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
          /* Set the Jacobian function for SNES */
          SNESSetJacobian(snes, A, B, SNESComputeJacobianDefault, NULL);

        } else if (context.precon_matrix_type == "colored_fd") {

          int stencil_width = 1;
          PetscOptionsGetInt( NULL,
                              NULL,
                              "-pc_matrix_colored_fd_stencil_width",
                              &stencil_width,
                              NULL );

          flag_mat_a = 1;
          MatCreateSNESMF(snes,&A);
          flag_mat_b = 1;
          MatCreateAIJ( MPI_COMM_WORLD,
                        context.ndofs,
                        context.ndofs,
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
                        2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
                        &B);
          MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
          if (!rank) {
            printf("PETSc:    Setting Jacobian non-zero pattern (stencil width %d).\n",
                    stencil_width );
          }
          PetscJacobianMatNonzeroEntriesImpl(B, stencil_width, &context);

          /* Set the Jacobian function for SNES */
          SNESSetJacobian(snes, A, B, SNESComputeJacobianDefaultColor, NULL);

        } else {

          if (!rank) {
            fprintf(  stderr,"Invalid input for \"-pc_matrix_type\": %s.\n", 
                      context.precon_matrix_type.c_str());
          }
          PetscFunctionReturn(0);

        }

        /* set PC side to right */
        SNESGetKSP(snes,&ksp);
        KSPSetPCSide(ksp, PC_RIGHT);

      } else {

        /* Matrix-free representation of the Jacobian */
        flag_mat_a = 1;
        MatCreateShell( MPI_COMM_WORLD,
                        context.ndofs,
                        context.ndofs,
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        &context,
                        &A);
        if ((!strcmp(snestype,SNESKSPONLY)) || (ptype == TS_LINEAR)) {
          /* linear problem */
          context.flag_is_linear = 1;
          MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunctionIMEX_Linear);
          SNESSetType(snes,SNESKSPONLY);
        } else {
          /* nonlinear problem */
          context.flag_is_linear = 0;
          context.jfnk_eps = 1e-7;
          PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps,NULL);
          MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunctionIMEX_JFNK);
        }
        MatSetUp(A);
        /* Set the RHSJacobian function for TS */
        TSSetIJacobian(ts,A,A,PetscIJacobianIMEX,&context);
        /* Set PC (preconditioner) to none */
        SNESGetKSP(snes,&ksp);
        KSPGetPC(ksp,&pc);
        PCSetType(pc,PCNONE);
      }

      /* read the implicit/explicit flags for each of the terms for IMEX schemes */
      /* default -> hyperbolic - explicit, parabolic and source - implicit       */
      PetscBool flag = PETSC_FALSE;
  
      context.flag_hyperbolic     = _EXPLICIT_; 
      context.flag_hyperbolic_f   = _EXPLICIT_; 
      context.flag_hyperbolic_df  = _IMPLICIT_; 
      context.flag_parabolic      = _IMPLICIT_; 
      context.flag_source         = _IMPLICIT_; 
  
      if (!strcmp(sim[0].solver.SplitHyperbolicFlux,"yes")) {
  
        flag = PETSC_FALSE; 
        PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_f_explicit",&flag,nullptr);
        if (flag == PETSC_TRUE) context.flag_hyperbolic_f = _EXPLICIT_; 
        flag = PETSC_FALSE; 
        PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_f_implicit",&flag,nullptr);
        if (flag == PETSC_TRUE) context.flag_hyperbolic_f = _IMPLICIT_; 
  
        flag = PETSC_FALSE; 
        PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_df_explicit",&flag,nullptr);
        if (flag == PETSC_TRUE) context.flag_hyperbolic_df = _EXPLICIT_; 
        flag = PETSC_FALSE; 
        PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_df_implicit",&flag,nullptr);
        if (flag == PETSC_TRUE) context.flag_hyperbolic_df = _IMPLICIT_; 
  
      } else {
  
        flag = PETSC_FALSE; 
        PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_explicit",&flag,nullptr);
        if (flag == PETSC_TRUE) context.flag_hyperbolic = _EXPLICIT_; 
        flag = PETSC_FALSE; 
        PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_implicit",&flag,nullptr);
        if (flag == PETSC_TRUE) context.flag_hyperbolic = _IMPLICIT_; 
  
      }
  
      flag = PETSC_FALSE; 
      PetscOptionsGetBool(nullptr,nullptr,"-parabolic_explicit",&flag,nullptr);
      if (flag == PETSC_TRUE) context.flag_parabolic = _EXPLICIT_; 
      flag = PETSC_FALSE; 
      PetscOptionsGetBool(nullptr,nullptr,"-parabolic_implicit",&flag,nullptr);
      if (flag == PETSC_TRUE) context.flag_parabolic = _IMPLICIT_; 
  
      flag = PETSC_FALSE; 
      PetscOptionsGetBool(nullptr,nullptr,"-source_explicit",&flag,nullptr);
      if (flag == PETSC_TRUE) context.flag_source = _EXPLICIT_; 
      flag = PETSC_FALSE; 
      PetscOptionsGetBool(nullptr,nullptr,"-source_implicit",&flag,nullptr);
      if (flag == PETSC_TRUE) context.flag_source = _IMPLICIT_; 
  
      flag = PETSC_FALSE;
      PetscOptionsGetBool(nullptr,nullptr,"-ts_arkimex_fully_implicit",&flag,nullptr);
      if (flag == PETSC_TRUE) {
        context.flag_hyperbolic_f   = _IMPLICIT_;
        context.flag_hyperbolic_df  = _IMPLICIT_;
        context.flag_hyperbolic     = _IMPLICIT_;
        context.flag_parabolic      = _IMPLICIT_;
        context.flag_source         = _IMPLICIT_;
      }
  
      /* print out a summary of the treatment of each term */
      if (!rank) {
        printf("Implicit-Explicit time-integration:-\n");
        if (!strcmp(sim[0].solver.SplitHyperbolicFlux,"yes")) {
          if (context.flag_hyperbolic_f == _EXPLICIT_)  printf("Hyperbolic (f-df) term: Explicit\n");
          else                                          printf("Hyperbolic (f-df) term: Implicit\n");
          if (context.flag_hyperbolic_df == _EXPLICIT_) printf("Hyperbolic (df)   term: Explicit\n");
          else                                          printf("Hyperbolic (df)   term: Implicit\n");
        } else {
          if (context.flag_hyperbolic == _EXPLICIT_)    printf("Hyperbolic        term: Explicit\n");
          else                                          printf("Hyperbolic        term: Implicit\n");
        }
        if (context.flag_parabolic == _EXPLICIT_)       printf("Parabolic         term: Explicit\n");
        else                                            printf("Parabolic         term: Implicit\n");
        if (context.flag_source    == _EXPLICIT_)       printf("Source            term: Explicit\n");
        else                                            printf("Source            term: Implicit\n");
      }
  
    } else if (     (!strcmp(time_scheme,TSEULER)) 
                ||  (!strcmp(time_scheme,TSRK   )) 
                ||  (!strcmp(time_scheme,TSSSP  )) ) {
  
      /* Explicit time integration */    
      TSSetRHSFunction(ts,nullptr,PetscRHSFunctionExpl,&context);
  
    } else if (     (!strcmp(time_scheme,TSCN)) 
                ||  (!strcmp(time_scheme,TSBEULER )) ) {
  
  
      /* Implicit time integration */
  
      TSSetIFunction(ts,nullptr,PetscIFunctionImpl,&context);
  
      SNES     snes;
      KSP      ksp;
      PC       pc;
      SNESType snestype;
      TSGetSNES(ts,&snes);
      SNESGetType(snes,&snestype);
  
#ifdef with_librom
      if (context.rom_mode == _ROM_MODE_INITIAL_GUESS_) {
        SNESSetComputeInitialGuess(snes, PetscSetInitialGuessROM, &context);
      }
#endif
  
      context.flag_use_precon = 0;
      PetscOptionsGetBool(  nullptr,
                            nullptr,
                            "-with_pc",
                            (PetscBool*)(&context.flag_use_precon),
                            nullptr );
  
      char precon_mat_type_c_st[_MAX_STRING_SIZE_] = "default";
      PetscOptionsGetString(  nullptr,
                              nullptr,
                              "-pc_matrix_type",
                              precon_mat_type_c_st,
                              _MAX_STRING_SIZE_,
                              nullptr );
      context.precon_matrix_type = std::string(precon_mat_type_c_st);
  
      if (context.flag_use_precon) {

        if (context.precon_matrix_type == "default") {

          /* Matrix-free representation of the Jacobian */
          flag_mat_a = 1;
          MatCreateShell( MPI_COMM_WORLD,
                          context.ndofs,
                          context.ndofs,
                          PETSC_DETERMINE,
                          PETSC_DETERMINE,
                          &context,
                          &A);
          if ((!strcmp(snestype,SNESKSPONLY)) || (ptype == TS_LINEAR)) {
            /* linear problem */
            context.flag_is_linear = 1;
            MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunction_Linear);
            SNESSetType(snes,SNESKSPONLY);
          } else {
            /* nonlinear problem */
            context.flag_is_linear = 0;
            context.jfnk_eps = 1e-7;
            PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps,NULL);
            MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunction_JFNK);
          }
          MatSetUp(A);
          /* check if Jacobian of the physical model is defined */
          for (int ns = 0; ns < nsims; ns++) {
            if ((!sim[ns].solver.JFunction) && (!sim[ns].solver.KFunction)) {
              if (!rank) {
                fprintf(stderr,"Error in SolvePETSc(): solver->JFunction  or solver->KFunction ");
                fprintf(stderr,"(point-wise Jacobians for hyperbolic or parabolic terms) must ");
                fprintf(stderr,"be defined for preconditioning.\n");
              }
              PetscFunctionReturn(1);
            }
          }
          /* Set up preconditioner matrix */
          flag_mat_b = 1;
          MatCreateAIJ( MPI_COMM_WORLD,
                        context.ndofs, 
                        context.ndofs,
                        PETSC_DETERMINE, 
                        PETSC_DETERMINE,
                        (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
                        2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
                        &B );
          MatSetBlockSize(B,sim[0].solver.nvars);
          /* Set the IJacobian function for TS */
          TSSetIJacobian(ts,A,B,PetscIJacobian,&context);

        } else if (context.precon_matrix_type == "fd") {

          flag_mat_a = 1;
          MatCreateSNESMF(snes,&A);
          flag_mat_b = 1;
          MatCreateAIJ( MPI_COMM_WORLD,
                        context.ndofs,
                        context.ndofs,
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
                        2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
                        &B);
          MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
          /* Set the Jacobian function for SNES */
          SNESSetJacobian(snes, A, B, SNESComputeJacobianDefault, NULL);

        } else if (context.precon_matrix_type == "colored_fd") {

          int stencil_width = 1;
          PetscOptionsGetInt( NULL,
                              NULL,
                              "-pc_matrix_colored_fd_stencil_width",
                              &stencil_width,
                              NULL );

          flag_mat_a = 1;
          MatCreateSNESMF(snes,&A);
          flag_mat_b = 1;
          MatCreateAIJ( MPI_COMM_WORLD,
                        context.ndofs,
                        context.ndofs,
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
                        2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
                        &B);
          MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
          if (!rank) {
            printf("PETSc:    Setting Jacobian non-zero pattern (stencil width %d).\n",
                    stencil_width );
          }
          PetscJacobianMatNonzeroEntriesImpl(B, stencil_width, &context);

          /* Set the Jacobian function for SNES */
          SNESSetJacobian(snes, A, B, SNESComputeJacobianDefaultColor, NULL);

        } else {

          if (!rank) {
            fprintf(  stderr,"Invalid input for \"-pc_matrix_type\": %s.\n", 
                      context.precon_matrix_type.c_str());
          }
          PetscFunctionReturn(0);

        }

        /* set PC side to right */
        SNESGetKSP(snes,&ksp);
        KSPSetPCSide(ksp, PC_RIGHT);

      } else {

        /* Matrix-free representation of the Jacobian */
        flag_mat_a = 1;
        MatCreateShell( MPI_COMM_WORLD,
                        context.ndofs,
                        context.ndofs,
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        &context,
                        &A);
        if ((!strcmp(snestype,SNESKSPONLY)) || (ptype == TS_LINEAR)) {
          /* linear problem */
          context.flag_is_linear = 1;
          MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunction_Linear);
          SNESSetType(snes,SNESKSPONLY);
        } else {
          /* nonlinear problem */
          context.flag_is_linear = 0;
          context.jfnk_eps = 1e-7;
          PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps,NULL);
          MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunction_JFNK);
        }
        MatSetUp(A);
        /* Set the RHSJacobian function for TS */
        TSSetIJacobian(ts,A,A,PetscIJacobian,&context);
        /* Set PC (preconditioner) to none */
        SNESGetKSP(snes,&ksp);
        KSPGetPC(ksp,&pc);
        PCSetType(pc,PCNONE);
      }
  
    } else {

      if (!rank) {
        fprintf(stderr, "Time integration type %s is not yet supported.\n", time_scheme);
      }
      PetscFunctionReturn(0);

    }
  
    /* Set pre/post-stage and post-timestep function */
    TSSetPreStep (ts,PetscPreTimeStep );
    TSSetPreStage(ts,PetscPreStage    );
    TSSetPostStage(ts,PetscPostStage  );
    TSSetPostStep(ts,PetscPostTimeStep);
    /* Set solution vector for TS */
    TSSetSolution(ts,Y);
    /* Set it all up */
    TSSetUp(ts);
    /* Set application context */
    TSSetApplicationContext(ts,&context);
  
    if (!rank) {
      if (context.flag_is_linear) printf("SolvePETSc(): Problem type is linear.\n");
      else                        printf("SolvePETSc(): Problem type is nonlinear.\n");
    }
  
    if (!rank) printf("** Starting PETSc time integration **\n");
    context.ti_runtime = 0.0;
    TSSolve(ts,Y);
    if (!rank) {
      printf("** Completed PETSc time integration (Final time: %f), total wctime: %f (seconds) **\n",
              context.waqt, context.ti_runtime );
    }
  
    /* Get the number of time steps */
    for (int ns = 0; ns < nsims; ns++) {
      TSGetStepNumber(ts,&(sim[ns].solver.n_iter));
    }
  
    /* get and write to file any auxiliary solutions */
    char aux_fname_root[4] = "ts0";
    TSGetSolutionComponents(ts,&iAuxSize,NULL);
    if (iAuxSize) {
      if (iAuxSize > 10) iAuxSize = 10;
      if (!rank) printf("Number of auxiliary solutions from time integration: %d\n",iAuxSize);
      VecDuplicate(Y,&Z);
      for (i=0; i<iAuxSize; i++) {
        TSGetSolutionComponents(ts,&i,&Z);
        for (int ns = 0; ns < nsims; ns++) {
          TransferVecFromPETSc(sim[ns].solver.u,Z,&context,ns,context.offsets[ns]);
          WriteArray( sim[ns].solver.ndims,
                      sim[ns].solver.nvars,
                      sim[ns].solver.dim_global,
                      sim[ns].solver.dim_local,
                      sim[ns].solver.ghosts,
                      sim[ns].solver.x,
                      sim[ns].solver.u,
                      &(sim[ns].solver),
                      &(sim[ns].mpi),
                      aux_fname_root );
        }
        aux_fname_root[2]++;
      }
      VecDestroy(&Z);
    }
  
    /* if available, get error estimates */
    PetscTimeError(ts);
  
    /* copy final solution from PETSc's vector */
    for (int ns = 0; ns < nsims; ns++) {
      TransferVecFromPETSc(sim[ns].solver.u,Y,&context,ns,context.offsets[ns]);
    }
  
    /* clean up */
    VecDestroy(&Y);
    if (flag_mat_a) { MatDestroy(&A); }
    if (flag_mat_b) { MatDestroy(&B); }
    if (flag_fdcoloring) { MatFDColoringDestroy(&fdcoloring); }
    TSDestroy(&ts);
    DMDestroy(&dm);
  
    /* write a final solution file, if last iteration did not write one */
    if (context.tic) { 
      for (int ns = 0; ns < nsims; ns++) {
        HyPar* solver = &(sim[ns].solver);
        MPIVariables* mpi = &(sim[ns].mpi);
        if (solver->PhysicsOutput) {
          solver->PhysicsOutput(solver,mpi, context.waqt);
        }
        CalculateError(solver,mpi);
      }
      OutputSolution(sim, nsims, context.waqt); 
    }
    /* calculate error if exact solution has been provided */
    for (int ns = 0; ns < nsims; ns++) {
      CalculateError(&(sim[ns].solver), &(sim[ns].mpi));
    }
  
    PetscCleanup(&context);

#ifdef with_librom
    context.op_times_arr.push_back(context.waqt);

    for (int ns = 0; ns < nsims; ns++) {
      ResetFilenameIndex( sim[ns].solver.filename_index, 
                          sim[ns].solver.index_length );
    }
  
    if (((libROMInterface*)context.rom_interface)->mode() == _ROM_MODE_TRAIN_) {
  
      ((libROMInterface*)context.rom_interface)->train();
      if (!rank) printf("libROM: total training wallclock time: %f (seconds).\n", 
                        ((libROMInterface*)context.rom_interface)->trainWallclockTime() );

      double total_rom_predict_time = 0;
      for (int iter = 0; iter < context.op_times_arr.size(); iter++) {
  
        double waqt = context.op_times_arr[iter];
  
        ((libROMInterface*)context.rom_interface)->predict(sim, waqt);
        if (!rank) printf(  "libROM: Predicted solution at time %1.4e using ROM, wallclock time: %f.\n", 
                            waqt, ((libROMInterface*)context.rom_interface)->predictWallclockTime() );
        total_rom_predict_time += ((libROMInterface*)context.rom_interface)->predictWallclockTime();
  
        /* calculate diff between ROM and PDE solutions */
        if (iter == (context.op_times_arr.size()-1)) {
          if (!rank) printf("libROM:   Calculating diff between PDE and ROM solutions.\n");
          for (int ns = 0; ns < nsims; ns++) {
            CalculateROMDiff(  &(sim[ns].solver),
                               &(sim[ns].mpi) );
          }
        }
        /* write the ROM solution to file */
        OutputROMSolution(sim, nsims, waqt); 
  
      }

      if (!rank) {
        printf( "libROM: total prediction/query wallclock time: %f (seconds).\n",
                total_rom_predict_time );
      }

      ((libROMInterface*)context.rom_interface)->saveROM();

    } else {

      for (int ns = 0; ns < nsims; ns++) {
        sim[ns].solver.rom_diff_norms[0]
          = sim[ns].solver.rom_diff_norms[1]
          = sim[ns].solver.rom_diff_norms[2]
          = -1;
      }

    }

  } else if (context.rom_mode == _ROM_MODE_PREDICT_) {

    for (int ns = 0; ns < nsims; ns++) {
      sim[ns].solver.rom_diff_norms[0]
        = sim[ns].solver.rom_diff_norms[1]
        = sim[ns].solver.rom_diff_norms[2]
        = -1;
      strcpy(sim[ns].solver.ConservationCheck,"no");
    }

    ((libROMInterface*)context.rom_interface)->loadROM();
    ((libROMInterface*)context.rom_interface)->projectInitialSolution(sim);

    {
      int start_iter = sim[0].solver.restart_iter;
      int n_iter = sim[0].solver.n_iter;
      double dt = sim[0].solver.dt;
  
      double cur_time = start_iter * dt;
      context.op_times_arr.push_back(cur_time);
  
      for (int iter = start_iter; iter < n_iter; iter++) {
        cur_time += dt;
        if (    ( (iter+1)%sim[0].solver.file_op_iter == 0)
            &&  ( (iter+1) < n_iter) ) {
          context.op_times_arr.push_back(cur_time);
        }
      }
  
      double t_final = n_iter*dt;
      context.op_times_arr.push_back(t_final);
    }

    double total_rom_predict_time = 0;
    for (int iter = 0; iter < context.op_times_arr.size(); iter++) {
  
      double waqt = context.op_times_arr[iter];
  
      ((libROMInterface*)context.rom_interface)->predict(sim, waqt);
      if (!rank) printf(  "libROM: Predicted solution at time %1.4e using ROM, wallclock time: %f.\n", 
                          waqt, ((libROMInterface*)context.rom_interface)->predictWallclockTime() );
      total_rom_predict_time += ((libROMInterface*)context.rom_interface)->predictWallclockTime();
  
      /* write the solution to file */
      for (int ns = 0; ns < nsims; ns++) {
        if (sim[ns].solver.PhysicsOutput) {
          sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
                                        &(sim[ns].mpi),
                                        waqt );
        }
      }
      OutputSolution(sim, nsims, waqt); 
  
    }

    /* calculate error if exact solution has been provided */
    for (int ns = 0; ns < nsims; ns++) {
      CalculateError(&(sim[ns].solver),
                     &(sim[ns].mpi) );
    }

    if (!rank) {
      printf( "libROM: total prediction/query wallclock time: %f (seconds).\n",
              total_rom_predict_time );
    }

  }

  delete ((libROMInterface*)context.rom_interface);
#endif

  PetscFunctionReturn(0);
}

int Solve	(	void *	s,
		int	nsims,
		int	rank,
		int	nproc
	)

Solve the PDE - time-integration

This function integrates the semi-discrete ODE (obtained from discretizing the PDE in space) using natively implemented time integration methods. It initializes the time integration object, iterates the simulation for the required number of time steps, and calculates the errors. After the specified number of iterations, it writes out some information to the screen and the solution to a file.

Parameters

s	Array of simulation objects of type SimulationObject
nsims	number of simulation objects
rank	MPI rank of this process
nproc	Number of MPI processes

Definition at line 37 of file Solve.cpp.

{
  SimulationObject* sim = (SimulationObject*) s;

  /* make sure none of the simulation objects sent in the array 
   * are "barebones" type */
  for (int ns = 0; ns < nsims; ns++) {
    if (sim[ns].is_barebones == 1) {
      fprintf(stderr, "Error in Solve(): simulation object %d on rank %d is barebones!\n",
              ns, rank );
      return 1;
    }
  }

  /* write out iblank to file for visualization */
  for (int ns = 0; ns < nsims; ns++) {
    if (sim[ns].solver.flag_ib) {

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

      WriteArray( sim[ns].solver.ndims,
                  1,
                  sim[ns].solver.dim_global,
                  sim[ns].solver.dim_local,
                  sim[ns].solver.ghosts,
                  sim[ns].solver.x,
                  sim[ns].solver.iblank,
                  &(sim[ns].solver),
                  &(sim[ns].mpi),
                  fname_root );
    }
  }

#ifdef with_librom
  if (!rank) printf("Setting up libROM interface.\n");
  libROMInterface rom_interface( sim, nsims, rank, nproc, sim[0].solver.dt );
  const std::string& rom_mode( rom_interface.mode() );
  std::vector<double> op_times_arr(0);
#endif

#ifdef with_librom
  if ((rom_mode == _ROM_MODE_TRAIN_) || (rom_mode == _ROM_MODE_NONE_)) {
#endif
    /* Define and initialize the time-integration object */
    TimeIntegration TS;
    if (!rank) printf("Setting up time integration.\n");
    TimeInitialize(sim, nsims, rank, nproc, &TS);
    double ti_runtime = 0.0;

    if (!rank) printf("Solving in time (from %d to %d iterations)\n",TS.restart_iter,TS.n_iter);
    for (TS.iter = TS.restart_iter; TS.iter < TS.n_iter; TS.iter++) {
  
      /* Write initial solution to file if this is the first iteration */
      if (!TS.iter) { 
        for (int ns = 0; ns < nsims; ns++) {
          if (sim[ns].solver.PhysicsOutput) {
            sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
                                          &(sim[ns].mpi),
                                          TS.waqt );
          }
        }
        OutputSolution(sim, nsims, TS.waqt); 
#ifdef with_librom
        op_times_arr.push_back(TS.waqt);
#endif
      }
  
#ifdef with_librom
      if ((rom_mode == _ROM_MODE_TRAIN_) && (TS.iter%rom_interface.samplingFrequency() == 0)) {
        rom_interface.takeSample( sim, TS.waqt );
      }
#endif
  
      /* Call pre-step function */
      TimePreStep (&TS);
#ifdef compute_rhs_operators
      /* compute and write (to file) matrix operators representing the right-hand side */
//      if (((TS.iter+1)%solver->file_op_iter == 0) || (!TS.iter)) 
//        { ComputeRHSOperators(solver,mpi,TS.waqt);
#endif
  
      /* Step in time */
      TimeStep (&TS);
  
      /* Call post-step function */
      TimePostStep (&TS);
  
      ti_runtime += TS.iter_wctime;
  
      /* Print information to screen */
      TimePrintStep(&TS);
  
      /* Write intermediate solution to file */
      if (      ((TS.iter+1)%sim[0].solver.file_op_iter == 0) 
            &&  ((TS.iter+1) < TS.n_iter) ) { 
        for (int ns = 0; ns < nsims; ns++) {
          if (sim[ns].solver.PhysicsOutput) {
            sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
                                          &(sim[ns].mpi),
                                          TS.waqt );
          }
        }
        OutputSolution(sim, nsims, TS.waqt); 
#ifdef with_librom
        op_times_arr.push_back(TS.waqt);
#endif
      }
  
    }
  
    double t_final = TS.waqt;
    TimeCleanup(&TS);

    if (!rank) {
      printf( "Completed time integration (Final time: %f), total wctime: %f (seconds).\n",
              t_final, ti_runtime );
      if (nsims > 1) printf("\n");
    }

    /* calculate error if exact solution has been provided */
    for (int ns = 0; ns < nsims; ns++) {
      CalculateError(&(sim[ns].solver),
                     &(sim[ns].mpi) );
    }

    /* write a final solution file */
    for (int ns = 0; ns < nsims; ns++) {
      if (sim[ns].solver.PhysicsOutput) {
        sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
                                      &(sim[ns].mpi),
                                      t_final );
      }
    }
    OutputSolution(sim, nsims, t_final); 

#ifdef with_librom
    op_times_arr.push_back(TS.waqt);

    for (int ns = 0; ns < nsims; ns++) {
      ResetFilenameIndex( sim[ns].solver.filename_index, 
                          sim[ns].solver.index_length );
    }
  
    if (rom_interface.mode() == _ROM_MODE_TRAIN_) {
  
      rom_interface.train();
      if (!rank) printf("libROM: total training wallclock time: %f (seconds).\n", 
                        rom_interface.trainWallclockTime() );

      double total_rom_predict_time = 0;
      for (int iter = 0; iter < op_times_arr.size(); iter++) {
  
        double waqt = op_times_arr[iter];
  
        rom_interface.predict(sim, waqt);
        if (!rank) printf(  "libROM: Predicted solution at time %1.4e using ROM, wallclock time: %f.\n", 
                            waqt, rom_interface.predictWallclockTime() );
        total_rom_predict_time += rom_interface.predictWallclockTime();
  
        /* calculate diff between ROM and PDE solutions */
        if (iter == (op_times_arr.size()-1)) {
          if (!rank) printf("libROM:   Calculating diff between PDE and ROM solutions.\n");
          for (int ns = 0; ns < nsims; ns++) {
            CalculateROMDiff(  &(sim[ns].solver),
                               &(sim[ns].mpi) );
          }
        }
        /* write the ROM solution to file */
        OutputROMSolution(sim, nsims,waqt); 
  
      }

      if (!rank) {
        printf( "libROM: total prediction/query wallclock time: %f (seconds).\n",
                total_rom_predict_time );
      }

      rom_interface.saveROM();

    } else {

      for (int ns = 0; ns < nsims; ns++) {
        sim[ns].solver.rom_diff_norms[0]
          = sim[ns].solver.rom_diff_norms[1]
          = sim[ns].solver.rom_diff_norms[2]
          = -1;
      }

    }

  } else if (rom_mode == _ROM_MODE_PREDICT_) {

    for (int ns = 0; ns < nsims; ns++) {
      sim[ns].solver.rom_diff_norms[0]
        = sim[ns].solver.rom_diff_norms[1]
        = sim[ns].solver.rom_diff_norms[2]
        = -1;
      strcpy(sim[ns].solver.ConservationCheck,"no");
    }

    rom_interface.loadROM();
    rom_interface.projectInitialSolution(sim);

    {
      int start_iter = sim[0].solver.restart_iter;
      int n_iter = sim[0].solver.n_iter;
      double dt = sim[0].solver.dt;
  
      double cur_time = start_iter * dt;
      op_times_arr.push_back(cur_time);
  
      for (int iter = start_iter; iter < n_iter; iter++) {
        cur_time += dt;
        if (    ( (iter+1)%sim[0].solver.file_op_iter == 0)
            &&  ( (iter+1) < n_iter) ) {
          op_times_arr.push_back(cur_time);
        }
      }
  
      double t_final = n_iter*dt;
      op_times_arr.push_back(t_final);
    }

    double total_rom_predict_time = 0;
    for (int iter = 0; iter < op_times_arr.size(); iter++) {
  
      double waqt = op_times_arr[iter];
  
      rom_interface.predict(sim, waqt);
      if (!rank) printf(  "libROM: Predicted solution at time %1.4e using ROM, wallclock time: %f.\n", 
                          waqt, rom_interface.predictWallclockTime() );
      total_rom_predict_time += rom_interface.predictWallclockTime();
  
      /* write the solution to file */
      for (int ns = 0; ns < nsims; ns++) {
        if (sim[ns].solver.PhysicsOutput) {
          sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
                                        &(sim[ns].mpi),
                                        waqt );
        }
      }
      OutputSolution(sim, nsims, waqt); 
  
    }

    /* calculate error if exact solution has been provided */
    for (int ns = 0; ns < nsims; ns++) {
      CalculateError(&(sim[ns].solver),
                     &(sim[ns].mpi) );
    }

    if (!rank) {
      printf( "libROM: total prediction/query wallclock time: %f (seconds).\n",
              total_rom_predict_time );
    }

  }
#endif

  return 0;
}