sph.h

用SPH方法编写的激波管程序

#include"iostream.h"
#include"fstream.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include"iomanip.h"
class SPH
{
public:
//-----------------param-------------------------------------------------------------------------
      int dim;//dim : Dimension of the problem (1, 2 or 3)
      int maxn;
	  int max_interaction;//maxn    : Maximum number of particles;max_interation : Maximum number of interaction pairs

      double x_maxgeom,x_mingeom,y_maxgeom,y_mingeom,z_maxgeom,z_mingeom;//Parameters for the computational geometry,x_maxgeom : Upper limit of allowed x-regime 
                                                                                   //x_mingeom : Lower limit of allowed x-regime,y_maxgeom : Upper limit of allowed y-regime 
                                                                                   //y_mingeom : Lower limit of allowed y-regime,z_maxgeom : Upper limit of allowed z-regime 
                                                                                   //z_mingeom : Lower limit of allowed z-regime 

      int pa_sph; //SPH algorithm for particle approximation (pa_sph);
                  //pa_sph = 1 : (e.g. (p(i)+p(j))/(rho(i)*rho(j))
                  //         2 : (e.g. (p(i)/rho(i)**2+p(j)/rho(j)**2)
      int nnps;//     Nearest neighbor particle searching (nnps) method
               //     nnps = 1 : Simplest and direct searching
               //            2 : Sorting grid linked list
               //            3 : Tree algorithm 

      int sle; //     Smoothing length evolution (sle) algorithm
               //     sle = 0 : Keep unchanged,
               //           1 : h = fac * (m/rho)^(1/dim)
               //           2 : dh/dt = (-1/dim)*(h/rho)*(drho/dt)
               //           3 : Other approaches (e.g. h = h_0 * (rho_0/rho)**(1/dim) ) 
      int skf; //     Smoothing kernel function 
               //     skf = 1, cubic spline kernel by W4 - Spline (Monaghan 1985)
               //         = 2, Gauss kernel   (Gingold and Monaghan 1981) 
               //         = 3, Quintic kernel (Morris 1997)

      bool summation_density, average_velocity, config_input,virtual_part, vp_input, visc;
	  bool ex_force, heat_artificial,visc_artificial, self_gravity, nor_density;
	  //     Switches for different senarios
      //     summation_density = .TRUE. : Use density summation model in the code, 
      //                        .FALSE.: Use continuiity equation
      //     average_velocity = .TRUE. : Monaghan treatment on average velocity,
      //                       .FALSE.: No average treatment.
      //     config_input = .TRUE. : Load initial configuration data,
      //                   .FALSE.: Generate initial configuration.
      //     virtual_part = .TRUE. : Use vritual particle,
      //                   .FALSE.: No use of vritual particle.
      //     vp_input = .TRUE. : Load virtual particle information,
      //               .FALSE.: Generate virtual particle information.
      //     visc = .true. : Consider viscosity,
      //           .false.: No viscosity.
      //     ex_force =.true. : Consider external force,
      //               .false.: No external force.
      //     visc_artificial = .true. : Consider artificial viscosity,
      //                      .false.: No considering of artificial viscosity.
      //     heat_artificial = .true. : Consider artificial heating,
      //                      .false.: No considering of artificial heating.
      //     self_gravity = .true. : Considering self_gravity,
      //                    .false.: No considering of self_gravity
      //     nor_density =  .true. : Density normalization by using CSPM,
      //                    .false.: No normalization.
   
      int    nsym;//     Symmetry of the problem
                  //     nsym = 0 : no symmetry,
                  //          = 1 : axis symmetry,
                  //          = 2 : center symmetry.  
//     Control parameters for output 
      bool int_stat;//     int_stat = .true. : Print statistics about SPH particle interactions.
                    //                        including virtual particle information
      int print_step, save_step, moni_particle;//     print_step: Print Timestep (On Screen)
                                               //     save_step : Save Timestep    (To Disk File)
                                               //     moni_particle: The particle number for information monitoring.
       
      double  pi;//changed

	  //     Simulation cases
      bool shocktube, shearcavity;//     shocktube = .true. : carry out shock tube simulation
                                  //     shearcavity = .true. : carry out shear cavity simulation
//------------------------------------------------------------------------------------------------------
//--------------------------------general variables------------------------------
 
//-------------------------------------------------------------------------------




//????????????????????以下所有变量均没有在解析函数中定义初始值。
//--------------------------------------input-----------------------------------
double *p;
//-------------------------------------------------------------------------------

	  
	  
//--------------------------------------kernel-----------------------------------
//double r;
//double *dw1dx;//程序中出现一个单参数和一个双参数的dwdx，为了区别特改为dw1dx
double *dx;
//double w1;//程序中出现一个整型和一个单参数的w，为了区别特改为w1
//-------------------------------------------------------------------------------
//-----------------------------------h_upgrade--------------------------------------
double dt;//???????????/是局部变量还是全局变量？？？源程序中是放在局部变量的位置
//------------------------------------------------------------------------------

//----------------------------------------direct_find----------------------------
int itimestep;
int *countiac;
//-------------------------------------------------------------------------------	  
//---------------------------------art_heat------------------------------------------------------

   int ntotal,niac;
   int *pair_i;
   int *pair_j;
   double *hsml;
   double *mass;
   double **x;
   double **vx;
   double *rho;
   double *u;
   double *c;//在人工热量中表示声速，在人工粘性中表示温度，因此需要区分？？？？？？
   double *w;
   double **dwdx;
//   double *dedt;	  
//--------------------------------------------------------------------------------------------------
//----------------------------------art_vsic--------------------------------
//   double **dvxdt;
//--------------------------------------------------------------------------
//----------------------------------av_vel-----------------------------------
  double **av;
//---------------------------------------------------------------------------------
//----------------------------------sum_densty-----------------------------
  int *itype;
//----------------------------------------------------------------------------
//----------------------------------con_densty-----------------------------
  double *drhodt;
//----------------------------------------------------------------------------
//-----------------------------------ext_force--------------------------------------

//------------------------------------------------------------------------------
//-----------------------------------int_force--------------------------------------
//此程序中出现两个itype的定义，请查明原因
double *eta;//好像程序的另一个地方也定义了这个变量？？？？？？？？？
double *t;
double *tdsdt;
    
//------------------------------------------------------------------------------
//-----------------------------------single_step--------------------------------------
double *s;//此变量本程序中没有使用
double *du;
double *ds;
double *drho;
double **dvx;
//------------------------------------------------------------------------------
//-----------------------------------time_integration--------------------------------------
//double *e;
int maxtimestep;
//------------------------------------------------------------------------------
//-----------------------------------virt_part--------------------------------------
int nvirt;
//------------------------------------------------------------------------------
//-----------------------------------viscosity--------------------------------------

//------------------------------------------------------------------------------
//p_gas   p_art_water临时设置的变量，以后应拿掉
//double p1;
//double c1;

 
public:
	SPH(int,int);
	virtual ~SPH();
public:

    void direct_find(int);
	void kernel(double r,double *dx,double ,double &w1,double *(&dw1dx));    
    void art_visc(int,double **(&dvx1dt),double *(&de1dt));
   
    void sum_density(int);
 
	
    void int_force(int ntotal,double **(&dvxdt),double *tdsdt,double *(&dedt));
    void single_step();
    void time_integration();

	void input();
    void output();                                                
   void friend main();

};