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
#define NOTECPLOT
#endif
 
double config[N_CONF]; 
 
#define CV_INIT_MEM(v, N)                       \
    do {                                \
    for(k = 0; k < N; ++k)                      \
    {                               \
        CV[k].v = (double **)malloc(n_x * sizeof(double *));    \
        if(CV[k].v == NULL)                     \
        {                           \
            printf("NOT enough memory! CV[%d].%s\n", k, #v);    \
            retval = 5;                     \
            goto return_NULL;                   \
        }                           \
        for(j = 0; j < n_x; ++j)                    \
        {                           \
            CV[k].v[j] = (double *)malloc(n_y * sizeof(double)); \
            if(CV[k].v[j] == NULL)              \
            {                       \
                printf("NOT enough memory! CV[%d].%s[%d]\n", k, #v, j); \
                retval = 5;                 \
                goto return_NULL;               \
            }                       \
        }                           \
    }                               \
    } while (0)
 
int main(int argc, char *argv[])
{
  int k, i, 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)2; // Dimensionality = 2
 
  // 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 (n_x*n_y) variables of type double.
     * The (n_x*n_y) array elements of these variables are the initial value.
     */
  struct flu_var FV0 = initialize_2D(argv[1], &N, &N_plot, &time_plot); // Structure of initial data array pointer.
    /* 
     * (n_x*n_y) is the number of initial value as well as the number of grids.
     * As (n_x*n_y) 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 n_x = (int)config[13], n_y = (int)config[14];
  const double h_x = config[10], h_y = config[11], gamma = config[6];
  const int order = (int)config[9];
  const _Bool dim_split = (_Bool)config[33]; // Dimensional splitting?
 
  // Structure of fluid variables in computational cells array pointer.
  struct cell_var_stru * CV = (struct cell_var_stru*)malloc(N * sizeof(struct cell_var_stru));
  double ** X, ** Y;
  double * cpu_time = (double *)malloc(N * sizeof(double));
  X = (double **)malloc((n_x+1) * sizeof(double *));
  Y = (double **)malloc((n_x+1) * sizeof(double *));
  if(cpu_time == NULL)
      {
      printf("NOT enough memory! CPU_time\n");
      retval = 5;
      goto return_NULL;
      }
  if(X == NULL || Y == NULL)
      {
      printf("NOT enough memory! X or Y\n");
      retval = 5;
      goto return_NULL;
      }
  for(j = 0; j <= n_x; ++j)
  {
    X[j] = (double *)malloc((n_y+1) * sizeof(double));
    Y[j] = (double *)malloc((n_y+1) * sizeof(double));
    if(X[j] == NULL || Y[j] == NULL)
    {
      printf("NOT enough memory! X[%d] or Y[%d]\n", j, j);
      retval = 5;
      goto return_NULL;
    }
  }
  if(CV == NULL)
      {
      printf("NOT enough memory! Cell Variables\n");
      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(V, N);
  CV_INIT_MEM(P, N);
  CV_INIT_MEM(E, N);
  // Initialize the values of energy in computational cells and (x,y)-coordinate of the cell interfaces.
  for(j = 0; j <= n_x; ++j)
      for(i = 0; i <= n_y; ++i) 
      {
          X[j][i] = j * h_x;
          Y[j][i] = i * h_y;
      }
  for(j = 0; j < n_x; ++j)
      for(i = 0; i < n_y; ++i)  
      {
          CV[0].RHO[j][i] = FV0.RHO[i*n_x + j];
          CV[0].U[j][i]   = FV0.U[i*n_x + j];
          CV[0].V[j][i]   = FV0.V[i*n_x + j];
          CV[0].P[j][i]   = FV0.P[i*n_x + j];
          CV[0].E[j][i]   = 0.5*CV[0].U[j][i]*CV[0].U[j][i] + CV[0].P[j][i]/(gamma - 1.0)/CV[0].RHO[j][i];
          CV[0].E[j][i]  += 0.5*CV[0].V[j][i]*CV[0].V[j][i];
      }
 
  if (strcmp(argv[4],"EUL") == 0) // Use GRP/Godunov scheme to solve it on Eulerian coordinate.
      {
      config[8] = (double)0;
      switch(order)
          {
          case 1:
          config[41] = 0.0; // alpha = 0.0
          case 2:
          if (dim_split)
              GRP_solver_2D_split_EUL_source(n_x, n_y, CV, X, Y, cpu_time, argv[2], N, &N_plot, time_plot);
          else
              GRP_solver_2D_EUL_source(n_x, n_y, CV, X, Y, cpu_time, argv[2], N, &N_plot, 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_2D_write(n_x, n_y, N_plot, CV, X, Y, cpu_time, argv[2], time_plot);
#endif
#ifdef HDF5PLOT
  file_2D_write_HDF5(n_x, n_y, N_plot, CV, X, Y, cpu_time, argv[2], time_plot);
#endif
#ifndef NOTECPLOT
  file_2D_write_POINT_TEC(n_x, n_y, 1, CV + N_plot-1, X, Y, cpu_time, argv[2], time_plot + N_plot-1);
#endif
 
 return_NULL:
  free(FV0.RHO);
  free(FV0.U);
  free(FV0.V);
  free(FV0.P);
  FV0.RHO = NULL;
  FV0.U   = NULL;
  FV0.V   = NULL;
  FV0.P   = NULL;
  for(k = 0; k < N; ++k)
  {
    for(j = 0; j < n_x; ++j)
    {
            free(CV[k].RHO[j]);
            free(CV[k].U[j]);
            free(CV[k].V[j]);
            free(CV[k].P[j]);
            free(CV[k].E[j]);
            CV[k].RHO[j] = NULL;
            CV[k].U[j]   = NULL;
            CV[k].V[j]   = NULL;
            CV[k].P[j]   = NULL;
            CV[k].E[j]   = NULL;
    }
    free(CV[k].RHO);
    free(CV[k].U);
    free(CV[k].V);
    free(CV[k].P);
    free(CV[k].E);
    CV[k].RHO = NULL;
    CV[k].U   = NULL;
    CV[k].V   = NULL;
    CV[k].P   = NULL;
    CV[k].E   = NULL;
  }
  free(CV);
  CV = NULL;
  for(j = 0; j <= n_x; ++j)
  {
      free(X[j]);
      free(Y[j]);
      X[j] = NULL;
      Y[j] = NULL;
  }
  free(X);
  free(Y);
  X = NULL;
  Y = NULL; 
  free(cpu_time);
  cpu_time = NULL;
  
  return retval;
}