<?xml version="1.0" encoding="utf-8"?>
4 问题
HiFiLES
Table of Contents
- 1. 类
- 2. 函数
- 2.1. CalcResidual(FlowSol.iniiter+isteps, i, &FlowSol)
- 2.1.1. DONE eles::extrapolatesolution(int indisuuptsfrom)
- 2.1.2. DONE void eles::evaluateinvFlux(int indisuuptsfrom)
- 2.1.3. DONE void intinters::calculatecommoninvFlux(void)
- 2.1.4. DONE viscous
- 2.1.5. DONE extrapolatetotalFlux()
- 2.1.6. DONE calculatedivergence(indivtconfuptsto)
- 2.1.7. HOLD viscous
- 2.1.8. DONE calculatecorrecteddivergence(indivtconfuptsto)
- 2.1. CalcResidual(FlowSol.iniiter+isteps, i, &FlowSol)
- 3. 程序变动笔记
- 4. 问题
[0/3]
1 类
1.1 Array
1.1.1 members
public: // #### constructors #### // default constructor array(); // constructor 1 array(int in_dim_0, int in_dim_1=1, int in_dim_2=1, int in_dim_3=1); // copy constructor array(const array<T>& in_array); // assignment array<T>& operator=(const array<T>& in_array); ~array(); protected: int dim_0; int dim_1; int dim_2; int dim_3; T* cpu_data; T* gpu_data; int cpu_flag; int gpu_flag;
1.1.2 methods
// setup void setup(int in_dim_0, int in_dim_1=1, int in_dim_2=1, int in_dim_3=1); // access/set 1d T& operator() (int in_pos_0); //reload operator () // access/set 2d T& operator() (int in_pos_0, int in_pos_1); // access/set 3d T& operator() (int in_pos_0, int in_pos_1, int in_pos_2); // access/set 4d T& operator() (int in_pos_0, int in_pos_1, int in_pos_2, int in_pos_3); // return pointer T* get_ptr_cpu(void); T* get_ptr_gpu(void); // return pointer T* get_ptr_cpu(int in_pos_0, int in_pos_1=0, int in_pos_2=0, int in_pos_3=0); T* get_ptr_gpu(int in_pos_0, int in_pos_1=0, int in_pos_2=0, int in_pos_3=0); // return dimension int get_dim(int in_dim); // method to get maximum value of array T get_max(void); // method to get minimum value of array T get_min(void); // print void print(void); // move data from cpu to gpu void mv_cpu_gpu(void); // copy data from cpu to gpu void cp_cpu_gpu(void); // move data from gpu to cpu void mv_gpu_cpu(void); // copy data from gpu to cpu void cp_gpu_cpu(void); // remove data from cpu void rm_cpu(void); /*! Initialize array to zero - Valid for numeric data types (int, float, double) */ void initialize_to_zero(); /*! Initialize array to given value */ void initialize_to_value(const T val);
1.2 Mesh
1.2.1 members
/** Basic parameters of mesh */ //unsigned long int n_eles, n_verts, n_dims, n_verts_global, n_cells_global; int iter; /** arrays which define the basic mesh geometry */ array<double> xv_0;//, xv; array< array<double> > xv; array<int> c2v,c2n_v,ctype,bctype_c,ic2icg,iv2ivg,ic2loc_c, f2c,f2loc_f,c2f,c2e,f2v,f2n_v,e2v,v2n_e; array<array<int> > v2e; /** #### Boundary information #### */ int n_bnds, n_faces; array<int> nBndPts; array<int> v2bc; /** vertex id = boundpts(bc_id)(ivert) */ array<array<int> > boundPts; /** Store motion flag for each boundary (currently 0=fixed, 1=moving, -1=volume) */ array<int> bound_flags; /** HiFiLES 'bcflag' for each boundary */ array<int> bc_list; /** replacing get_bc_name() from geometry.cpp */ map<string,int> bc_name; /** inverse of bc_name */ map<int,string> bc_flag; // nBndPts.setup(n_bnds); boundPts.setup(nBnds,nPtsPerBnd); array<double> vel_old,vel_new, xv_new; array< array<double> > grid_vel;
1.2.2 methods
public: // #### constructors #### /** default constructor */ mesh(void); /** default destructor */ ~mesh(void); // #### methods #### /** Mesh motion wrapper */ void move(int _iter, int in_rk_step, solution *FlowSol); /** peform prescribed mesh motion using linear elasticity method */ void deform(solution* FlowSol); /** peform prescribed mesh motion using rigid translation/rotation */ void rigid_move(solution *FlowSol); /** Perturb the mesh points (test case for freestream preservation) */ void perturb(solution *FlowSol); /** update grid velocity & apply to eles */ void set_grid_velocity(solution *FlowSol, double dt); /** update the mesh: re-set spts, transforms, etc. */ void update(solution *FlowSol); /** setup information for boundary motion */ //void setup_boundaries(array<int> bctype); /** write out mesh to file */ void write_mesh(int mesh_type, double sim_time); /** write out mesh in Gambit .neu format */ void write_mesh_gambit(double sim_time); /** write out mesh in Gmsh .msh format */ void write_mesh_gmsh(double sim_time); void setup(solution *in_FlowSol, array<double> &in_xv, array<int> &in_c2v, array<int> &in_c2n_v, array<int> &in_iv2ivg, array<int> &in_ctype); private: bool start; array<double> xv_nm1, xv_nm2, xv_nm3;//, xv_new, vel_old, vel_new; /** Global stiffness matrix for linear elasticity solution */ CSysMatrix StiffnessMatrix; CSysVector LinSysRes; CSysVector LinSysSol; /// Copy of the pointer to the Flow Solution structure struct solution *FlowSol; /** global stiffness psuedo-matrix for linear-elasticity mesh motion */ array<array<double> > stiff_mat; unsigned long LinSolIters; int failedIts; double min_vol, min_length, solver_tolerance; double time, rk_time; int rk_step; // Coefficients for LS-RK45 time-stepping array<double> RK_a, RK_b, RK_c; /** create individual-element stiffness matrix - triangles */ bool set_2D_StiffMat_ele_tri(array<double> &stiffMat_ele,int ele_id); /** create individual-element stiffness matrix - quadrilaterals */ bool set_2D_StiffMat_ele_quad(array<double> &stiffMat_ele,int ele_id); /** create individual-element stiffness matrix - tetrahedrons */ //bool set_2D_StiffMat_ele_tet(array<double> &stiffMat_ele,int ele_id, solution *FlowSol); /** create individual-element stiffness matrix - hexahedrons */ //bool set_2D_StiffMat_ele_hex(array<double> &stiffMat_ele,int ele_id, solution *FlowSol); /** * transfrom single-element stiffness matrix to nodal contributions in order to * add to global stiffness matrix */ void add_StiffMat_EleTri(array<double> StiffMatrix_Elem, int id_pt_0, int id_pt_1, int id_pt_2); void add_StiffMat_EleQuad(array<double> StiffMatrix_Elem, int id_pt_0, int id_pt_1, int id_pt_2, int id_pt_3); /** Set given/known displacements of vertices on moving boundaries in linear system */ void set_boundary_displacements(solution *FlowSol); /** meant to check for any inverted cells (I think) and return minimum volume */ double check_grid(solution *FlowSol); /** transfer solution from LinSysSol to xv_new */ void update_grid_coords(void); /** find minimum length in mesh */ void set_min_length(void); /* --- Linear-Elasticy Mesh Deformation Routines Taken from SU^2, 3/26/2014 --- * This version here has been stripped down to the bare essentials needed for HiFiLES * For more details (and additional awesome features), see su2.stanford.edu */ /* These functions will (eventually) be replaced with pre-existing HiFiLES functions, but * for now, a HUGE THANKS to the SU^2 dev team for making this practically plug & play! */ /*! * \brief Add the stiffness matrix for a 2-D triangular element to the global stiffness matrix for the entire mesh (node-based). * \param[in] StiffMatrix_Elem - Element stiffness matrix to be filled. * \param[in] PointCornders - Element vertex ID's */ void add_FEA_stiffMat(array<double> &stiffMat_ele, array<int> &PointCorners); /*! * \brief Build the stiffness matrix for a 3-D hexahedron element. The result will be placed in StiffMatrix_Elem. * \param[in] geometry - Geometrical definition of the problem. * \param[in] StiffMatrix_Elem - Element stiffness matrix to be filled. * \param[in] CoordCorners[8][3] - Index value for Node 1 of the current hexahedron. */ void set_stiffmat_ele_3d(array<double> &stiffMat_ele, int ic, double scale); /*! * \brief Build the stiffness matrix for a 3-D hexahedron element. The result will be placed in StiffMatrix_Elem. * \param[in] geometry - Geometrical definition of the problem. * \param[in] StiffMatrix_Elem - Element stiffness matrix to be filled. * \param[in] CoordCorners[8][3] - Index value for Node 1 of the current hexahedron. */ void set_stiffmat_ele_2d(array<double> &stiffMat_ele, int ic, double scale); /*! * \brief Shape functions and derivative of the shape functions * \param[in] Xi - Local coordinates. * \param[in] Eta - Local coordinates. * \param[in] Mu - Local coordinates. * \param[in] CoordCorners[8][3] - Coordiantes of the corners. * \param[in] shp[8][4] - Shape function information */ double ShapeFunc_Hexa(double Xi, double Eta, double Mu, double CoordCorners[8][3], double DShapeFunction[8][4]); /*! * \brief Shape functions and derivative of the shape functions * \param[in] Xi - Local coordinates. * \param[in] Eta - Local coordinates. * \param[in] Mu - Local coordinates. * \param[in] CoordCorners[8][3] - Coordiantes of the corners. * \param[in] shp[8][4] - Shape function information */ double ShapeFunc_Tetra(double Xi, double Eta, double Mu, double CoordCorners[8][3], double DShapeFunction[8][4]); /*! * \brief Shape functions and derivative of the shape functions * \param[in] Xi - Local coordinates. * \param[in] Eta - Local coordinates. * \param[in] Mu - Local coordinates. * \param[in] CoordCorners[8][3] - Coordiantes of the corners. * \param[in] shp[8][4] - Shape function information */ double ShapeFunc_Pyram(double Xi, double Eta, double Mu, double CoordCorners[8][3], double DShapeFunction[8][4]); /*! * \brief Shape functions and derivative of the shape functions * \param[in] Xi - Local coordinates. * \param[in] Eta - Local coordinates. * \param[in] Mu - Local coordinates. * \param[in] CoordCorners[8][3] - Coordiantes of the corners. * \param[in] shp[8][4] - Shape function information */ double ShapeFunc_Wedge(double Xi, double Eta, double Mu, double CoordCorners[8][3], double DShapeFunction[8][4]); /*! * \brief Shape functions and derivative of the shape functions * \param[in] Xi - Local coordinates. * \param[in] Eta - Local coordinates. * \param[in] Mu - Local coordinates. * \param[in] CoordCorners[8][3] - Coordiantes of the corners. * \param[in] shp[8][4] - Shape function information */ double ShapeFunc_Triangle(double Xi, double Eta, double CoordCorners[8][3], double DShapeFunction[8][4]); /*! * \brief Shape functions and derivative of the shape functions * \param[in] Xi - Local coordinates. * \param[in] Eta - Local coordinates. * \param[in] Mu - Local coordinates. * \param[in] CoordCorners[8][3] - Coordiantes of the corners. * \param[in] shp[8][4] - Shape function information */ double ShapeFunc_Rectangle(double Xi, double Eta, double CoordCorners[8][3], double DShapeFunction[8][4]); ;
1.3 solution
#pragma once //这是一个预编译命令，因为这个头文件要被很多源文件include，所以只加载一次 #include "array.h" //array 主要定义了数组和向量，因为要用到各个各样，不同长度和维度的向量，这里提供了一个模板 #include <string> #include "input.h" //定义了input类，这个主要用于输入，就是获得各种配置参量和网格文件 #include "eles.h" //这个类主要是用于计算Res项，它中定义了各种计算的函数。 #include "eles_tris.h" //这个类是eles的派生类，主要处理三角形 #include "eles_quads.h" //这个类是eles的派生类，主要处理四边形 #include "eles_hexas.h" //这个类是eles的派生类，主要处理六面体 #include "eles_tets.h" //这个类是eles的派生类，主要处理四面体 #include "eles_pris.h" //这个类是eles的派生类，主要处理棱柱 #include "int_inters.h" //int_inters类是inters类的一个派生类，主要处理interface的东西 #include "bdy_inters.h" //bdy_inters类是inters类的一个派生类，主要处理boundary的东西 #ifdef _MPI #include "mpi.h" #include "mpi_inters.h" #endif class int_inters; /*!< Forwards declaration */ class bdy_inters; /*!< Forwards declaration */ class mpi_inters; /*!< Forwards declaration */ struct solution { int viscous; double time; double ene_hist; double grad_ene_hist; array<int> num_f_per_c; int n_ele_types; int n_dims; int num_eles; int num_verts; int num_edges; int num_inters; int num_cells_global; int n_steps; int adv_type; int plot_freq; int restart_dump_freq; int ini_iter; int write_type; array<eles*> mesh_eles; eles_quads mesh_eles_quads; eles_tris mesh_eles_tris; eles_hexas mesh_eles_hexas; eles_tets mesh_eles_tets; eles_pris mesh_eles_pris; int n_int_inter_types; int n_bdy_inter_types; array<int_inters> mesh_int_inters; array<bdy_inters> mesh_bdy_inters; int rank; /*! No-slip wall flux point coordinates for wall models. */ array< array<double> > loc_noslip_bdy; /*! Diagnostic output quantities. */ array<double> body_force; array<double> inv_force; array<double> vis_force; array<double> norm_residual; array<double> integral_quantities; double coeff_lift; double coeff_drag; /*! Plotting resolution. */ int p_res; #ifdef _MPI int nproc; int n_mpi_inter_types; array<mpi_inters> mesh_mpi_inters; array<int> error_states; int n_mpi_inters; /*! No-slip wall flux point coordinates for wall models. */ array< array<double> > loc_noslip_bdy_global; #endif };
1.4 eles
1.4.1 members
protected: // #### members #### /// flag to avoid re-setting-up transform arrays bool first_time; /*! mesh motion flag */ int motion; /*! viscous flag */ int viscous; /*! LES flag */ int LES; /*! SGS model */ int sgs_model; /*! LES filter flag */ int filter; /*! near-wall model */ int wall_model; /*! number of elements */ int n_eles; /*! number of elements that have a boundary face*/ int n_bdy_eles; /*! number of dimensions */ int n_dims; /*! number of prognostic fields */ int n_fields; /*! number of diagnostic fields */ int n_diagnostic_fields; /*! number of time averaged diagnostic fields */ int n_average_fields; /*! order of solution polynomials */ int order; /*! order of interface cubature rule */ int inters_cub_order; /*! order of interface cubature rule */ int volume_cub_order; /*! order of solution polynomials in restart file*/ int order_rest; /*! number of solution points per element */ int n_upts_per_ele; /*! number of solution points per element */ int n_upts_per_ele_rest; /*! number of flux points per element */ int n_fpts_per_ele; /*! number of vertices per element */ int n_verts_per_ele; array<int> connectivity_plot; /*! plotting resolution */ int p_res; /*! solution point type */ int upts_type; /*! flux point type */ int fpts_type; /*! number of plot points per element */ int n_ppts_per_ele; /*! number of plot elements per element */ int n_peles_per_ele; /*! Global cell number of element */ array<int> ele2global_ele; /*! Global cell number of element */ array<int> bdy_ele2ele; /*! Boundary condition type of faces */ array<int> bctype; /*! number of shape points per element */ array<int> n_spts_per_ele; /*! transformed normal at flux points */ array<double> tnorm_fpts; /*! transformed normal at flux points */ array< array<double> > tnorm_inters_cubpts; /*! location of solution points in standard element */ array<double> loc_upts; /*! location of solution points in standard element */ array<double> loc_upts_rest; /*! location of flux points in standard element */ array<double> tloc_fpts; /*! location of interface cubature points in standard element */ array< array<double> > loc_inters_cubpts; /*! weight of interface cubature points in standard element */ array< array<double> > weight_inters_cubpts; /*! location of volume cubature points in standard element */ array<double> loc_volume_cubpts; /*! weight of cubature points in standard element */ array<double> weight_volume_cubpts; /*! transformed normal at cubature points */ array< array<double> > tnorm_cubpts; /*! location of plot points in standard element */ array<double> loc_ppts; /*! location of shape points in standard element (simplex elements only)*/ array<double> loc_spts; /*! number of interfaces per element */ int n_inters_per_ele; /*! number of flux points per interface */ array<int> n_fpts_per_inter; /*! number of cubature points per interface */ array<int> n_cubpts_per_inter; /*! number of cubature points per interface */ int n_cubpts_per_ele; /*! element type (0=>quad,1=>tri,2=>tet,3=>pri,4=>hex) */ int ele_type; /*! order of polynomials defining shapes */ int s_order; /*! maximum number of shape points used by any element */ int max_n_spts_per_ele; /*! position of shape points (mesh vertices) in static-physical domain */ array<double> shape; /*! position of shape points (mesh vertices) in dynamic-physical domain */ array<double> shape_dyn; /*! Description: Mesh velocity at shape points \n indexing: (in_ele)(in_spt, in_dim) \n */ array<double> vel_spts; /*! Description: Mesh velocity at flux points (interpolated using shape basis funcs) \n indexing: (in_dim, in_fpt, in_ele) \n */ array<double> grid_vel_upts, grid_vel_fpts, vel_ppts; /*! nodal shape basis contributions at flux points */ array<double> nodal_s_basis_fpts; /*! nodal shape basis contributions at solution points */ array<double> nodal_s_basis_upts; /*! nodal shape basis contributions at output plot points */ array<double> nodal_s_basis_ppts; /*! nodal shape basis contributions at output plot points */ array<double> nodal_s_basis_vol_cubpts; /*! nodal shape basis contributions at output plot points */ array<array<double> > nodal_s_basis_inters_cubpts; /*! nodal shape basis derivative contributions at flux points */ array<double> d_nodal_s_basis_fpts; /*! nodal shape basis derivative contributions at solution points */ array<double> d_nodal_s_basis_upts; /*! nodal shape basis contributions at output plot points */ array<double> d_nodal_s_basis_vol_cubpts; /*! nodal shape basis contributions at output plot points */ array<array<double> > d_nodal_s_basis_inters_cubpts; /*! temporary solution storage at a single solution point */ array<double> temp_u; /*! temporary grid velocity storage at a single solution point */ array<double> temp_v; /*! temporary grid velocity storage at a single solution point (transformed to static frame) */ array<double> temp_v_ref; /*! temporary flux storage for GCL at a single solution point (transformed to static frame) */ array<double> temp_f_GCL; /*! temporary flux storage for GCL at a single solution point (transformed to static frame) */ array<double> temp_f_ref_GCL; /*! constansts for RK time-stepping */ array<double> RK_a, RK_b, RK_c; /*! temporary solution gradient storage */ array<double> temp_grad_u; /*! Matrix of filter weights at solution points */ array<double> filter_upts; /*! extra arrays for similarity model: Leonard tensors, velocity/energy products */ array<double> Lu, Le, uu, ue; /*! temporary flux storage */ array<double> temp_f; /*! temporary flux storage */ array<double> temp_f_ref; /*! temporary subgrid-scale flux storage */ array<double> temp_sgsf; /*! temporary subgrid-scale flux storage for dynamic->static transformation */ array<double> temp_sgsf_ref; /*! storage for distance of solution points to nearest no-slip boundary */ array<double> wall_distance; array<double> wall_distance_mag; array<double> twall; /*! number of storage levels for time-integration scheme */ int n_adv_levels; /*! determinant of Jacobian (transformation matrix) at solution points * (J = |G|) */ array<double> detjac_upts; /*! determinant of Jacobian (transformation matrix) at flux points * (J = |G|) */ array<double> detjac_fpts; /*! determinant of jacobian at volume cubature points. TODO: what is this really? */ array< array<double> > vol_detjac_inters_cubpts; /*! determinant of volume jacobian at cubature points. TODO: what is this really? */ array< array<double> > vol_detjac_vol_cubpts; /*! Full vector-transform matrix from static physical->computational frame, at solution points * [Determinant of Jacobian times inverse of Jacobian] [J*G^-1] */ array<double> JGinv_upts; /*! Full vector-transform matrix from static physical->computational frame, at flux points * [Determinant of Jacobian times inverse of Jacobian] [J*G^-1] */ array<double> JGinv_fpts; /*! Magnitude of transformed face-area normal vector from computational -> static-physical frame * [magntiude of (normal dot inverse static transformation matrix)] [ |J*(G^-1)*(n*dA)| ] */ array<double> tdA_fpts; /*! determinant of interface jacobian at flux points */ array< array<double> > inter_detjac_inters_cubpts; /*! normal at flux points*/ array<double> norm_fpts; /*! static-physical coordinates at flux points*/ array<double> pos_fpts; /*! static-physical coordinates at solution points*/ array<double> pos_upts; /*! normal at interface cubature points*/ array< array<double> > norm_inters_cubpts; /*! determinant of dynamic jacobian at solution points ( |G| ) */ array<double> J_dyn_upts; /*! determinant of dynamic jacobian at flux points ( |G| ) */ array<double> J_dyn_fpts; /*! Dynamic transformation matrix at solution points ( |G|*G^-1 ) */ array<double> JGinv_dyn_upts; /*! Dynamic->Static transformation matrix at flux points ( |G|*G^-1 ) */ array<double> JGinv_dyn_fpts; /*! Static->Dynamic transformation matrix at flux points ( G/|G| ) */ array<double> JinvG_dyn_fpts; /*! transformed gradient of determinant of dynamic jacobian at solution points */ array<double> tgrad_J_dyn_upts; /*! transformed gradient of determinant of dynamic jacobian at flux points */ array<double> tgrad_J_dyn_fpts; /*! normal at flux points in dynamic mesh */ array<double> norm_dyn_fpts; /*! physical coordinates at flux points in dynamic mesh */ array<double> dyn_pos_fpts, dyn_pos_upts; /*! magnitude of transformed face-area normal vector from static-physical -> dynamic-physical frame * [magntiude of (normal dot inverse dynamic transformation matrix)] [ |J*(G^-1)*(n*dA)| ] */ array<double> ndA_dyn_fpts; /*! description: transformed discontinuous solution at the solution points indexing: \n matrix mapping: */ array< array<double> > disu_upts; /*! running time-averaged diagnostic fields at solution points */ array<double> disu_average_upts; /*! time (in secs) until start of time average period for above diagnostic fields */ double spinup_time; /*! filtered solution at solution points for similarity and SVV LES models */ array<double> disuf_upts; /*! position at the plot points */ array< array<double> > pos_ppts; /*! description: transformed discontinuous solution at the flux points \n indexing: (in_fpt, in_field, in_ele) \n matrix mapping: (in_fpt || in_field, in_ele) */ array<double> disu_fpts; /*! description: transformed discontinuous flux at the solution points \n indexing: (in_upt, in_dim, in_field, in_ele) \n matrix mapping: (in_upt, in_dim || in_field, in_ele) */ array<double> tdisf_upts; /*! description: subgrid-scale flux at the solution points \n indexing: (in_upt, in_dim, in_field, in_ele) \n matrix mapping: (in_upt, in_dim || in_field, in_ele) */ array<double> sgsf_upts; /*! description: subgrid-scale flux at the flux points \n indexing: (in_fpt, in_dim, in_field, in_ele) \n matrix mapping: (in_fpt, in_dim || in_field, in_ele) */ array<double> sgsf_fpts; /*! normal transformed discontinuous flux at the flux points indexing: \n matrix mapping: */ array<double> norm_tdisf_fpts; /*! normal transformed continuous flux at the flux points indexing: \n matrix mapping: */ array<double> norm_tconf_fpts; /*! divergence of transformed continuous flux at the solution points indexing: \n matrix mapping: */ array< array<double> > div_tconf_upts; /*! delta of the transformed discontinuous solution at the flux points */ array<double> delta_disu_fpts; /*! gradient of discontinuous solution at solution points */ array<double> grad_disu_upts; /*! gradient of discontinuous solution at flux points */ array<double> grad_disu_fpts; /*! transformed discontinuous viscous flux at the solution points */ //array<double> tdisvisf_upts; /*! normal transformed discontinuous viscous flux at the flux points */ //array<double> norm_tdisvisf_fpts; /*! normal transformed continuous viscous flux at the flux points */ //array<double> norm_tconvisf_fpts; /*! transformed gradient of determinant of jacobian at solution points */ array<double> tgrad_detjac_upts; /*! source term for SA turbulence model at solution points */ array<double> src_upts; array<double> d_nodal_s_basis; // TODO: change naming (comments) to reflect reuse /*! * description: discontinuous solution for GCL at the solution points \n * indexing: (in_upt, in_ele) \n */ array< array<double> > Jbar_upts; /*! * description: discontinuous solution for GCL at the flux points \n * indexing: (in_fpt, in_ele) \n */ array< array<double> > Jbar_fpts; /*! * description: transformed discontinuous flux for GCL at the solution points \n * indexing: (in_upt, in_dim, in_ele) \n */ array<double> tdisf_GCL_upts; /*! * description: transformed discontinuous flux for GCL at the flux points \n * indexing: (in_fpt, in_dim, in_ele) \n */ array<double> tdisf_GCL_fpts; /*! * normal transformed discontinuous GCL flux at the flux points \n * indexing: (in_fpt, in_ele) \n */ array<double> norm_tdisf_GCL_fpts; /*! * normal transformed continuous GCL flux at the flux points \n * indexing: (in_fpt, in_ele) \n */ array<double> norm_tconf_GCL_fpts; /*! * divergence of transformed continuous GCL flux at the solution points * indexing: (in_upt, in_ele) \n */ array< array<double> > div_tconf_GCL_upts; #ifdef _GPU cusparseHandle_t handle; #endif /*! operator to go from transformed discontinuous solution at the solution points to transformed discontinuous solution at the flux points */ array<double> opp_0; array<double> opp_0_data; array<int> opp_0_cols; array<int> opp_0_b; array<int> opp_0_e; int opp_0_sparse; #ifdef _GPU array<double> opp_0_ell_data; array<int> opp_0_ell_indices; int opp_0_nnz_per_row; #endif /*! operator to go from transformed discontinuous inviscid flux at the solution points to divergence of transformed discontinuous inviscid flux at the solution points */ array< array<double> > opp_1; array< array<double> > opp_1_data; array< array<int> > opp_1_cols; array< array<int> > opp_1_b; array< array<int> > opp_1_e; int opp_1_sparse; #ifdef _GPU array< array<double> > opp_1_ell_data; array< array<int> > opp_1_ell_indices; array<int> opp_1_nnz_per_row; #endif /*! operator to go from transformed discontinuous inviscid flux at the solution points to normal transformed discontinuous inviscid flux at the flux points */ array< array<double> > opp_2; array< array<double> > opp_2_data; array< array<int> > opp_2_cols; array< array<int> > opp_2_b; array< array<int> > opp_2_e; int opp_2_sparse; #ifdef _GPU array< array<double> > opp_2_ell_data; array< array<int> > opp_2_ell_indices; array<int> opp_2_nnz_per_row; #endif /*! operator to go from normal correction inviscid flux at the flux points to divergence of correction inviscid flux at the solution points*/ array<double> opp_3; array<double> opp_3_data; array<int> opp_3_cols; array<int> opp_3_b; array<int> opp_3_e; int opp_3_sparse; #ifdef _GPU array<double> opp_3_ell_data; array<int> opp_3_ell_indices; int opp_3_nnz_per_row; #endif /*! operator to go from transformed solution at solution points to transformed gradient of transformed solution at solution points */ array< array<double> > opp_4; array< array<double> > opp_4_data; array< array<int> > opp_4_cols; array< array<int> > opp_4_b; array< array<int> > opp_4_e; int opp_4_sparse; #ifdef _GPU array< array<double> > opp_4_ell_data; array< array<int> > opp_4_ell_indices; array< int > opp_4_nnz_per_row; #endif /*! operator to go from transformed solution at flux points to transformed gradient of transformed solution at solution points */ array< array<double> > opp_5; array< array<double> > opp_5_data; array< array<int> > opp_5_cols; array< array<int> > opp_5_b; array< array<int> > opp_5_e; int opp_5_sparse; #ifdef _GPU array< array<double> > opp_5_ell_data; array< array<int> > opp_5_ell_indices; array<int> opp_5_nnz_per_row; #endif /*! operator to go from transformed solution at solution points to transformed gradient of transformed solution at flux points */ array<double> opp_6; array<double> opp_6_data; array<int> opp_6_cols; array<int> opp_6_b; array<int> opp_6_e; int opp_6_sparse; #ifdef _GPU array<double> opp_6_ell_data; array<int> opp_6_ell_indices; int opp_6_nnz_per_row; #endif /*! operator to go from discontinuous solution at the solution points to discontinuous solution at the plot points */ array<double> opp_p; array< array<double> > opp_inters_cubpts; array<double> opp_volume_cubpts; /*! operator to go from discontinuous solution at the restart points to discontinuous solution at the solutoin points */ array<double> opp_r; /*! dimensions for blas calls */ int Arows, Acols; int Brows, Bcols; int Astride, Bstride, Cstride; /*! general settings for mkl sparse blas */ char matdescra[6]; /*! transpose setting for mkl sparse blas */ char transa; /*! zero for mkl sparse blas */ double zero; /*! one for mkl sparse blas */ double one; /*! number of fields multiplied by number of elements */ int n_fields_mul_n_eles; /*! number of dimensions multiplied by number of solution points per element */ int n_dims_mul_n_upts_per_ele; int rank; int nproc; /*! mass flux through inlet */ double mass_flux; /*! reference element length */ array<double> h_ref; /*! element local timestep */ array<double> dt_local; double dt_local_new; array<double> dt_local_mpi; /*! Artificial Viscosity variables */ array<double> vandermonde; array<double> inv_vandermonde; array<double> vandermonde2D; array<double> inv_vandermonde2D; array<double> area_coord_upts; array<double> area_coord_fpts; array<double> epsilon; array<double> epsilon_upts; array<double> epsilon_fpts; array<double> concentration_array; array<double> sensor; array<double> sigma; array<double> min_dt_local; /*! Global cell number of element as in the code */ array<int> ele2global_ele_code;
1.4.2 methods
// #### methods #### /*! setup */ void setup(int in_n_eles, int in_max_s_spts_per_ele); /*! setup initial conditions */ void set_ics(double& time); /*! read data from restart file */ void read_restart_data(ifstream& restart_file); /*! write data to restart file */ void write_restart_data(ofstream& restart_file); /*! write extra restart file containing x,y,z of solution points instead of solution data */ void write_restart_mesh(ofstream& restart_file); /*! move all to from cpu to gpu */ void mv_all_cpu_gpu(void); /*! move wall distance array to from cpu to gpu */ void mv_wall_distance_cpu_gpu(void); /*! move wall distance magnitude array to from cpu to gpu */ void mv_wall_distance_mag_cpu_gpu(void); /*! copy transformed discontinuous solution at solution points to cpu */ void cp_disu_upts_gpu_cpu(void); /*! copy transformed discontinuous solution at solution points to gpu */ void cp_disu_upts_cpu_gpu(void); void cp_grad_disu_upts_gpu_cpu(void); /*! copy determinant of jacobian at solution points to cpu */ void cp_detjac_upts_gpu_cpu(void); /*! copy divergence at solution points to cpu */ void cp_div_tconf_upts_gpu_cpu(void); /*! copy local time stepping reference length at solution points to cpu */ void cp_h_ref_gpu_cpu(void); /*! copy source term at solution points to cpu */ void cp_src_upts_gpu_cpu(void); /*! copy elemental sensor values to cpu */ void cp_sensor_gpu_cpu(void); /*! copy AV co-eff values at solution points to cpu */ void cp_epsilon_upts_gpu_cpu(void); /*! remove transformed discontinuous solution at solution points from cpu */ void rm_disu_upts_cpu(void); /*! remove determinant of jacobian at solution points from cpu */ void rm_detjac_upts_cpu(void); /*! calculate the discontinuous solution at the flux points */ void extrapolate_solution(int in_disu_upts_from); /*! Calculate terms for some LES models */ void calc_sgs_terms(int in_disu_upts_from); /*! calculate transformed discontinuous inviscid flux at solution points */ void evaluate_invFlux(int in_disu_upts_from); /*! calculate divergence of transformed discontinuous flux at solution points */ void calculate_divergence(int in_div_tconf_upts_to); /*! calculate normal transformed discontinuous flux at flux points */ void extrapolate_totalFlux(void); /*! calculate subgrid-scale flux at flux points */ void evaluate_sgsFlux(void); /*! calculate divergence of transformed continuous flux at solution points */ void calculate_corrected_divergence(int in_div_tconf_upts_to); /*! calculate uncorrected transformed gradient of the discontinuous solution at the solution points */ void calculate_gradient(int in_disu_upts_from); /*! calculate corrected gradient of the discontinuous solution at solution points */ void correct_gradient(void); /*! calculate corrected gradient of the discontinuous solution at flux points */ void extrapolate_corrected_gradient(void); /*! calculate corrected gradient of solution at flux points */ //void extrapolate_corrected_gradient(void); /*! calculate transformed discontinuous viscous flux at solution points */ void evaluate_viscFlux(int in_disu_upts_from); /*! calculate divergence of transformed discontinuous viscous flux at solution points */ //void calc_div_tdisvisf_upts(int in_div_tconinvf_upts_to); /*! calculate normal transformed discontinuous viscous flux at flux points */ //void calc_norm_tdisvisf_fpts(void); /*! calculate divergence of transformed continuous viscous flux at solution points */ //void calc_div_tconvisf_upts(int in_div_tconinvf_upts_to); /*! calculate source term for SA turbulence model at solution points */ void calc_src_upts_SA(int in_disu_upts_from); /*! advance solution using a runge-kutta scheme */ void AdvanceSolution(int in_step, int adv_type); /*! Calculate element local timestep */ double calc_dt_local(int in_ele); /*! get number of elements */ int get_n_eles(void); // get number of ppts_per_ele int get_n_ppts_per_ele(void); // get number of peles_per_ele int get_n_peles_per_ele(void); // get number of verts_per_ele int get_n_verts_per_ele(void); /*! get number of solution points per element */ int get_n_upts_per_ele(void); /*! get element type */ int get_ele_type(void); /*! get number of dimensions */ int get_n_dims(void); /*! get number of fields */ int get_n_fields(void); /*! set shape */ void set_shape(int in_max_n_spts_per_ele); /*! set shape node */ void set_shape_node(int in_spt, int in_ele, array<double>& in_pos); /*! Set new position of shape node in dynamic domain */ void set_dynamic_shape_node(int in_spt, int in_ele, array<double> &in_pos); /*! set bc type */ void set_bctype(int in_ele, int in_inter, int in_bctype); /*! set bc type */ void set_bdy_ele2ele(void); /*! set number of shape points */ void set_n_spts(int in_ele, int in_n_spts); /*! set global element number */ void set_ele2global_ele(int in_ele, int in_global_ele); /*! get a pointer to the transformed discontinuous solution at a flux point */ double* get_disu_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_field, int in_ele); /*! get a pointer to the normal transformed continuous flux at a flux point */ double* get_norm_tconf_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_field, int in_ele); /*! get a pointer to the determinant of the jacobian at a flux point (static->computational) */ double* get_detjac_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_ele); /*! get a pointer to the determinant of the jacobian at a flux point (dynamic->static) */ double* get_detjac_dyn_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_ele); /*! get pointer to the equivalent of 'dA' (face area) at a flux point in static physical space */ double* get_tdA_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_ele); /*! get pointer to the equivalent of 'dA' (face area) at a flux point in dynamic physical space */ double* get_ndA_dyn_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_ele); /*! get a pointer to the normal at a flux point */ double* get_norm_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_dim, int in_ele); /*! get a pointer to the normal at a flux point in dynamic space */ double* get_norm_dyn_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_dim, int in_ele); /*! get a CPU pointer to the coordinates at a flux point */ double* get_loc_fpts_ptr_cpu(int in_inter_local_fpt, int in_ele_local_inter, int in_dim, int in_ele); /*! get a GPU pointer to the coordinates at a flux point */ double* get_loc_fpts_ptr_gpu(int in_inter_local_fpt, int in_ele_local_inter, int in_dim, int in_ele); /*! get a CPU pointer to the dynamic physical coordinates at a flux point */ double* get_pos_dyn_fpts_ptr_cpu(int in_inter_local_fpt, int in_ele_local_inter, int in_dim, int in_ele); /*! get a pointer to delta of the transformed discontinuous solution at a flux point */ double* get_delta_disu_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_field, int in_ele); /*! get a pointer to gradient of discontinuous solution at a flux point */ double* get_grad_disu_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_dim, int in_field, int in_ele); /*! get a pointer to gradient of discontinuous solution at a flux point */ double* get_normal_disu_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_field, int in_ele, array<double> temp_loc, double temp_pos[3]); /*! get a pointer to the normal transformed continuous viscous flux at a flux point */ //double* get_norm_tconvisf_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_field, int in_ele); /*! get a pointer to the subgrid-scale flux at a flux point */ double* get_sgsf_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_field, int in_dim, int in_ele); /*! set opp_0 */ void set_opp_0(int in_sparse); /*! set opp_1 */ void set_opp_1(int in_sparse); /*! set opp_2 */ void set_opp_2(int in_sparse); /*! set opp_3 */ void set_opp_3(int in_sparse); /*! set opp_4 */ void set_opp_4(int in_sparse); /*! set opp_5 */ void set_opp_5(int in_sparse); /*! set opp_6 */ void set_opp_6(int in_sparse); /*! set opp_p */ void set_opp_p(void); /*! set opp_p */ void set_opp_inters_cubpts(void); /*! set opp_p */ void set_opp_volume_cubpts(void); /*! set opp_r */ void set_opp_r(void); /*! calculate position of the plot points */ void calc_pos_ppts(int in_ele, array<double>& out_pos_ppts); void set_rank(int in_rank); virtual void set_connectivity_plot()=0; void set_disu_upts_to_zero_other_levels(void); array<int> get_connectivity_plot(); /*! calculate solution at the plot points */ void calc_disu_ppts(int in_ele, array<double>& out_disu_ppts); /*! calculate gradient of solution at the plot points */ void calc_grad_disu_ppts(int in_ele, array<double>& out_grad_disu_ppts); /*! calculate sensor at the plot points */ void calc_sensor_ppts(int in_ele, array<double>& out_sensor_ppts); /*! calculate AV-co-efficients at the plot points */ void calc_epsilon_ppts(int in_ele, array<double>& out_epsilon_ppts); /*! calculate time-averaged diagnostic fields at the plot points */ void calc_time_average_ppts(int in_ele, array<double>& out_disu_average_ppts); /*! calculate diagnostic fields at the plot points */ void calc_diagnostic_fields_ppts(int in_ele, array<double>& in_disu_ppts, array<double>& in_grad_disu_ppts, array<double>& in_sensor_ppts, array<double> &in_epsilon_ppts, array<double>& out_diag_field_ppts, double& time); /*! calculate position of a solution point */ void calc_pos_upt(int in_upt, int in_ele, array<double>& out_pos); /*! get physical position of a flux point */ void calc_pos_fpt(int in_fpt, int in_ele, array<double>& out_pos); /*! returns position of a solution point */ double get_loc_upt(int in_upt, int in_dim); /*! set transforms */ void set_transforms(void); /*! set transforms at the interface cubature points */ void set_transforms_inters_cubpts(void); /*! set transforms at the volume cubature points */ void set_transforms_vol_cubpts(void); /*! Calculate distance of solution points to no-slip wall */ void calc_wall_distance(int n_seg_noslip_inters, int n_tri_noslip_inters, int n_quad_noslip_inters, array< array<double> > loc_noslip_bdy); /*! Calculate distance of solution points to no-slip wall in parallel */ void calc_wall_distance_parallel(array<int> n_seg_noslip_inters, array<int> n_tri_noslip_inters, array<int> n_quad_noslip_inters, array< array<double> > loc_noslip_bdy_global, int nproc); /*! calculate position */ void calc_pos(array<double> in_loc, int in_ele, array<double>& out_pos); /*! calculate derivative of position */ void calc_d_pos(array<double> in_loc, int in_ele, array<double>& out_d_pos); /*! calculate derivative of position at a solution point (using pre-computed gradients) */ void calc_d_pos_upt(int in_upt, int in_ele, array<double>& out_d_pos); /*! calculate derivative of position at a flux point (using pre-computed gradients) */ void calc_d_pos_fpt(int in_fpt, int in_ele, array<double>& out_d_pos); // #### virtual methods #### virtual void setup_ele_type_specific()=0; /*! prototype for element reference length calculation */ virtual double calc_h_ref_specific(int in_eles) = 0; virtual int read_restart_info(ifstream& restart_file)=0; virtual void write_restart_info(ofstream& restart_file)=0; /*! Compute interface jacobian determinant on face */ virtual double compute_inter_detjac_inters_cubpts(int in_inter, array<double> d_pos)=0; /*! evaluate nodal basis */ virtual double eval_nodal_basis(int in_index, array<double> in_loc)=0; /*! evaluate nodal basis for restart file*/ virtual double eval_nodal_basis_restart(int in_index, array<double> in_loc)=0; /*! evaluate derivative of nodal basis */ virtual double eval_d_nodal_basis(int in_index, int in_cpnt, array<double> in_loc)=0; virtual void fill_opp_3(array<double>& opp_3)=0; /*! evaluate divergence of vcjh basis */ //virtual double eval_div_vcjh_basis(int in_index, array<double>& loc)=0; /*! evaluate nodal shape basis */ virtual double eval_nodal_s_basis(int in_index, array<double> in_loc, int in_n_spts)=0; /*! evaluate derivative of nodal shape basis */ virtual void eval_d_nodal_s_basis(array<double> &d_nodal_s_basis, array<double> in_loc, int in_n_spts)=0; /*! Calculate SGS flux */ void calc_sgsf_upts(array<double>& temp_u, array<double>& temp_grad_u, double& detjac, int ele, int upt, array<double>& temp_sgsf); /*! rotate velocity components to surface*/ array<double> calc_rotation_matrix(array<double>& norm); /*! calculate wall shear stress using LES wall model*/ void calc_wall_stress(double rho, array<double>& urot, double ene, double mu, double Pr, double gamma, double y, array<double>& tau_wall, double q_wall); /*! Wall function calculator for Breuer-Rodi wall model */ double wallfn_br(double yplus, double A, double B, double E, double kappa); /*! Calculate element volume */ virtual double calc_ele_vol(double& detjac)=0; double compute_res_upts(int in_norm_type, int in_field); /*! calculate body forcing at solution points */ void evaluate_body_force(int in_file_num); /*! Compute volume integral of diagnostic quantities */ void CalcIntegralQuantities(int n_integral_quantities, array <double>& integral_quantities); /*! Compute time-average diagnostic quantities */ void CalcTimeAverageQuantities(double& time); void compute_wall_forces(array<double>& inv_force, array<double>& vis_force, double& temp_cl, double& temp_cd, ofstream& coeff_file, bool write_forces); array<double> compute_error(int in_norm_type, double& time); array<double> get_pointwise_error(array<double>& sol, array<double>& grad_sol, array<double>& loc, double& time, int in_norm_type); /*! calculate position of a point in physical (dynamic) space from (r,s,t) coordinates*/ void calc_pos_dyn(array<double> in_loc, int in_ele, array<double> &out_pos); /*! * Calculate dynamic position of solution point * \param[in] in_upt - ID of solution point within element to evaluate at * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_pos_dyn_fpt(int in_fpt, int in_ele, array<double> &out_pos); /*! * Calculate dynamic position of flux point * \param[in] in_upt - ID of solution point within element to evaluate at * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_pos_dyn_upt(int in_upt, int in_ele, array<double> &out_pos); /*! * Calculate dynamic position of plot point * \param[in] in_ppt - ID of plot point within element to evaluate at * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_pos_dyn_ppt(int in_ppt, int in_ele, array<double> &out_pos); /*! * Calculate dynamic position of volume cubature point * \param[in] in_cubpt - ID of cubature point within element to evaluate at * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_pos_dyn_vol_cubpt(int in_cubpt, int in_ele, array<double> &out_pos); /*! * Calculate dynamic position of interface cubature point * \param[in] in_cubpt - ID of cubature point on element face to evaluate at * \param[in] in_face - local face ID within element * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_pos_dyn_inters_cubpt(int in_cubpt, int in_face, int in_ele, array<double> &out_pos); /*! * Calculate derivative of dynamic position wrt reference (initial,static) position * \param[in] in_loc - position of point in computational space * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_d_pos_dyn(array<double> in_loc, int in_ele, array<double> &out_d_pos); /*! * Calculate derivative of dynamic position wrt reference (initial,static) position at fpt * \param[in] in_fpt - ID of flux point within element to evaluate at * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_d_pos_dyn_fpt(int in_fpt, int in_ele, array<double> &out_d_pos); /*! * Calculate derivative of dynamic position wrt reference (initial,static) position at upt * \param[in] in_upt - ID of solution point within element to evaluate at * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_d_pos_dyn_upt(int in_upt, int in_ele, array<double> &out_d_pos); /*! * Calculate derivative of dynamic position wrt reference (initial,static) position at volume cubature point * \param[in] in_cubpt - ID of cubature point within element to evaluate at * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_d_pos_dyn_vol_cubpt(int in_cubpt, int in_ele, array<double> &out_d_pos); /*! * Calculate derivative of dynamic position wrt reference (initial,static) position at volume cubature point * \param[in] in_cubpt - ID of cubature point on face to evaluate at * \param[in] in_face - local face ID within element * \param[in] in_ele - local element ID * \param[out] out_d_pos - array of size (n_dims,n_dims); (i,j) = dx_i / dX_j */ void calc_d_pos_dyn_inters_cubpt(int in_cubpt, int in_face, int in_ele, array<double> &out_d_pos); /*! pre-computing shape basis contributions at flux points for more efficient access */ void store_nodal_s_basis_fpts(void); /*! pre-computing shape basis contributions at solution points for more efficient access */ void store_nodal_s_basis_upts(void); /*! pre-computing shape basis contributions at plot points for more efficient access */ void store_nodal_s_basis_ppts(void); /*! pre-computing shape basis contributions at plot points for more efficient access */ void store_nodal_s_basis_vol_cubpts(void); /*! pre-computing shape basis contributions at plot points for more efficient access */ void store_nodal_s_basis_inters_cubpts(void); /*! pre-computing shape basis deriavative contributions at flux points for more efficient access */ void store_d_nodal_s_basis_fpts(void); /*! pre-computing shape basis derivative contributions at solution points for more efficient access */ void store_d_nodal_s_basis_upts(void); /*! pre-computing shape basis derivative contributions at solution points for more efficient access */ void store_d_nodal_s_basis_vol_cubpts(void); /*! pre-computing shape basis derivative contributions at solution points for more efficient access */ void store_d_nodal_s_basis_inters_cubpts(void); /*! initialize arrays for storing grid velocities */ void initialize_grid_vel(int in_max_n_spts_per_ele); /*! set grid velocity on element shape points */ void set_grid_vel_spt(int in_ele, int in_spt, array<double> in_vel); /*! interpolate grid velocity from shape points to flux points */ void set_grid_vel_fpts(int in_rk_step); /*! interpolate grid velocity from shape points to solution points */ void set_grid_vel_upts(int in_rk_step); /*! interpolate grid velocity from shape points to plot points */ void set_grid_vel_ppts(void); /*! Get array of grid velocity at all plot points */ array<double> get_grid_vel_ppts(void); /*! Get pointer to grid velocity at a flux point */ double *get_grid_vel_fpts_ptr(int in_ele, int in_ele_local_inter, int in_inter_local_fpt, int in_dim); /*! Set the transformation variables for dynamic-physical -> static-physical frames */ void set_transforms_dynamic(void); /* --- Geometric Conservation Law (GCL) Funcitons --- */ /*! Update the dynamic transformation variables with the GCL-corrected Jacobian determinant */ void correct_dynamic_transforms(void); /*! GCL Residual-Calculation Steps */ void evaluate_GCL_flux(int in_disu_upts_from); void extrapolate_GCL_solution(int in_disu_upts_from); void extrapolate_GCL_flux(void); void calculate_divergence_GCL(int in_div_tconf_upts_to); void calculate_corrected_divergence_GCL(int in_div_tconf_upts_to); double *get_disu_GCL_fpts_ptr(int in_inter_local_fpt, int in_ele_local_inter, int in_ele); /* --------------------------------------------------- */ /*! Set the time step for the current iteration */ void set_dt(int in_step, int adv_type);
1.5 elestris
1.5.1 members
protected: // methods void set_vandermonde(); void set_vandermonde_restart(); // members array<double> vandermonde; array<double> vandermonde_rest; array<double> inv_vandermonde; array<double> inv_vandermonde_rest; array<double> loc_1d_fpts;
1.5.2 methods
public: // #### constructors #### // default constructor eles_tris(); // #### methods #### /*! set shape */ //void set_shape(array<int> &in_n_spts_per_ele); void set_connectivity_plot(); /*! set location of solution points */ void set_loc_upts(void); /*! set location of flux points */ void set_tloc_fpts(void); /*! set location and weights of interface cubature points */ void set_inters_cubpts(void); /*! set location and weights of volume cubature points */ void set_volume_cubpts(void); /*! set location of plot points */ void set_loc_ppts(void); /*! set location of shape points */ // void set_loc_spts(void); /*! set transformed normals at flux points */ void set_tnorm_fpts(void); //#### helper methods #### void setup_ele_type_specific(void); /*! read restart info */ int read_restart_info(ifstream& restart_file); /*! write restart info */ void write_restart_info(ofstream& restart_file); /*! Compute interface jacobian determinant on face */ double compute_inter_detjac_inters_cubpts(int in_inter, array<double> d_pos); /*! evaluate nodal basis */ double eval_nodal_basis(int in_index, array<double> in_loc); /*! evaluate nodal basis for restart file*/ double eval_nodal_basis_restart(int in_index, array<double> in_loc); /*! evaluate derivative of nodal basis */ double eval_d_nodal_basis(int in_index, int in_cpnt, array<double> in_loc); /*! evaluate divergence of vcjh basis */ //double eval_div_vcjh_basis(int in_index, array<double>& loc); void fill_opp_3(array<double>& opp_3); /*! evaluate nodal shape basis */ double eval_nodal_s_basis(int in_index, array<double> in_loc, int in_n_spts); /*! evaluate derivative of nodal shape basis */ void eval_d_nodal_s_basis(array<double> &d_nodal_s_basis, array<double> in_loc, int in_n_spts); /*! Compute the filter matrix for subgrid-scale models */ void compute_filter_upts(void); /*! Calculate element volume */ double calc_ele_vol(double& detjac); /*! Element reference length calculation */ double calc_h_ref_specific(int in_ele);
[ ]void setlocupts(void)
1.6 inters
1.6.1 members
protected: // #### members #### int inters_type; // segment, quad or tri 界面的类型可能是线段（二维）、三角形或四边形（三维） int order; int viscous; int LES; int wall_model; int n_inters; //interfaces的个数，我的一个目的就是弄清楚如何得到这个。 int n_fpts_per_inter; //每个Interface上有多少个Flux points。 int n_fields; //流场的个数，二维的话是4个，三维是五个。其实就是U的个数。 int n_dims; //向量的维数，如果是二维的话，通量有两个分量，三维的话有三个分量。 int motion; //!< Mesh motion flag 网格是否运动的flag array<double*> disu_fpts_l; array<double*> delta_disu_fpts_l; array<double*> norm_tconf_fpts_l; //array<double*> norm_tconvisf_fpts_l; array<double*> detjac_fpts_l; array<double*> tdA_fpts_l; array<double*> norm_fpts; array<double*> pos_fpts; array<double*> pos_dyn_fpts; array<double> pos_disu_fpts_l; array<double*> grad_disu_fpts_l; array<double*> normal_disu_fpts_l; array<double> temp_u_l; array<double> temp_u_r; array<double> temp_grad_u_l; array<double> temp_grad_u_r; array<double> temp_normal_u_l; array<double> temp_pos_u_l; array<double> temp_f_l; array<double> temp_f_r; array<double> temp_fn_l; array<double> temp_fn_r; array<double> temp_f; array<double> temp_loc; // LES and wall model quantities array<double*> sgsf_fpts_l; array<double*> sgsf_fpts_r; array<double> temp_sgsf_l; array<double> temp_sgsf_r; array<int> lut; array<double> v_l, v_r, um, du; // Dynamic grid variables: // Note: grid velocity is continuous across interfaces array<double*> ndA_dyn_fpts_l; array<double*> norm_dyn_fpts; array<double*> J_dyn_fpts_l; array<double*> grid_vel_fpts; array<double*> disu_GCL_fpts_l; array<double*> norm_tconf_GCL_fpts_l; double temp_u_GCL_l; double temp_f_GCL_l; array<double> temp_v; array<double> temp_fn_ref_l; array<double> temp_fn_ref_r;
1.6.2 methods
public: // #### constructors #### // default constructor inters(); // default destructor ~inters(); // #### methods #### /*! setup inters */ void setup_inters(int in_n_inters, int in_inter_type); /*! Set normal flux to be normal * f_r */ void right_flux(array<double> &f_r, array<double> &norm, array<double> &fn, int n_dims, int n_fields, double gamma); /*! Compute common inviscid flux using Rusanov flux */ void rusanov_flux(array<double> &u_l, array<double> &u_r, array<double> &v_g, array<double> &f_l, array<double> &f_r, array<double> &norm, array<double> &fn, int n_dims, int n_fields, double gamma); /*! Compute common inviscid flux using Roe flux */ void roe_flux(array<double> &u_l, array<double> &u_r, array<double> &v_g, array<double> &norm, array<double> &fn, int n_dims, int n_fields, double gamma); /*! Compute common inviscid flux using Lax-Friedrich flux (works only for wave equation) */ void lax_friedrich(array<double> &u_l, array<double> &u_r, array<double> &norm, array<double> &fn, int n_dims, int n_fields, double lambda, array<double>& wave_speed); /*! Compute common viscous flux using LDG formulation */ void ldg_flux(int flux_spec, array<double> &u_l, array<double> &u_r, array<double> &f_l, array<double> &f_r, array<double> &norm, array<double> &fn, int n_dims, int n_fields, double tau, double pen_fact); /*! Compute common solution using LDG formulation */ void ldg_solution(int flux_spec, array<double> &u_l, array<double> &u_r, array<double> &u_c, double pen_fact, array<double>& norm); /*! get look up table for flux point connectivity based on rotation tag */ void get_lut(int in_rot_tag); /*! Compute common flux at boundaries using convective flux formulation */ void convective_flux_boundary(array<double> &f_l, array<double> &f_r, array<double> &norm, array<double> &fn, int n_dims, int n_fields);
1.7 intinters
1.7.1 members
1.7.2 methods
1.8 Out
1.8.1 members
没有members
1.8.2 methods
/*! write an output file in Tecplot ASCII format */ void write_tec(int in_file_num, struct solution* FlowSol); /*! write an output file in VTK ASCII format */ void write_vtu(int in_file_num, struct solution* FlowSol); /*! writing a restart file */ void write_restart(int in_file_num, struct solution* FlowSol); /*! compute forces on wall faces*/ void CalcForces(int in_file_num, struct solution* FlowSol); /*! compute integral diagnostic quantities */ void CalcIntegralQuantities(int in_file_num, struct solution* FlowSol); /*! Calculate time averaged diagnostic quantities */ void CalcTimeAverageQuantities(struct solution* FlowSol); /*! compute error */ void compute_error(int in_file_num, struct solution* FlowSol); /*! calculate residual */ void CalcNormResidual(struct solution* FlowSol); /*! monitor convergence of residual */ void HistoryOutput(int in_file_num, clock_t init, ofstream *write_hist, struct solution* FlowSol); /*! check if the solution is bounded !*/ void check_stability(struct solution* FlowSol);
2 函数
2.1 CalcResidual(FlowSol.iniiter+isteps, i, &FlowSol)
2.1.1 DONE eles::extrapolatesolution(int indisuuptsfrom)
- 代码
void eles::extrapolate_solution(int in_disu_upts_from) { if (n_eles!=0) { /*! Performs C = (alpha*A*B) + (beta*C) where: \n alpha = 1.0 \n beta = 0.0 \n A = opp_0 \n B = disu_upts(in_disu_upts_from) \n C = disu_fpts opp_0 is the polynomial extrapolation matrix; has dimensions n_f_pts_per_ele by n_upts_per_ele Recall: opp_0(j,i) = value of the ith nodal basis at the jth flux point location in the reference domain (vector of solution values at flux points) = opp_0 * (vector of solution values at nodes) */ Arows = n_fpts_per_ele; //插值矩阵行数等于每个单元上的Flux Point的个数 Acols = n_upts_per_ele;//插值矩阵的列数等于每个单元上Solution Point的个数 Brows = Acols; //B矩阵的行数等于插值矩阵的列数，因为这是每个Solution Point的值 Bcols = n_fields*n_eles;//B矩阵的列数单元的个数 Astride = Arows; //这三个我不大理解 Bstride = Brows; Cstride = Arows; #ifdef _CPU if(opp_0_sparse==0) // dense { #if defined _ACCELERATE_BLAS || defined _MKL_BLAS || defined _STANDARD_BLAS cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,Arows,Bcols,Acols,1.0,opp_0.get_ptr_cpu(),Astride,disu_upts(in_disu_upts_from).get_ptr_cpu(),Bstride,0.0,disu_fpts.get_ptr_cpu(),Cstride); #elif defined _NO_BLAS dgemm(Arows,Bcols,Acols,1.0,0.0,opp_0.get_ptr_cpu(),disu_upts(in_disu_upts_from).get_ptr_cpu(),disu_fpts.get_ptr_cpu()); #endif } else if(opp_0_sparse==1) // mkl blas four-array csr format { #if defined _MKL_BLAS mkl_dcsrmm(&transa,&n_fpts_per_ele,&n_fields_mul_n_eles,&n_upts_per_ele,&one,matdescra,opp_0_data.get_ptr_cpu(),opp_0_cols.get_ptr_cpu(),opp_0_b.get_ptr_cpu(),opp_0_e.get_ptr_cpu(),disu_upts(in_disu_upts_from).get_ptr_cpu(),&n_upts_per_ele,&zero,disu_fpts.get_ptr_cpu(),&n_fpts_per_ele); #endif } else { cout << "ERROR: Unknown storage for opp_0 ... " << endl; } #endif #ifdef _GPU if(opp_0_sparse==0) { cublasDgemm('N','N',Arows,Bcols,Acols,1.0,opp_0.get_ptr_gpu(),Astride,disu_upts(in_disu_upts_from).get_ptr_gpu(),Bstride,0.0,disu_fpts.get_ptr_gpu(),Cstride); } else if (opp_0_sparse==1) { bespoke_SPMV(n_fpts_per_ele,n_upts_per_ele,n_fields,n_eles,opp_0_ell_data.get_ptr_gpu(),opp_0_ell_indices.get_ptr_gpu(),opp_0_nnz_per_row,disu_upts(in_disu_upts_from).get_ptr_gpu(),disu_fpts.get_ptr_gpu(),ele_type,order,0); } else { cout << "ERROR: Unknown storage for opp_0 ... " << endl; } #endif } }
2.1.2 DONE void eles::evaluateinvFlux(int indisuuptsfrom)
- 代码
void eles::evaluate_invFlux(int in_disu_upts_from) { if (n_eles!=0) { #ifdef _CPU int i,j,k,l,m; for(i=0;i<n_eles;i++) { for(j=0;j<n_upts_per_ele;j++) { for(k=0;k<n_fields;k++) { temp_u(k)=disu_upts(in_disu_upts_from)(j,i,k); //temp_u(k) 是存储第i个单元上第j个点上的solution的临时变量 } if (motion) { // Transform solution from static frame to dynamic frame for (k=0; k<n_fields; k++) { temp_u(k) /= J_dyn_upts(j,i); } // Get mesh velocity in dynamic frame for (k=0; k<n_dims; k++) { temp_v(k) = grid_vel_upts(j,i,k); } // Temporary flux vector for dynamic->static transformation temp_f_ref.setup(n_fields,n_dims); }else{ temp_v.initialize_to_zero(); } if(n_dims==2) { calc_invf_2d(temp_u,temp_f);//第j个点上的无粘通量，并保存到二维数组temp_f(k,0)中 if (motion) calc_alef_2d(temp_u, temp_v, temp_f); } else if(n_dims==3) { calc_invf_3d(temp_u,temp_f); if (motion) calc_alef_3d(temp_u, temp_v, temp_f); } else { FatalError("Invalid number of dimensions!"); } // Transform from dynamic-physical space to static-physical space if (motion) { for(k=0; k<n_fields; k++) { for(l=0; l<n_dims; l++) { temp_f_ref(k,l)=0.; for(m=0; m<n_dims; m++) { temp_f_ref(k,l) += JGinv_dyn_upts(l,m,j,i)*temp_f(k,m); } } } // Copy Static-Physical Domain flux back to temp_f for (k=0; k<n_fields; k++) { for (l=0; l<n_dims; l++) { temp_f(k,l) = temp_f_ref(k,l); } } } // Transform from static physical space to computational space，将物理域转化为计算域，也就是乘以|J|/J^{-1},最终得到tdisf_upts(j,i,k,l),注意此处的m是为了区分通量的分量 for(k=0;k<n_fields;k++) { for(l=0;l<n_dims;l++) { tdisf_upts(j,i,k,l)=0.; for(m=0;m<n_dims;m++) { tdisf_upts(j,i,k,l) += JGinv_upts(l,m,j,i)*temp_f(k,m);//JGinv_upts(j,i,l,m)*temp_f(k,m); } } } } } #endif } }
2.1.3 DONE void intinters::calculatecommoninvFlux(void)
- 代码
void int_inters::calculate_common_invFlux(void) { #ifdef _CPU array<double> norm(n_dims), fn(n_fields); //viscous array<double> u_c(n_fields); for(int i=0;i<n_inters;i++) { for(int j=0;j<n_fpts_per_inter;j++) { // calculate discontinuous solution at flux points for(int k=0;k<n_fields;k++) { temp_u_l(k)=(*disu_fpts_l(j,i,k)); temp_u_r(k)=(*disu_fpts_r(j,i,k)); } if (motion) { // Transform solution to dynamic space for (int k=0; k<n_fields; k++) { temp_u_l(k) /= (*J_dyn_fpts_l(j,i)); temp_u_r(k) /= (*J_dyn_fpts_r(j,i)); } // Get mesh velocity for (int k=0; k<n_dims; k++) { temp_v(k)=(*grid_vel_fpts(j,i,k)); } }else{ temp_v.initialize_to_zero(); } // Interface unit-normal vector if (motion) { for (int m=0;m<n_dims;m++) norm(m) = *norm_dyn_fpts(j,i,m); }else{ for (int m=0;m<n_dims;m++) norm(m) = *norm_fpts(j,i,m); } // Calling Riemann solver if (run_input.riemann_solve_type==0) // Rusanov { // calculate flux from discontinuous solution at flux points if(n_dims==2) { calc_invf_2d(temp_u_l,temp_f_l); calc_invf_2d(temp_u_r,temp_f_r); if (motion) { calc_alef_2d(temp_u_l,temp_v,temp_f_l); calc_alef_2d(temp_u_r,temp_v,temp_f_r); } } else if(n_dims==3) { calc_invf_3d(temp_u_l,temp_f_l); calc_invf_3d(temp_u_r,temp_f_r); if (motion) { calc_alef_3d(temp_u_l,temp_v,temp_f_l); calc_alef_3d(temp_u_r,temp_v,temp_f_r); } } else FatalError("ERROR: Invalid number of dimensions ... "); rusanov_flux(temp_u_l,temp_u_r,temp_v,temp_f_l,temp_f_r,norm,fn,n_dims,n_fields,run_input.gamma); } else if (run_input.riemann_solve_type==1) { // Lax-Friedrich lax_friedrich(temp_u_l,temp_u_r,norm,fn,n_dims,n_fields,run_input.lambda,run_input.wave_speed); } else if (run_input.riemann_solve_type==2) { // ROE roe_flux(temp_u_l,temp_u_r,temp_v,norm,fn,n_dims,n_fields,run_input.gamma); } else FatalError("Riemann solver not implemented"); // Transform back to computational space from dynamic physical space if (motion) { for(int k=0; k<n_fields; k++) { (*norm_tconf_fpts_l(j,i,k)) = fn(k)*(*ndA_dyn_fpts_l(j,i))*(*tdA_fpts_l(j,i)); (*norm_tconf_fpts_r(j,i,k)) =-fn(k)*(*ndA_dyn_fpts_r(j,i))*(*tdA_fpts_r(j,i)); } } else { // Transform back to reference space from static physical space for(int k=0;k<n_fields;k++) { (*norm_tconf_fpts_l(j,i,k))= fn(k)*(*tdA_fpts_l(j,i)); (*norm_tconf_fpts_r(j,i,k))=-fn(k)*(*tdA_fpts_r(j,i)); } } if(viscous) { // Calling viscous riemann solver if (run_input.vis_riemann_solve_type==0) ldg_solution(0,temp_u_l,temp_u_r,u_c,run_input.pen_fact,norm); else FatalError("Viscous Riemann solver not implemented"); if (motion) // include transformation back to static space { for(int k=0;k<n_fields;k++) { *delta_disu_fpts_l(j,i,k) = (u_c(k) - temp_u_l(k))*(*J_dyn_fpts_l(j,i)); *delta_disu_fpts_r(j,i,k) = (u_c(k) - temp_u_r(k))*(*J_dyn_fpts_r(j,i)); } } else { for(int k=0;k<n_fields;k++) { *delta_disu_fpts_l(j,i,k) = (u_c(k) - temp_u_l(k)); *delta_disu_fpts_r(j,i,k) = (u_c(k) - temp_u_r(k)); } } } } } #endif } }
2.1.4 DONE viscous
- DONE calculategradient(indisuuptsfrom)
- 代码
void eles::calculate_gradient(int in_disu_upts_from) { if (n_eles!=0) { /*! Performs C = (alpha*A*B) + (beta*C) where: \n alpha = 1.0 \n beta = 0.0 \n A = opp_4 \n B = disu_upts \n C = grad_disu_upts opp_4 is the polynomial gradient matrix; has dimensions n_upts_per_ele by n_upts_per_ele Recall: opp_4(i)(k,j) = eval_d_nodal_basis(j,i,loc); = derivative of the jth nodal basis at the kth nodal (solution) point location in the reference domain for the ith dimension (vector of gradient values at solution points) = opp_4 * (vector of solution values at solution points in all elements of the same type) */ Arows = n_upts_per_ele; Acols = n_upts_per_ele; Brows = Acols; Bcols = n_fields*n_eles; Astride = Arows; Bstride = Brows; Cstride = Arows; #ifdef _CPU if(opp_4_sparse==0) // dense { #if defined _ACCELERATE_BLAS || defined _MKL_BLAS || defined _STANDARD_BLAS for (int i=0;i<n_dims;i++) { cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,Arows,Bcols,Acols,1.0,opp_4(i).get_ptr_cpu(),Astride,disu_upts(in_disu_upts_from).get_ptr_cpu(),Bstride,0.0,grad_disu_upts.get_ptr_cpu(0,0,0,i),Cstride); } #elif defined _NO_BLAS for (int i=0;i<n_dims;i++) { dgemm(Arows,Bcols,Acols,1.0,0.0,opp_4(i).get_ptr_cpu(),disu_upts(in_disu_upts_from).get_ptr_cpu(),grad_disu_upts.get_ptr_cpu(0,0,0,i)); } #endif } else if(opp_4_sparse==1) // mkl blas four-array csr format { #if defined _MKL_BLAS // implement #endif } else { cout << "ERROR: Unknown storage for opp_4 ... " << endl; } #endif #ifdef _GPU if (opp_4_sparse==0) { for (int i=0;i<n_dims;i++) { cublasDgemm('N','N',Arows,Bcols,Acols,1.0,opp_4(i).get_ptr_gpu(),Astride,disu_upts(in_disu_upts_from).get_ptr_gpu(),Bstride,0.0,grad_disu_upts.get_ptr_gpu(0,0,0,i),Cstride); } } else if (opp_4_sparse==1) { for (int i=0;i<n_dims;i++) { bespoke_SPMV(Arows,Acols,n_fields,n_eles,opp_4_ell_data(i).get_ptr_gpu(),opp_4_ell_indices(i).get_ptr_gpu(),opp_4_nnz_per_row(i),disu_upts(in_disu_upts_from).get_ptr_gpu(),grad_disu_upts.get_ptr_gpu(0,0,0,i),ele_type,order,0); } } #endif } /* cout << "OUTPUT" << endl; #ifdef _GPU grad_disu_upts.cp_gpu_cpu(); #endif for (int i=0;i<n_upts_per_ele;i++) for (int j=0;j<n_eles;j++) for (int k=0;k<n_fields;k++) for (int m=0;m<n_dims;m++) { if (ele2global_ele(j)==53) cout << setprecision(10) << i << " " << ele2global_ele(j) << " " << k << " " << m << " " << grad_disu_upts(i,j,k,m) << endl; } */ }
- 代码
- DONE correctgradient()
- 代码
void eles::correct_gradient(void) { if (n_eles!=0) { Arows = n_upts_per_ele; Acols = n_fpts_per_ele; Brows = Acols; Bcols = n_fields*n_eles; Astride = Arows; Bstride = Brows; Cstride = Arows; #ifdef _CPU if(opp_5_sparse==0) // dense { #if defined _ACCELERATE_BLAS || defined _MKL_BLAS || defined _STANDARD_BLAS for (int i=0;i<n_dims;i++) { cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,Arows,Bcols,Acols,1.0,opp_5(i).get_ptr_cpu(),Astride,delta_disu_fpts.get_ptr_cpu(),Bstride,1.0,grad_disu_upts.get_ptr_cpu(0,0,0,i),Cstride); } #elif defined _NO_BLAS for (int i=0;i<n_dims;i++) { dgemm(Arows,Bcols,Acols,1.0,1.0,opp_5(i).get_ptr_cpu(),delta_disu_fpts.get_ptr_cpu(),grad_disu_upts.get_ptr_cpu(0,0,0,i)); } #endif } else if(opp_5_sparse==1) // mkl blas four-array csr format { #if defined _MKL_BLAS // impelement #endif } else { cout << "ERROR: Unknown storage for opp_5 ... " << endl; } // Transform to physical space double detjac; double inv_detjac; double rx,ry,rz,sx,sy,sz,tx,ty,tz; double Xx,Xy,Xz,Yx,Yy,Yz,Zx,Zy,Zz; double ur,us,ut,uX,uY,uZ; for (int i=0;i<n_eles;i++) { for (int j=0;j<n_upts_per_ele;j++) { // Transform to static-physical domain detjac = detjac_upts(j,i); inv_detjac = 1.0/detjac; rx = JGinv_upts(0,0,j,i)*inv_detjac; ry = JGinv_upts(0,1,j,i)*inv_detjac; sx = JGinv_upts(1,0,j,i)*inv_detjac; sy = JGinv_upts(1,1,j,i)*inv_detjac; //physical gradient if(n_dims==2) { for(int k=0;k<n_fields;k++) { ur = grad_disu_upts(j,i,k,0); us = grad_disu_upts(j,i,k,1); grad_disu_upts(j,i,k,0) = ur*rx + us*sx; grad_disu_upts(j,i,k,1) = ur*ry + us*sy; } } if (n_dims==3) { rz = JGinv_upts(0,2,j,i)*inv_detjac; sz = JGinv_upts(1,2,j,i)*inv_detjac; tx = JGinv_upts(2,0,j,i)*inv_detjac; ty = JGinv_upts(2,1,j,i)*inv_detjac; tz = JGinv_upts(2,2,j,i)*inv_detjac; for (int k=0;k<n_fields;k++) { ur = grad_disu_upts(j,i,k,0); us = grad_disu_upts(j,i,k,1); ut = grad_disu_upts(j,i,k,2); grad_disu_upts(j,i,k,0) = ur*rx + us*sx + ut*tx; grad_disu_upts(j,i,k,1) = ur*ry + us*sy + ut*ty; grad_disu_upts(j,i,k,2) = ur*rz + us*sz + ut*tz; } } if (motion) { // Transform to dynamic-physical domain detjac = J_dyn_upts(j,i); inv_detjac = 1.0/detjac; Xx = JGinv_dyn_upts(0,0,j,i)*inv_detjac; Xy = JGinv_dyn_upts(0,1,j,i)*inv_detjac; Yx = JGinv_dyn_upts(1,1,j,i)*inv_detjac; Yy = JGinv_dyn_upts(1,1,j,i)*inv_detjac; //physical gradient if(n_dims==2) { for(int k=0;k<n_fields;k++) { uX = grad_disu_upts(j,i,k,0); uY = grad_disu_upts(j,i,k,1); grad_disu_upts(j,i,k,0) = uX*Xx + uY*Yx; grad_disu_upts(j,i,k,1) = uX*Xy + uY*Yy; } } if (n_dims==3) { Xz = JGinv_dyn_upts(j,i,0,2)*inv_detjac; Yz = JGinv_dyn_upts(j,i,1,2)*inv_detjac; Zx = JGinv_dyn_upts(j,i,2,0)*inv_detjac; Zy = JGinv_dyn_upts(j,i,2,1)*inv_detjac; Zz = JGinv_dyn_upts(j,i,2,2)*inv_detjac; for (int k=0;k<n_fields;k++) { uX = grad_disu_upts(j,i,k,0); uY = grad_disu_upts(j,i,k,1); uZ = grad_disu_upts(j,i,k,2); grad_disu_upts(j,i,k,0) = uX*Xx + uY*Yx + uZ*Zx; grad_disu_upts(j,i,k,1) = uX*Xy + uY*Yy + uZ*Zy; grad_disu_upts(j,i,k,2) = uX*Xz + uY*Yz + uZ*Zz; } } } } } #endif }
- 代码
- DONE extrapolatecorrectedgradient()
- 代码
void eles::extrapolate_corrected_gradient(void) { if (n_eles!=0) { Arows = n_fpts_per_ele; Acols = n_upts_per_ele; Brows = Acols; Bcols = n_fields*n_eles; Astride = Arows; Bstride = Brows; Cstride = Arows; #ifdef _CPU if(opp_6_sparse==0) // dense { #if defined _ACCELERATE_BLAS || defined _MKL_BLAS || defined _STANDARD_BLAS for (int i=0;i<n_dims;i++) { cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,Arows,Bcols,Acols,1.0,opp_6.get_ptr_cpu(),Astride,grad_disu_upts.get_ptr_cpu(0,0,0,i),Bstride,0.0,grad_disu_fpts.get_ptr_cpu(0,0,0,i),Cstride); } #elif defined _NO_BLAS for (int i=0;i<n_dims;i++) { dgemm(Arows,Bcols,Acols,1.0,0.0,opp_6.get_ptr_cpu(),grad_disu_upts.get_ptr_cpu(0,0,0,i),grad_disu_fpts.get_ptr_cpu(0,0,0,i)); } #endif } else if(opp_6_sparse==1) // mkl blas four-array csr format { #if defined _MKL_BLAS // implement #endif } else { cout << "ERROR: Unknown storage for opp_6 ... " << endl; } #endif }
- 代码
- DONE evaluateviscFlux(indisuuptsfrom)
- 代码
void eles::evaluate_viscFlux(int in_disu_upts_from) { if (n_eles!=0) { #ifdef _CPU int i,j,k,l,m; double detjac; for(i=0;i<n_eles;i++) { // Calculate viscous flux for(j=0;j<n_upts_per_ele;j++) { detjac = detjac_upts(j,i); // solution in static-physical domain for(k=0;k<n_fields;k++) { temp_u(k)=disu_upts(in_disu_upts_from)(j,i,k); // gradient in dynamic-physical domain for (m=0;m<n_dims;m++) { temp_grad_u(k,m) = grad_disu_upts(j,i,k,m); } } // Transform to dynamic-physical domain if (motion) { for (k=0; k<n_fields; k++) { temp_u(k) /= J_dyn_upts(j,i); } } if(n_dims==2) { calc_visf_2d(temp_u,temp_grad_u,temp_f); } else if(n_dims==3) { calc_visf_3d(temp_u,temp_grad_u,temp_f); } else { cout << "ERROR: Invalid number of dimensions ... " << endl; } // If LES or wall model, calculate SGS viscous flux if(LES != 0 || wall_model != 0) { calc_sgsf_upts(temp_u,temp_grad_u,detjac,i,j,temp_sgsf); // Add SGS or wall flux to viscous flux for(k=0;k<n_fields;k++) for(l=0;l<n_dims;l++) temp_f(k,l) += temp_sgsf(k,l); } // If LES, add SGS flux to global array (needed for interface flux calc) if(LES > 0) { // Transfer back to static-phsycial domain if (motion) { temp_sgsf_ref.initialize_to_zero(); for(k=0;k<n_fields;k++) { for(l=0;l<n_dims;l++) { for(m=0;m<n_dims;m++) { temp_sgsf_ref(k,l)+=JGinv_dyn_upts(l,m,j,i)*temp_sgsf(k,m); } } } // Copy back to original flux array for (k=0; k<n_fields; k++) { for(l=0; l<n_dims; l++) { temp_sgsf(k,l) = temp_sgsf_ref(k,l); } } } // Transfer back to computational domain for(k=0;k<n_fields;k++) { for(l=0;l<n_dims;l++) { sgsf_upts(j,i,k,l) = 0.0; for(m=0;m<n_dims;m++) { sgsf_upts(j,i,k,l)+=JGinv_upts(l,m,j,i)*temp_sgsf(k,m); } } } } // Transfer back to static-phsycial domain if (motion) { temp_f_ref.initialize_to_zero(); for(k=0;k<n_fields;k++) { for(l=0;l<n_dims;l++) { for(m=0;m<n_dims;m++) { temp_f_ref(k,l)+=JGinv_dyn_upts(l,m,j,i)*temp_f(k,m); } } } // Copy back to original flux array for(l=0; l<n_dims; l++) { for (k=0; k<n_fields; k++) { temp_f(k,l) = temp_f_ref(k,l); } } } // Transform viscous flux for(k=0;k<n_fields;k++) { for(l=0;l<n_dims;l++) { for(m=0;m<n_dims;m++) { tdisf_upts(j,i,k,l)+=JGinv_upts(l,m,j,i)*temp_f(k,m); } } } } } #endif } }
- 代码
2.1.6 DONE calculatedivergence(indivtconfuptsto)
- 代码
void eles::calculate_divergence(int in_div_tconf_upts_to) { if (n_eles!=0) { #ifdef _CPU if(opp_2_sparse==0) // dense { #if defined _ACCELERATE_BLAS || defined _MKL_BLAS || defined _STANDARD_BLAS cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n_upts_per_ele,n_fields*n_eles,n_upts_per_ele,1.0,opp_2(0).get_ptr_cpu(),n_upts_per_ele,tdisf_upts.get_ptr_cpu(0,0,0,0),n_upts_per_ele,0.0,div_tconf_upts(in_div_tconf_upts_to).get_ptr_cpu(),n_upts_per_ele); for (int i=1;i<n_dims;i++) { cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n_upts_per_ele,n_fields*n_eles,n_upts_per_ele,1.0,opp_2(i).get_ptr_cpu(),n_upts_per_ele,tdisf_upts.get_ptr_cpu(0,0,0,i),n_upts_per_ele,1.0,div_tconf_upts(in_div_tconf_upts_to).get_ptr_cpu(),n_upts_per_ele); } #elif defined _NO_BLAS dgemm(n_upts_per_ele,n_fields*n_eles,n_upts_per_ele,1.0,0.0,opp_2(0).get_ptr_cpu(),tdisf_upts.get_ptr_cpu(0,0,0,0),div_tconf_upts(in_div_tconf_upts_to).get_ptr_cpu()); for (int i=1;i<n_dims;i++) { dgemm(n_upts_per_ele,n_fields*n_eles,n_upts_per_ele,1.0,1.0,opp_2(i).get_ptr_cpu(),tdisf_upts.get_ptr_cpu(0,0,0,i),div_tconf_upts(in_div_tconf_upts_to).get_ptr_cpu()); } #endif } else if(opp_2_sparse==1) // mkl blas four-array csr format { #if defined _MKL_BLAS mkl_dcsrmm(&transa,&n_upts_per_ele,&n_fields_mul_n_eles,&n_upts_per_ele,&one,matdescra,opp_2_data(0).get_ptr_cpu(),opp_2_cols(0).get_ptr_cpu(),opp_2_b(0).get_ptr_cpu(),opp_2_e(0).get_ptr_cpu(),tdisf_upts.get_ptr_cpu(0,0,0,0),&n_upts_per_ele,&zero,div_tconf_upts(in_div_tconf_upts_to).get_ptr_cpu(),&n_upts_per_ele); for (int i=1;i<n_dims;i++) { mkl_dcsrmm(&transa,&n_upts_per_ele,&n_fields_mul_n_eles,&n_upts_per_ele,&one,matdescra,opp_2_data(i).get_ptr_cpu(),opp_2_cols(i).get_ptr_cpu(),opp_2_b(i).get_ptr_cpu(),opp_2_e(i).get_ptr_cpu(),tdisf_upts.get_ptr_cpu(0,0,0,i),&n_upts_per_ele,&one,div_tconf_upts(in_div_tconf_upts_to).get_ptr_cpu(),&n_upts_per_ele); } #endif } else { cout << "ERROR: Unknown storage for opp_2 ... " << endl; } } /* for (int j=0;j<n_eles;j++) for (int i=0;i<n_upts_per_ele;i++) //for (int k=0;k<n_fields;k++) cout << scientific << setw(16) << setprecision(12) << div_tconf_upts(0)(i,j,0) << endl; */ }
2.1.7 HOLD viscous
- DONE calculatecommonviscFlux()
- 代码
void int_inters::calculate_common_viscFlux(void) { #ifdef _CPU array<double> norm(n_dims), fn(n_fields); for(int i=0;i<n_inters;i++) { for(int j=0;j<n_fpts_per_inter;j++) { // obtain discontinuous solution at flux points if (motion) { // Transform to dynamic-physical domain for(int k=0;k<n_fields;k++) { temp_u_l(k)=(*disu_fpts_l(j,i,k))/(*J_dyn_fpts_l(j,i)); temp_u_r(k)=(*disu_fpts_r(j,i,k))/(*J_dyn_fpts_r(j,i)); } } else { for(int k=0;k<n_fields;k++) { temp_u_l(k)=(*disu_fpts_l(j,i,k)); temp_u_r(k)=(*disu_fpts_r(j,i,k)); } } // obtain physical gradient of discontinuous solution at flux points for(int k=0;k<n_dims;k++) { for(int l=0;l<n_fields;l++) { temp_grad_u_l(l,k) = *grad_disu_fpts_l(j,i,l,k); temp_grad_u_r(l,k) = *grad_disu_fpts_r(j,i,l,k); } } // calculate flux from discontinuous solution at flux points if(n_dims==2) { calc_visf_2d(temp_u_l,temp_grad_u_l,temp_f_l); calc_visf_2d(temp_u_r,temp_grad_u_r,temp_f_r); } else if(n_dims==3) { calc_visf_3d(temp_u_l,temp_grad_u_l,temp_f_l); calc_visf_3d(temp_u_r,temp_grad_u_r,temp_f_r); } else FatalError("ERROR: Invalid number of dimensions ... "); // If LES, get SGS flux and add to viscous flux if(LES) { for(int k=0;k<n_dims;k++) { for(int l=0;l<n_fields;l++) { // pointers to subgrid-scale fluxes temp_sgsf_l(l,k) = *sgsf_fpts_l(j,i,l,k); temp_sgsf_r(l,k) = *sgsf_fpts_r(j,i,l,k); // Add SGS fluxes to viscous fluxes temp_f_l(l,k) += temp_sgsf_l(l,k); temp_f_r(l,k) += temp_sgsf_r(l,k); } } } // storing normal components if (motion) { for (int m=0;m<n_dims;m++) norm(m) = *norm_dyn_fpts(j,i,m); } else { for (int m=0;m<n_dims;m++) norm(m) = *norm_fpts(j,i,m); } // Calling viscous riemann solver if (run_input.vis_riemann_solve_type==0) ldg_flux(0,temp_u_l,temp_u_r,temp_f_l,temp_f_r,norm,fn,n_dims,n_fields,run_input.tau,run_input.pen_fact); else FatalError("Viscous Riemann solver not implemented"); // Transform back to reference space if (motion) { for(int k=0;k<n_fields;k++) { (*norm_tconf_fpts_l(j,i,k))+= fn(k)*(*tdA_fpts_l(j,i))*(*ndA_dyn_fpts_l(j,i)); (*norm_tconf_fpts_r(j,i,k))+= -fn(k)*(*tdA_fpts_r(j,i))*(*ndA_dyn_fpts_r(j,i)); } } else { for(int k=0;k<n_fields;k++) { (*norm_tconf_fpts_l(j,i,k))+= fn(k)*(*tdA_fpts_l(j,i)); (*norm_tconf_fpts_r(j,i,k))+= -fn(k)*(*tdA_fpts_r(j,i)); } } } } #endif }
- TODO evaluateboundaryConditionsviscFlux(FlowSol->time)
2.1.8 DONE calculatecorrecteddivergence(indivtconfuptsto)
- 代码
void eles::calculate_corrected_divergence(int in_div_tconf_upts_to) { if (n_eles!=0) { #ifdef _CPU #if defined _ACCELERATE_BLAS || defined _MKL_BLAS || defined _STANDARD_BLAS cblas_daxpy(n_eles*n_fields*n_fpts_per_ele,-1.0,norm_tdisf_fpts.get_ptr_cpu(),1,norm_tconf_fpts.get_ptr_cpu(),1); #elif defined _NO_BLAS daxpy(n_eles*n_fields*n_fpts_per_ele,-1.0,norm_tdisf_fpts.get_ptr_cpu(),norm_tconf_fpts.get_ptr_cpu()); #endif if(opp_3_sparse==0) // dense { #if defined _ACCELERATE_BLAS || defined _MKL_BLAS || defined _STANDARD_BLAS cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n_upts_per_ele,n_fields*n_eles,n_fpts_per_ele,1.0,opp_3.get_ptr_cpu(),n_upts_per_ele,norm_tconf_fpts.get_ptr_cpu(),n_fpts_per_ele,1.0,div_tconf_upts(in_div_tconf_upts_to).get_ptr_cpu(),n_upts_per_ele); #elif defined _NO_BLAS dgemm(n_upts_per_ele,n_fields*n_eles,n_fpts_per_ele,1.0,1.0,opp_3.get_ptr_cpu(),norm_tconf_fpts.get_ptr_cpu(),div_tconf_upts(in_div_tconf_upts_to).get_ptr_cpu()); #endif } else if(opp_3_sparse==1) // mkl blas four-array csr format { #if defined _MKL_BLAS mkl_dcsrmm(&transa,&n_upts_per_ele,&n_fields_mul_n_eles,&n_fpts_per_ele,&one,matdescra,opp_3_data.get_ptr_cpu(),opp_3_cols.get_ptr_cpu(),opp_3_b.get_ptr_cpu(),opp_3_e.get_ptr_cpu(),norm_tconf_fpts.get_ptr_cpu(),&n_fpts_per_ele,&one,div_tconf_upts(in_div_tconf_upts_to).get_ptr_cpu(),&n_upts_per_ele); #endif } else { cout << "ERROR: Unknown storage for opp_3 ... " << endl; } for (int i=0;i<n_upts_per_ele;i++) for (int j=0;j<n_eles;j++) for (int k=0;k<n_fields;k++) if (isnan(div_tconf_upts(in_div_tconf_upts_to)(j,i,k))) FatalError("NaN in residual, exiting."); #endif } }
3 程序变动笔记
3.1 DONE 第一次更改程序
4 问题[0/3]
[ ]eles计算程序最终从mesh网格中获取了什么？[ ]Array数组是如何实现的？[ ]搞明白那些是再reference Space计算的？
扫一扫
234
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
