InitialSolution.c File Reference

Read in initial solution from file. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <basic.h>
#include <common.h>
#include <arrayfunctions_gpu.h>
#include <io.h>
#include <mpivars.h>
#include <simulation_object.h>

Go to the source code of this file.

Functions
int	VolumeIntegral (double *, double *, void *, void *)
int	InitialSolution (void *s, int nsims)

Detailed Description

Read in initial solution from file.

Author

Debojyoti Ghosh, Youngdae Kim

Definition in file InitialSolution.c.

Function Documentation

int VolumeIntegral	(	double *	VolumeIntegral,
		double *	u,
		void *	s,
		void *	m
	)

Compute the volume integral of the solution.

Parameters

VolumeIntegral	The computed volume integral
u	Solution
s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 14 of file VolumeIntegral.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;

  int ndims   = solver->ndims;
  int nvars   = solver->nvars;
  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int d,v;

  /* calculate local volume integral of the solution */
  int index[ndims], done = 0;
  double *local_integral = (double*) calloc (nvars,sizeof(double));
  _ArraySetValue_(local_integral,nvars,0.0);
  _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
    double dxinv[ndims];
    for (d=0; d<ndims; d++) { _GetCoordinate_(d,index[d],dim,ghosts,solver->dxinv,dxinv[d]); }
    double dV = 1.0; for (d=0; d<ndims; d++) dV *= (1.0/dxinv[d]);
    for (v=0; v<nvars; v++) local_integral[v] += (u[p*nvars+v]*dV);
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }
  /* sum over all processors to get global integral of the solution */
  IERR MPISum_double(VolumeIntegral,local_integral,nvars,&mpi->world); CHECKERR(ierr);
  free(local_integral);

  return(0);
}

int InitialSolution	(	void *	s,
		int	nsims
	)

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;
}