hydrocode.c

浏览该文件的文档.

 
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
 
#include "../include/var_struc.h"
#include "../include/file_io.h"
#include "../include/finite_volume.h"
 
 
#ifdef DOXYGEN_PREDEFINED
#define NODATPLOT
#define HDF5PLOT
#endif
 
double config[N_CONF]; 
 
#define CV_INIT_MEM(v, N)                       \
    do {                                \
    CV.v = (double **)malloc(N * sizeof(double *));         \
    if(CV.v == NULL)                        \
        {                               \
        printf("NOT enough memory! %s\n", #v);          \
        retval = 5;                     \
        goto return_NULL;                   \
        }                               \
    CV.v[0] = FV0.v;                        \
    for(k = 1; k < N; ++k)                      \
        {                               \
        CV.v[k] = (double *)malloc(m * sizeof(double));     \
        if(CV.v[k] == NULL)                 \
            {                           \
            printf("NOT enough memory! %s[%d]\n", #v, k);   \
            retval = 5;                 \
            goto return_NULL;               \
            }                           \
        }                               \
    } while (0)
 
int main(int argc, char *argv[])
{
  int k, j, retval = 0;
  // Initialize configuration data array
  for(k = 1; k < N_CONF; k++)
      config[k] = INFINITY;
 
  char * scheme = NULL; // Riemann_exact(Godunov), GRP
  arg_preprocess(4, argc, argv, scheme);
 
  // Set dimension.
  config[0] = (double)1; // Dimensionality = 1
 
  // The number of times steps of the fluid data stored for plotting.
  int N, N_plot;
  double * time_plot;
    /* 
     * We read the initial data files.
     * The function initialize return a point pointing to the position
     * of a block of memory consisting (m) variables of type double.
     * The (m) array elements of these variables are the initial value.
     */
  struct flu_var FV0 = initialize_1D(argv[1], &N, &N_plot, &time_plot); // Structure of initial data array pointer.
    /* 
     * (m) is the number of initial value as well as the number of grids.
     * As (m) is frequently use to represent the number of grids,
     * we do not use the name such as num_cell here to correspond to
     * notation in the math theory.
     */
  const int m = (int)config[3];
  const double h = config[10], gamma = config[6];
  const int order = (int)config[9];
 
  // Structure of fluid variables in computational cells array pointer.
  struct cell_var_stru CV = {NULL};
  double ** X = NULL;
  double * cpu_time = (double *)malloc(N * sizeof(double));
  X = (double **)malloc(N * sizeof(double *));
  if(cpu_time == NULL)
      {
      printf("NOT enough memory! CPU_time\n");
      retval = 5;
      goto return_NULL;
      }
  if(X == NULL)
      {
      printf("NOT enough memory! X\n");
      retval = 5;
      goto return_NULL;
      }
  for(k = 0; k < N; ++k)
  {
    X[k] = (double *)malloc((m+1) * sizeof(double));
    if(X[k] == NULL)
    {
      printf("NOT enough memory! X[%d]\n", k);
      retval = 5;
      goto return_NULL;
    }
  }
  // Initialize arrays of fluid variables in cells.
  CV_INIT_MEM(RHO, N);
  CV_INIT_MEM(U, N);
  CV_INIT_MEM(P, N);
  CV.E = (double **)malloc(N * sizeof(double *));
  if(CV.E == NULL)
      {
      printf("NOT enough memory! E\n");
      retval = 5;
      goto return_NULL;
      }
  for(k = 0; k < N; ++k)
  {
    CV.E[k] = (double *)malloc(m * sizeof(double));
    if(CV.E[k] == NULL)
    {
      printf("NOT enough memory! E[%d]\n", k);
      retval = 5;
      goto return_NULL;
    }
  }
  // Initialize the values of energy in computational cells and x-coordinate of the cell interfaces.
  for(j = 0; j <= m; ++j)
      X[0][j] = h * j;
  for(j = 0; j < m; ++j)
      CV.E[0][j] = 0.5*CV.U[0][j]*CV.U[0][j] + CV.P[0][j]/(gamma - 1.0)/CV.RHO[0][j];
 
  if (strcmp(argv[4],"LAG") == 0) // Use GRP/Godunov scheme to solve it on Lagrangian coordinate.
      {
      config[8] = (double)1;
      switch(order)
          {
          case 1:
          Godunov_solver_LAG_source(m, CV, X, cpu_time, &N, time_plot);
          break;
          case 2:
          GRP_solver_LAG_source(m, CV, X, cpu_time, &N, time_plot);
          break;
          default:
          printf("NOT appropriate order of the scheme! The order is %d.\n", order);
          retval = 4;
          goto return_NULL;
          }
      }
  else if (strcmp(argv[4],"EUL") == 0) // Use GRP/Godunov scheme to solve it on Eulerian coordinate.
      {
      config[8] = (double)0;
      for (k = 1; k < N; ++k)
          for (j = 0; j <= m; ++j)
          X[k][j] = X[0][j];
      switch(order)
          {
          case 1:
          Godunov_solver_EUL_source(m, CV, cpu_time, &N, time_plot);
          break;
          case 2:
          GRP_solver_EUL_source(m, CV, cpu_time, &N, time_plot);
          break;
          default:
          printf("NOT appropriate order of the scheme! The order is %d.\n", order);
          retval = 4;
          goto return_NULL;
          }
      }
  else
      {
      printf("NOT appropriate coordinate framework! The framework is %s.\n", argv[4]);
      retval = 4;
      goto return_NULL;
      }
 
  // Write the final data down.
#ifndef NODATPLOT
  file_1D_write(m, N, CV, X, cpu_time, argv[2], time_plot);
#endif
#ifdef HDF5PLOT
  file_1D_write_HDF5(m, N, CV, X, cpu_time, argv[2], time_plot);
#endif
 
 return_NULL:
  free(FV0.RHO);
  free(FV0.U);
  free(FV0.P);
  FV0.RHO = NULL;
  FV0.U   = NULL;
  FV0.P   = NULL;
  for(k = 1; k < N; ++k)
  {
    free(CV.E[k]);
    free(CV.RHO[k]);
    free(CV.U[k]);
    free(CV.P[k]);
    free(X[k]);
    CV.E[k]   = NULL;
    CV.RHO[k] = NULL;
    CV.U[k]   = NULL;
    CV.P[k]   = NULL;
    X[k] = NULL;
  }
  free(CV.E[0]);
  CV.E[0]   = NULL;
  CV.RHO[0] = NULL;
  CV.U[0]   = NULL;
  CV.P[0]   = NULL;
  free(CV.E);
  free(CV.RHO);
  free(CV.U);
  free(CV.P);
  CV.E   = NULL;
  CV.RHO = NULL;
  CV.U   = NULL;
  CV.P   = NULL;
  free(X);
  X = NULL;
  free(cpu_time);
  cpu_time = NULL;
  
  return retval;
}