#include "site-stuff.h" #include #include "MHD.h" #include "variables.h" #ifdef LC_UNSC #define IBOUND ibound_ #define JBOUND jbound_ #define KBOUND kbound_ #define METRIC metric_ #define INITAL inital_ #define MHD_FORTRAN_ARRAYS mhd_fortran_arrays_ #define ZERO zero_ #define MHD_HACK mhd_hack_ #define FLUX_FIX flux_fix_ #define FDIVCHECK fdivcheck_ #define RINGAV ringav_ #define HDFTAKE_CLOSE hdftake_close_ #elif LC_NOUNSC #define IBOUND ibound #define JBOUND jbound #define KBOUND kbound #define METRIC metric #define INITAL inital #define MHD_FORTRAN_ARRAYS mhd_fortran_arrays #define ZERO zero #define MHD_HACK mhd_hack #define FLUX_FIX flux_fix #define FDIVCHECK fdivcheck #define RINGAV ringav #define HDFTAKE_CLOSE hdftake_close #endif #define PI 3.1415926 #define REARTH 6.5e8 /// Defined in metric.F extern "C" void METRIC(); /// Defined in init-fortran.F extern "C" void INITAL(); /// Defined in hack.F extern "C" void MHD_HACK(); /// Defined in hack.F extern "C" void FLUX_FIX(FLoaT *); /// Defined in test-arrays.F extern "C" void FDIVCHECK(int *, int *, int*); /// Defined in therest.F extern "C" void RINGAV(); /// Defined in bound.F extern "C" void IBOUND(FLoaT *, FLoaT*, int *, int *, int *); /// Defined in bound.F extern "C" void JBOUND(FLoaT *, FLoaT*, int *, int *, int *); /// Defined in bound.F extern "C" void KBOUND(FLoaT *, FLoaT*, int *, int *, int *); /// Defined in test-arrays.F extern "C" void TEST_ARRAYS(); /// Defined in zero.F extern "C" void ZERO(FLoaT *); /// Defined in mhd_fortran_arrays.F extern "C" void MHD_FORTRAN_ARRAYS(FLoaT **, int*, int*, int*, int*, int*); /// Defined in hdftake-para.F extern "C" int HDFTAKE_CLOSE(); void get_local_fortran_info(int , int, int, int, int , int * ); /// Constructor - Initialize the arrays /** * "initial" does most of the dirty work * * @param ni * @param nj * @param nk * @param no */ MHD::MHD(int ni, int nj, int nk, int no) :X(ni+1,nj+1,nk+1), Y(ni+1,nj+1,nk+1),Z(ni+1,nj+1,nk+1), Volume(ni+1,nj+1,nk+1),Volume_inv(ni+1,nj+1,nk+1), Face_i(ni+1,nj+1,nk+1),Face_j(ni+1,nj+1,nk+1),Face_k(ni+1,nj+1,nk+1), Tran_i(ni+1,nj+1,nk+1,9),Tran_j(ni+1,nj+1,nk+1,9),Tran_k(ni+1,nj+1,nk+1,9), Edge_i(ni+1,nj+1,nk+1),Edge_j(ni+1,nj+1,nk+1),Edge_k(ni+1,nj+1,nk+1), Vec_i(ni+1,nj+1,nk+1,9),Vec_j(ni+1,nj+1,nk+1,9),Vec_k(ni+1,nj+1,nk+1,9), Cell_Length(ni+1,nj+1,nk+1), Rho(ni+1,nj+1,nk+1),Vx(ni+1,nj+1,nk+1),Vy(ni+1,nj+1,nk+1), Vz(ni+1,nj+1,nk+1),Csound(ni+1,nj+1,nk+1),Pressure(ni+1,nj+1,nk+1), Bx(ni+1,nj+1,nk+1),By(ni+1,nj+1,nk+1),Bz(ni+1,nj+1,nk+1), Bi(ni+1,nj+1,nk+1),Bj(ni+1,nj+1,nk+1),Bk(ni+1,nj+1,nk+1), Ei(ni+1,nj+1,nk+1),Ej(ni+1,nj+1,nk+1),Ek(ni+1,nj+1,nk+1), Rho_1(ni+1,nj+1,nk+1),Vx_1(ni+1,nj+1,nk+1),Vy_1(ni+1,nj+1,nk+1), Vz_1(ni+1,nj+1,nk+1),Csound_1(ni+1,nj+1,nk+1),Pressure_1(ni+1,nj+1,nk+1), Bx_1(ni+1,nj+1,nk+1),By_1(ni+1,nj+1,nk+1),Bz_1(ni+1,nj+1,nk+1), Bi_1(ni+1,nj+1,nk+1),Bj_1(ni+1,nj+1,nk+1),Bk_1(ni+1,nj+1,nk+1) { #ifdef DEBUG_MODE_ON printf("DEBUG: process (%d) in MHD.C::MHD::MHD(...)\n", RTIME.local_proc); #endif // set basic dimensioning information NI = ni; NJ = nj; NK = nk; NORDER=no; NO2 = no/2; // setBase calls will give the following effective Indexes Range I(-NO2+1,NI+12); // Fortran indexing (1 - NI) active Range J(-NO2+1,NJ+12); Range K(-NO2+1,NK+12); // reset Bases X.setBase(1); Y.setBase(1); Z.setBase(1); Volume.setBase(1); Volume_inv.setBase(1); Face_i.setBase(1);Face_j.setBase(1);Face_k.setBase(1); Tran_i.setBase(1);Tran_j.setBase(1);Tran_k.setBase(1); Edge_i.setBase(1);Edge_j.setBase(1);Edge_k.setBase(1); Vec_i.setBase(1);Vec_j.setBase(1);Vec_k.setBase(1); Cell_Length.setBase(1); // X2.setBase(1); Y2.setBase(1); Phi.setBase(1); Rho.setBase(1); Vx.setBase(1); Vy.setBase(1); Vz.setBase(1); Csound.setBase(1); Pressure.setBase(1); Bx.setBase(1); By.setBase(1); Bz.setBase(1); Bi.setBase(1); Bj.setBase(1); Bk.setBase(1); Ei.setBase(1); Ej.setBase(1); Ek.setBase(1); Rho_1.setBase(1); Vx_1.setBase(1); Vy_1.setBase(1); Vz_1.setBase(1); Csound_1.setBase(1); Pressure_1.setBase(1); Bx_1.setBase(1); By_1.setBase(1); Bz_1.setBase(1); Bi_1.setBase(1); Bj_1.setBase(1); Bk_1.setBase(1); } /***************************************************************************/ /// Constructor - Initialize the arrays /** * constructor for parallel arrays spread over proc_first to proc_last * @param ni * @param nj * @param nk * @param no * @param mhd_procs - number of processors */ MHD::MHD(int ni, int nj, int nk, int no, Range mhd_procs) :X(ni+1,nj+1,nk+1), Y(ni+1,nj+1,nk+1),Z(ni+1,nj+1,nk+1), Volume(ni+1,nj+1,nk+1),Volume_inv(ni+1,nj+1,nk+1), Face_i(ni+1,nj+1,nk+1),Face_j(ni+1,nj+1,nk+1),Face_k(ni+1,nj+1,nk+1), Tran_i(ni+1,nj+1,nk+1,9),Tran_j(ni+1,nj+1,nk+1,9),Tran_k(ni+1,nj+1,nk+1,9), Edge_i(ni+1,nj+1,nk+1),Edge_j(ni+1,nj+1,nk+1),Edge_k(ni+1,nj+1,nk+1), Vec_i(ni+1,nj+1,nk+1,9),Vec_j(ni+1,nj+1,nk+1,9),Vec_k(ni+1,nj+1,nk+1,9), Cell_Length(ni+1,nj+1,nk+1), Rho(ni+1,nj+1,nk+1),Vx(ni+1,nj+1,nk+1),Vy(ni+1,nj+1,nk+1), Vz(ni+1,nj+1,nk+1),Csound(ni+1,nj+1,nk+1),Pressure(ni+1,nj+1,nk+1), Bx(ni+1,nj+1,nk+1),By(ni+1,nj+1,nk+1),Bz(ni+1,nj+1,nk+1), Bi(ni+1,nj+1,nk+1),Bj(ni+1,nj+1,nk+1),Bk(ni+1,nj+1,nk+1), Ei(ni+1,nj+1,nk+1),Ej(ni+1,nj+1,nk+1),Ek(ni+1,nj+1,nk+1), Rho_1(ni+1,nj+1,nk+1),Vx_1(ni+1,nj+1,nk+1),Vy_1(ni+1,nj+1,nk+1), Vz_1(ni+1,nj+1,nk+1),Csound_1(ni+1,nj+1,nk+1),Pressure_1(ni+1,nj+1,nk+1), Bx_1(ni+1,nj+1,nk+1),By_1(ni+1,nj+1,nk+1),Bz_1(ni+1,nj+1,nk+1), Bi_1(ni+1,nj+1,nk+1),Bj_1(ni+1,nj+1,nk+1),Bk_1(ni+1,nj+1,nk+1), MHD_P(mhd_procs), Tran_P(mhd_procs) { #ifdef DEBUG_MODE_ON printf("DEBUG: process (%d) in MHD.C::MHD::MHD(...)\n", RTIME.local_proc); #endif // set basic dimensioning information #ifdef MSETUP /********* * FIXME * ********************************************************************** * In an email on January 21, 2008, Slava said: * * Some of the output (in MSETUP mode) is supposed to be used by * create-param.pl script to read in some dimensioning * information. The printf("inside the parallel code...") and * printf("On processor...") in MHD.C are absolutely necessary for * create-param.pl. You can see these by svn diffing revision 221 to * revision 217 of MHD.C. Perhaps, some other output statements are * important as well, I have not checked very thoroughly. * * This output should be written better to make more sense... * * Best solution: Eliminate create-param.pl code and use MSPHERE.C to * handle include file setup. **********************************************************************/ std::cerr << "***MHD.C*** To-Do: MSETUP dimension output for create-param.pl\n"; std::cerr << "Should be a cleaner way of outputting dimension info for MSETUP.\n"; std::cout << "inside the parallel code " << ni << " " << nj << " " << nk << " " << no << std::endl; #endif #ifdef DEBUG_MODE_ON std::cout << "MHD.C: Dimensions = <" << ni << ", " << nj << ", " << nk << ", " << no << ">\n"; #endif NI = ni; NJ = nj; NK = nk; NORDER=no; NO2 = no/2; // partioning definition for the MHD arrays // set k partitioning to false to simplify ringav /* Partitioning_Type MHD_P(part);*/ MHD_P.partitionAlongAxis(0,TRUE,NO2); MHD_P.partitionAlongAxis(1,TRUE,NO2); MHD_P.partitionAlongAxis(2,FALSE,NO2); /* Partitioning_Type Tran_P(part); */ Tran_P.partitionAlongAxis(0,TRUE,NO2); Tran_P.partitionAlongAxis(1,TRUE,NO2); Tran_P.partitionAlongAxis(2,FALSE,NO2); Tran_P.partitionAlongAxis(3,FALSE,0); // setBase calls will give the following effective Indexes Range I(1,NI+1); // Fortran indexing (1 - NI) active Range J(1,NJ+1); Range K(1,NK+1); // reset Bases X.setBase(1); Y.setBase(1); Z.setBase(1); Volume.setBase(1); Volume_inv.setBase(1); Face_i.setBase(1);Face_j.setBase(1);Face_k.setBase(1); Tran_i.setBase(1);Tran_j.setBase(1);Tran_k.setBase(1); Edge_i.setBase(1);Edge_j.setBase(1);Edge_k.setBase(1); Vec_i.setBase(1);Vec_j.setBase(1);Vec_k.setBase(1); Cell_Length.setBase(1); Rho.setBase(1); Vx.setBase(1); Vy.setBase(1); Vz.setBase(1); Csound.setBase(1); Pressure.setBase(1); Bx.setBase(1); By.setBase(1); Bz.setBase(1); Bi.setBase(1); Bj.setBase(1); Bk.setBase(1); Ei.setBase(1); Ej.setBase(1); Ek.setBase(1); Rho_1.setBase(1); Vx_1.setBase(1); Vy_1.setBase(1); Vz_1.setBase(1); Csound_1.setBase(1); Pressure_1.setBase(1); Bx_1.setBase(1); By_1.setBase(1); Bz_1.setBase(1); Bi_1.setBase(1); Bj_1.setBase(1); Bk_1.setBase(1); // partition the arrays throughout the processors X.partition(MHD_P); // MHD_P.display("MHD_P Partition"); #ifdef MSETUP X.displayPartitioning("wahoo Mountain Dew"); #endif Y.partition(MHD_P); Z.partition(MHD_P); Volume.partition(MHD_P); Volume_inv.partition(MHD_P); Face_i.partition(MHD_P);Face_j.partition(MHD_P);Face_k.partition(MHD_P); Tran_i.partition(Tran_P);Tran_j.partition(Tran_P);Tran_k.partition(Tran_P); Edge_i.partition(MHD_P);Edge_j.partition(MHD_P);Edge_k.partition(MHD_P); Vec_i.partition(Tran_P);Vec_j.partition(Tran_P);Vec_k.partition(Tran_P); Cell_Length.partition(MHD_P); Rho.partition(MHD_P); Vx.partition(MHD_P); Vy.partition(MHD_P); Vz.partition(MHD_P); Csound.partition(MHD_P); Pressure.partition(MHD_P); Bx.partition(MHD_P); By.partition(MHD_P); Bz.partition(MHD_P); Bi.partition(MHD_P); Bj.partition(MHD_P); Bk.partition(MHD_P); Ei.partition(MHD_P); Ej.partition(MHD_P); Ek.partition(MHD_P); Rho_1.partition(MHD_P); Vx_1.partition(MHD_P); Vy_1.partition(MHD_P); Vz_1.partition(MHD_P); Csound_1.partition(MHD_P); Pressure_1.partition(MHD_P); Bx_1.partition(MHD_P); By_1.partition(MHD_P); Bz_1.partition(MHD_P); Bi_1.partition(MHD_P); Bj_1.partition(MHD_P); Bk_1.partition(MHD_P); // Tran_P.display("Tran_P partition"); Tran_i.displayPartitioning("tran stuff"); #ifdef MSETUP MHD_P.display("MHD_P again"); #endif // done with initialization int klow = X.getLocalBase(2); int khigh = X.getLocalBound(2); int jlow = X.getLocalBase(1); int jhigh = X.getLocalBound(1); int ilow = X.getLocalBase(0); int ihigh = X.getLocalBound(0); #ifdef MSETUP /********* * FIXME * ********************************************************************** * In an email on January 21, 2008, Slava said: * * Some of the output (in MSETUP mode) is supposed to be used by * create-param.pl script to read in some dimensioning * information. The printf("inside the parallel code...") and * printf("On processor...") in MHD.C are absolutely necessary for * create-param.pl. You can see these by svn diffing revision 221 to * revision 217 of MHD.C. Perhaps, some other output statements are * important as well, I have not checked very thoroughly. * * This output should be written better to make more sense... * * Best solution: Eliminate create-param.pl code and use MSPHERE.C to * handle include file setup. **********************************************************************/ std::cerr << "***MHD.C*** FIXME: MSETUP dimension output for create-param.pl"; std::cerr << "Should be a cleaner way of outputting dimension info for MSETUP.\n"; int my_pe = Communication_Manager::My_Process_Number; printf("On processor %d :: %d %d %d %d %d %d\n",my_pe,ilow,ihigh, jlow,jhigh,klow,khigh); #endif #ifdef DEBUG_MODE_ON std:: cout << "MHD.C: On processor " << my_pe << ": \n"; std::cout << " = <" << ilow <<", " << ihigh << ">\n"; std::cout << " = <" << jlow <<", " << jhigh << ">\n"; std::cout << " = <" << klow <<", " << khigh << ">\n"; #endif } /* MHD::MHD(int ni, int nj, int nk, int no, Range mhd_procs) */ /***************************************************************************/ #ifndef MSETUP /// Initializes the MHD data for MSETUP /** * The constructors * - MHD(int, int, int, int ); * - MHD(int,int,int,int,Range); * initialize the (MHD data) arrays, but most of the dirty work is * actually done in MHD::initial(int,int,int). * * initial hands off the location of the local data arrays to fortran. * * @note initial is only defined with the preprocessor directive MSETUP. * * @fixme the "initial" function should probably be a subclass of MHD * rather than a method of MHD. For example, we could have an MSETUP * sub-class. This may reduce some of the confusing preprocessor * directives in MHD.C. */ void MHD::initial(int nmhd, int proc_first, int proc_last) { #ifdef DEBUG_MODE_ON printf("DEBUG: process (%d) in MHD.C::MHD::initial(...)\n", RTIME.local_proc); #endif FLoaT *variable[100]; int i_info[6], j_info[6], k_info[6]; int mhd_proc; // use variables.h to get the indices for the arrays variable[ _x] = X.getDataPointer(); variable[ _y] = Y.getDataPointer(); variable[ _z] = Z.getDataPointer(); variable[ _vol] = Volume.getDataPointer(); variable[ _volinv] = Volume_inv.getDataPointer(); variable[ _face_i] = Face_i.getDataPointer(); variable[ _tran_i] = Tran_i.getDataPointer(); variable[ _face_j] = Face_j.getDataPointer(); variable[ _tran_j] = Tran_j.getDataPointer(); variable[ _face_k] = Face_k.getDataPointer(); variable[ _tran_k] = Tran_k.getDataPointer(); variable[ _edge_i] = Edge_i.getDataPointer(); variable[ _vec_i] = Vec_i.getDataPointer(); variable[ _edge_j] = Edge_j.getDataPointer(); variable[ _vec_j] = Vec_j.getDataPointer(); variable[ _edge_k] = Edge_k.getDataPointer(); variable[ _vec_k] = Vec_k.getDataPointer(); variable[ _cell_length] = Cell_Length.getDataPointer(); variable[ _rho] = Rho.getDataPointer(); variable[ _vx] = Vx.getDataPointer(); variable[ _vy] = Vy.getDataPointer(); variable[ _vz] = Vz.getDataPointer(); variable[ _c] = Csound.getDataPointer(); variable[ _p] = Pressure.getDataPointer(); variable[ _bx] = Bx.getDataPointer(); variable[ _by] = By.getDataPointer(); variable[ _bz] = Bz.getDataPointer(); variable[ _bi] = Bi.getDataPointer(); variable[ _bj] = Bj.getDataPointer(); variable[ _bk] = Bk.getDataPointer(); variable[ _ei] = Ei.getDataPointer(); variable[ _ej] = Ej.getDataPointer(); variable[ _ek] = Ek.getDataPointer(); variable[ _rho1] = Rho_1.getDataPointer(); variable[ _vx1] = Vx_1.getDataPointer(); variable[ _vy1] = Vy_1.getDataPointer(); variable[ _vz1] = Vz_1.getDataPointer(); variable[ _c1] = Csound_1.getDataPointer(); variable[ _p1] = Pressure_1.getDataPointer(); variable[ _bx1] = Bx_1.getDataPointer(); variable[ _by1] = By_1.getDataPointer(); variable[ _bz1] = Bz_1.getDataPointer(); variable[ _bi1] = Bi_1.getDataPointer(); variable[ _bj1] = Bj_1.getDataPointer(); variable[ _bk1] = Bk_1.getDataPointer(); /* * The fortran call passes off the pointers to the data in the MHD * class through the "variable" pointer array. The local data array * range is passed off in the info arrays which contain [0] -> * local lower, [1] local high dimension, [2] -> local active high * dimension, [3] offset from local index to global index * */ int my_pe = Communication_Manager::My_Process_Number; get_local_fortran_info(X.getLocalBase(0), X.getLocalBound(0), X.getBase(0), X.getBound(0), NO2, i_info); get_local_fortran_info(X.getLocalBase(1), X.getLocalBound(1), X.getBase(1), X.getBound(1), NO2, j_info); get_local_fortran_info(X.getLocalBase(2), X.getLocalBound(2), X.getBase(2), X.getBound(2), NO2, k_info); #ifdef DEBUG_MODE_ON std::cout << "??? bounds JJJ on processor " << my_pe << ":\n" << " " << X.getLocalBase(2) << " " << X.getLocalBound(2) << " " << X.getBase(2) << " " << X.getBound(2) << " " << NO2 << "\n"; std::cout << "???? I,J,K info on processor " << my_pe << ":\n" << "i_info" << "\t" << "j_info" << "\t" << "k_info" << i_info[0] << "\t" << j_info[0] << k_info[0] << i_info[1] << "\t" << j_info[1] << k_info[1] << i_info[2] << "\t" << j_info[2] << k_info[2] << i_info[3] << "\t" << j_info[3] << k_info[3] << "\n"; #endif // Communication_Manager::Sync(); mhd_proc = FALSE; if (my_pe >= proc_first && my_pe <= proc_last ) { mhd_proc = TRUE; MHD_PROCESSOR = TRUE; } MHD_FORTRAN_ARRAYS(variable,i_info, j_info, k_info, &my_pe, &mhd_proc ); // zero out the data space where needed ZERO(X.getDataPointer()); ZERO(Y.getDataPointer()); ZERO(Z.getDataPointer()); ZERO(Rho.getDataPointer()); ZERO(Vx.getDataPointer()); ZERO(Vy.getDataPointer()); ZERO(Vz.getDataPointer()); ZERO(Csound.getDataPointer()); ZERO(Pressure.getDataPointer()); ZERO(Bx.getDataPointer()); ZERO(By.getDataPointer()); ZERO(Bz.getDataPointer()); ZERO(Bi.getDataPointer()); ZERO(Bj.getDataPointer()); ZERO(Bk.getDataPointer()); ZERO(Ei.getDataPointer()); ZERO(Ej.getDataPointer()); ZERO(Ek.getDataPointer()); #ifdef ARRAY_TEST TEST_ARRAYS(); #endif // Communication_Manager::Sync(); } /* MHD::initial(int nmhd, int proc_first, int proc_last) */ /***************************************************************************/ /// get_local_fortran information /** * The fortran call passes off the pointers to the data in the MHD * class through the "variable" pointer array. The local data array * range is passed off in the info arrays which contain [0] -> * local lower, [1] local high dimension, [2] -> local active high * dimension, [3] offset from local index to global index * * @param il - lowest cell recognized by P++ (local) * @param iup - highest cell recognized by P++ (local) * @param il_g - lower bound of global array * @param iup_g - upper bound of global array * * @param info[0] = local lowest guard cell * @param info[1] = local highest guard cell * @param info[2] = number of active cells * @param info[3] = offset in global array for cell 1 * @param info[4] = lowest active cell -globally ( = 1 ) * @param info[5] = upper bound on active cells - globally * * @see initial(int,int,int); * * @fixme this should be a private member function of class MHD. * * @fixme the "get_local_fortran_info" function should probably be a subclass of MHD * rather than a method of MHD. For example, we could have an MSETUP * sub-class. This may reduce some of the confusing preprocessor * directives in MHD.C. * * @note get_local_fortran_info is only defined with the preprocessor directive MSETUP. */ void get_local_fortran_info(int il, int iup, int il_g, int iup_g, int no2, int *info ) { #ifdef DEBUG_MODE_ON printf("DEBUG: process (%d) in MHD.C::MHD::get_local_fortran_info(...)\n", RTIME.local_proc); #endif // il - is lowest cell recognized by P++ (local) // iup - is highest cell recognized by P++ (local) // il_g - lower bound of global array // iup_g - upper bound of global array // info[0] = local lowest guard cell // info[1] = local highest guard cell // info[2] = number of active cells // info[3] = offset in global array for cell 1 // info[4] = lowest active cell -globally ( = 1 ) // info[5] = upper bound on active cells - globally info[4] = il_g; info[5] = iup_g -1; // P++ arrays dimensioned ni+1 // default -- center partition -- il is the lowest guard cell, iup -- is the highest guard cell info[0] = -no2 + 1; info[1]= iup - il + 1 - no2; info[2]= iup-il +1 - 2*no2; info[3] = il + no2-1; #ifndef PARA_PPP // no partitioning -- no guard cells evident for P++ info[0] = -no2+1; info[1] = iup-il+1+no2; info[2] = iup-il; info[3] = 0; #else // not partitioned in this direction if ( il == il_g && iup == iup_g ) { info[0] = -no2+1; info[1] = iup-il_g+1+no2; info[2] = iup-il; info[3] = 0; } else { // check for low side if ( il == il_g ) { info[0] = -no2+1; info[1] = iup-il+1; info[2] = iup-il+1-no2; info[3]=il-1; } // check for high side // our definitions have, e.g., ni+1 instead of ni as the size of arrays, but want to // classify the active cells on the basis of ni if (iup == iup_g ) { info[0] = -no2+1; info[1] = iup-il+1, info[2] = iup-il-no2; info[3] = il+ no2 -1; } } /* else */ #endif /* ifndef PARA_PPP ... else ... */ } /* get_local_fortran_info(int il, int iup, int il_g, int iup_g, int no2, int *info ) */ /***************************************************************************/ /// startup /** * Sets partitioning of BBX, BBY, BBZ on the processors. Updates ghost boundaries. * * @param filename * @note startup is only defined with the preprocessor directive MSETUP. * * @fixme the "startup" function should probably be a subclass of MHD * rather than a method of MHD. For example, we could have an MSETUP * sub-class. This may reduce some of the confusing preprocessor * directives in MHD.C. * */ void MHD::startup( char *filename) { #ifdef DEBUG_MODE_ON printf("DEBUG: process (%d) in MHD.C::MHD::startup(...)\n", RTIME.local_proc); #endif int i,j,k; // Dummy initialization of variables // would call hdftake here #ifdef DEBUG int mype = Communication_Manager::My_Process_Number; printf("in startup %d\n", mype); fflush(NULL); #endif TEST_ARRAYS(); int uno = 1; int NI2 = NI/2; Range I(1,NI+1); Range J(1,NJ+1); Range K(1,NK+1); Range Ilow(1,NI2); Range Ihigh(NI2+1,NI+1); // fill in boundaries for X,Y,Z FloatARRAY BBX(NI+1,NJ+1,NK+1); BBX.setBase(1); BBX.partition(MHD_P); FloatARRAY BBY(NI+1,NJ+1,NK+1); BBY.setBase(1); BBY.partition(MHD_P); FloatARRAY BBZ(NI+1,NJ+1,NK+1); BBZ.setBase(1); BBZ.partition(MHD_P); TEST_ARRAYS(); #ifdef DEBUG printf("before grid boundaries %d\n", mype); fflush(NULL); #endif Optimization_Manager::setOptimizedScalarIndexing(On); for ( k=1; k < NK+2; k++ ) { for (j=1; j < NJ+2; j++ ) { for (i=1; i < NO2+1; i++ ) { BBX(i,j,k) = X(1,j,k) - (X(i+1,j,k)-X(1,j,k)); BBY(i,j,k) = Y(1,j,k) - (Y(i+1,j,k)-Y(1,j,k)); BBZ(i,j,k) = Z(1,j,k) - (Z(i+1,j,k)-Z(1,j,k)); BBX(NI-NO2+i,j,k) = X(NI+1,j,k) - (X(NI+1-i,j,k)-X(NI+1,j,k)); BBY(NI-NO2+i,j,k) = Y(NI+1,j,k) - (Y(NI+1-i,j,k)-Y(NI+1,j,k)); BBZ(NI-NO2+i,j,k) = Z(NI+1,j,k) - (Z(NI+1-i,j,k)-Z(NI+1,j,k)); } } #ifdef DEBUG printf("in loop, after k %d %d\n",k, mype); fflush(NULL); #endif } #ifdef DEBUG printf("before bound calls %d\n", mype); fflush(NULL); #endif IBOUND(BBX.getDataPointer(),X.getDataPointer(),&uno,&uno,&uno); IBOUND(BBY.getDataPointer(),Y.getDataPointer(),&uno,&uno,&uno); IBOUND(BBZ.getDataPointer(),Z.getDataPointer(),&uno,&uno,&uno); #ifdef DEBUG printf("after loop %d\n", mype); fflush(NULL); #endif Range Kl(1,NK/2); Range Kh(NK/2+1,NK); Range J_1(1,NO2); Range J_2(NJ-NO2+1,NJ); Range J_P(1,NJ+1); TEST_ARRAYS(); #ifdef CARTESIAN for ( k=1; k < NK+2; k++ ) { for (j=1; j < NO2+1; j++ ) { for (i=1; i < NI+2; i++ ) { BBX(i,j,k) = X(i,1,k) - (X(i,j+1,k)-X(i,1,k)); BBY(i,j,k) = Y(i,1,k) - (Y(i,j+1,k)-Y(i,1,k)); BBZ(i,j,k) = Z(i,1,k) - (Z(i,j+1,k)-Z(i,1,k)); BBX(i,NJ-NO2+j,k) = X(i,NJ+1,k) - (X(i,NJ+1-j,k)-X(i,NJ+1,k)); BBY(i,NJ-NO2+j,k) = Y(i,NJ+1,k) - (Y(i,NJ+1-j,k)-Y(i,NJ+1,k)); BBZ(i,NJ-NO2+j,k) = Z(i,NJ+1,k) - (Z(i,NJ+1-j,k)-Z(i,NJ+1,k)); } } } #else /* #ifdef CARTESIAN */ BBX(I,J_1,Kl) = X(I,J_1+1,Kh); BBX(I,J_1,Kh) = X(I,J_1+1,Kl); BBX(I,J_2,Kl) = X(I,J_2,Kh); BBX(I,J_2,Kh) = X(I,J_2,Kl); BBX(I,NJ+1,K) = X(I,NJ+1,K); BBX(I,J_P,NK+1) = BBX(I,J_P,1); BBY(I,J_1,Kl) = Y(I,J_1+1,Kh); BBY(I,J_1,Kh) = Y(I,J_1+1,Kl); BBY(I,J_2,Kl) = Y(I,J_2,Kh); BBY(I,J_2,Kh) = Y(I,J_2,Kl); BBY(I,NJ+1,K) = Y(I,NJ+1,K); BBY(I,J_P,NK+1) = BBY(I,J_P,1); BBZ(I,J_1,Kl) = Z(I,J_1+1,Kh); BBZ(I,J_1,Kh) = Z(I,J_1+1,Kl); BBZ(I,J_2,Kl) = Z(I,J_2,Kh); BBZ(I,J_2,Kh) = Z(I,J_2,Kl); BBZ(I,NJ+1,K) = Z(I,NJ+1,K); BBZ(I,J_P,NK+1) = BBZ(I,J_P,1); #endif /* #ifdef CARTESIAN */ JBOUND(BBX.getDataPointer(),X.getDataPointer(),&uno,&uno,&uno); JBOUND(BBY.getDataPointer(),Y.getDataPointer(),&uno,&uno,&uno); JBOUND(BBZ.getDataPointer(),Z.getDataPointer(),&uno,&uno,&uno); Range JJ(1,NJ+1); Range K_0(1,NO2); Range K_1(NK-NO2+1,NK); TEST_ARRAYS(); #ifdef CARTESIAN for ( k=1; k < NO2+1; k++ ) { for (j=1; j < NJ+2; j++ ) { for (i=1; i < NI+1; i++ ) { BBX(i,j,k) = X(i,j,1) - (X(i,j,k+1)-X(i,j,1)); BBY(i,j,k) = Y(i,j,1) - (Y(i,j,k+1)-Y(i,j,1)); BBZ(i,j,k) = Z(i,j,1) - (Z(i,j,k+1)-Z(i,j,1)); BBX(i,j,NK-NO2+k) = X(i,j,NK+1) - (X(i,j,NK+1-k)-X(i,j,NK+1)); BBY(i,j,NK-NO2+k) = Y(i,j,NK+1) - (Y(i,j,NK+1-k)-Y(i,j,NK+1)); BBZ(i,j,NK-NO2+k) = Z(i,j,NK+1) - (Z(i,j,NK+1-k)-Z(i,j,NK+1)); } } } #else /* #ifdef CARTESIAN */ BBX(I,JJ,K_0) = X(I,JJ,K_1); BBX(I,JJ,K_1) = X(I,JJ,K_0+1); BBY(I,JJ,K_0) = Y(I,JJ,K_1); BBY(I,JJ,K_1) = Y(I,JJ,K_0+1); BBZ(I,JJ,K_0) = Z(I,JJ,K_1); BBZ(I,JJ,K_1) = Z(I,JJ,K_0+1); #endif /* ifdef CARTESIAN */ KBOUND(BBY.getDataPointer(), Y.getDataPointer(),&uno,&uno,&uno); KBOUND(BBZ.getDataPointer(), Z.getDataPointer(),&uno,&uno,&uno); KBOUND(BBX.getDataPointer(), X.getDataPointer(),&uno,&uno,&uno); Optimization_Manager::setOptimizedScalarIndexing(Off); #ifdef DEBUG printf("before update ghost %d\n", mype); fflush(NULL); #endif X.updateGhostBoundaries(); Y.updateGhostBoundaries(); Z.updateGhostBoundaries(); TEST_ARRAYS(); #ifdef DEBUG printf("before metric %d\n", mype); fflush(NULL); #endif // do all the good metric stuff METRIC(); #ifdef DEBUG printf("after metric %d\n", mype); fflush(NULL); #endif Volume.updateGhostBoundaries(); Volume_inv.updateGhostBoundaries(); Face_i.updateGhostBoundaries(); Face_j.updateGhostBoundaries(); Face_k.updateGhostBoundaries(); Tran_i.updateGhostBoundaries(); Tran_j.updateGhostBoundaries(); Tran_k.updateGhostBoundaries(); Edge_i.updateGhostBoundaries(); Edge_j.updateGhostBoundaries(); Edge_k.updateGhostBoundaries(); Vec_i.updateGhostBoundaries(); Vec_j.updateGhostBoundaries(); Vec_k.updateGhostBoundaries(); Cell_Length.updateGhostBoundaries(); INITAL(); // end of dummy initialization Communication_Manager::Sync(); this->boundary(); #ifdef DEBUG printf("leaving the this gusto %d\n", mype); fflush(stdout); #endif } /* void MHD::startup( char *filename) */ /***************************************************************************/ /// step /** * Make the magnetic fluxes accurate to double precision roundoff. * Call ringav() if DO_RINGAV preprocessor directive is set. * * @param Dt * @note step is only defined with the preprocessor directive MSETUP. * * @fixme the "startup" function should probably be a subclass of MHD * rather than a method of MHD. For example, we could have an MSETUP * sub-class. This may reduce some of the confusing preprocessor * directives in MHD.C. * * @fixme uses MHD_HACK() */ void MHD::step(FLoaT *Dt) { #ifdef DEBUG_MODE_ON printf("DEBUG: process (%d) in MHD.C::MHD::step(FLoaT *Dt)\n", RTIME.local_proc); #endif MHD_HACK(); #ifdef DO_RINGAV /* order is important here! ringav averages the cell centered quanitities, but only computes the electric field necessary to average the magnetic face fluxes. This is added to the E's already computed and, so, the actual averaging gets done in FLUX_FIX */ RINGAV(); #endif /* #ifdef DO_RINGAV */ // code to make the magnetic fluxes accurate to double precision roundoff Ei.updateGhostBoundaries(); Ej.updateGhostBoundaries(); Ek.updateGhostBoundaries(); FLUX_FIX(Dt); Bi.updateGhostBoundaries(); Bj.updateGhostBoundaries(); Bk.updateGhostBoundaries(); } /* void MHD::step(FLoaT *Dt) */ #endif /* #ifdef MSETUP */ // goes to below constructor -- stuff we don't need for setup calls /***************************************************************************/ /// dummy_start /** * Dummy initialization of variables * * Initializes X,Y,Z, Vx,Vy,Vz, rho, Bx,By,Bz and Csound * */ void MHD::dummy_start() { #ifdef DEBUG_MODE_ON printf("DEBUG: process (%d) in MHD.C::MHD::dummy_start(...)\n", RTIME.local_proc); #endif int i,j,k; FLoaT r,ct,st,cp,sp; // Dummy initialization of variables // would call hdftake here int NI2 = NI/2; Range I(1,NI+1); Range J(1,NJ+1); Range K(1,NK+1); Range Ilow(1,NI2); Range Ihigh(NI2+1,NI+1); Optimization_Manager::setOptimizedScalarIndexing(On); /* FIXME: What calculation is going on here? It looks like some * sort of transformation of variables... */ for ( i=1; i