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
#endif
#ifdef DOXYGEN_PREDEFINED
#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[])
{
    printf("\n");
    int k, i, j, retval = 0;
    for (k = 0; k < argc; k++)
    printf("%s ", argv[k]);
    printf("\n");
    printf("TEST:\n  %s\n", argv[1]);
    if(argc < 5)
    {
        printf("Test Beginning: ARGuments Counter %d is less than 5.\n", argc);
        return 4;
    }
    else
    printf("Test Beginning: ARGuments Counter = %d.\n", argc);
 
    // Initialize configuration data array
    for(k = 1; k < N_CONF; k++)
        config[k] = INFINITY;
 
    // Set dimension.
    int dim;
    dim = atoi(argv[3]);
    if (dim != 2)
    {
        printf("No appropriate dimension was entered!\n");
        return 4;
    }
    config[0] = (double)dim;
 
    printf("Configurating:\n");
    char * endptr;
    double conf_tmp;
    for (k = 6; k < argc; k++)
    {
        errno = 0;
        j = strtoul(argv[k], &endptr, 10);
        if (errno != ERANGE && *endptr == '=')
        {                           
            endptr++;
            errno = 0;
            conf_tmp = strtod(endptr, &endptr);
            if (errno != ERANGE && *endptr == '\0')
            {
                config[j] = conf_tmp;
                printf("%3d-th configuration: %g (ARGument)\n", j, conf_tmp);
            }
            else
            {
                printf("Configuration error in ARGument variable %d! ERROR after '='!\n", k);
                return 4;
            }
        }
        else
        {
            printf("Configuration error in ARGument variable %d! ERROR before '='!\n", k);
            return 4;
        }
    }
    
  // Set order and scheme.
  int order; // 1, 2
  char * scheme; // Riemann_exact(Godunov), GRP
  printf("Order[_Scheme]: %s\n",argv[4]);
  errno = 0;
  order = strtoul(argv[4], &scheme, 10);
  if (*scheme == '_')
      scheme++;
  else if (*scheme != '\0' || errno == ERANGE)
      {
      printf("No order or Wrog scheme!\n");
      return 4;
      }
  config[9] = (double)order;
 
    /* 
     * We read the initial data files.
     * The function initialize return a point pointing to the position
     * of a block of memory consisting (m+1) variables of type double.
     * The value of first array element of these variables is m.
     * The following m variables are the initial value.
     */
    struct flu_var FV0 = _2D_initialize(argv[1]); // 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_grid here to correspond to
     * notation in the math theory.
     */
  const int n_x = (int)FV0.RHO[1], n_y = (int)FV0.RHO[0];
  const double h_x = config[10], h_y = config[11], gamma = config[6];
  // The number of times steps of the fluid data stored for plotting.
  int N = 2; // (int)(config[5]) + 1;
  double time_plot[2];
 
  // Structure of fluid variables in computational cells array pointer.
  struct cell_var_stru * CV = malloc(N * sizeof(struct cell_var_stru));
  double ** X, ** Y;
  double * cpu_time = 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 + 2];
          CV[0].U[j][i]   = FV0.U[i*n_x + j + 2];
          CV[0].V[j][i]   = FV0.V[i*n_x + j + 2];
          CV[0].P[j][i]   = FV0.P[i*n_x + j + 2];
          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];
      }
 
  _Bool const dim_split = (_Bool)config[33]; // Dimensional splitting?
  if (strcmp(argv[5],"EUL") == 0) // Use GRP/Godunov scheme to solve it on Eulerian coordinate.
      {
      config[8] = (double)0;
      switch(order)
          {
          case 1:
          // Godunov_solver_2D_EUL_source(n_x, n_y, CV, cpu_time);
          config[41] = 0.0; // alpha = 0.0
          GRP_solver_2D_EUL_source(n_x, n_y, CV, cpu_time, time_plot);
          break;
          case 2:
          if (dim_split)
              GRP_solver_2D_split_EUL_source(n_x, n_y, CV, cpu_time, time_plot);
          else
              GRP_solver_2D_EUL_source(n_x, n_y, CV, cpu_time, 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[5]);
      retval = 4;
      goto return_NULL;
      }
 
  // Write the final data down.
#ifndef NODATPLOT
  _2D_file_write(n_x, n_y, N, CV, X, Y, cpu_time, argv[2], time_plot);
#endif
#ifndef NOTECPLOT
  _2D_TEC_file_write(n_x, n_y, N, CV, X, Y, cpu_time, argv[2], time_plot);
#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;
}