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"
 
 
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 + 1;                        \
    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[])
{
    printf("\n");
    int k, 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 != 1)
    {
        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 = _1D_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 m = (int)FV0.RHO[0];
  const double h = config[10], gamma = config[6];
  // The number of times steps of the fluid data stored for plotting.
  const int N = 2; // (int)(config[5]) + 1;
  double time_plot[2];
 
  struct cell_var_stru CV = {NULL}; // Structure of fluid variables in computational cells array pointer.
  double ** X = NULL;
  double * cpu_time = 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[5],"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, time_plot);
          break;
          case 2:
          GRP_solver_LAG_source(m, CV, X, 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 if (strcmp(argv[5],"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, time_plot);
          break;
          case 2:
          GRP_solver_EUL_source(m, 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.
  _1D_file_write(m, N, CV, X, cpu_time, argv[2], time_plot);
 
 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;
}