From 5da02d40dca8d10328cd39d340fa65e5070e222b Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Wed, 10 Jan 2024 19:53:48 -0500 Subject: [PATCH] Replace tabs with spaces (#2708) --- CITATION.md | 78 +-- Diagnostics/DustCollapse/main.cpp | 470 +++++++++--------- Diagnostics/Radiation/Radiation_utils.H | 216 ++++---- Diagnostics/Radiation/gaussian_pulse.cpp | 192 +++---- Diagnostics/Radiation/lgt_frnt1d.cpp | 186 +++---- Diagnostics/Radiation/rad_shock.cpp | 280 +++++------ Diagnostics/Radiation/rad_source.cpp | 114 ++--- Diagnostics/Radiation/rad_sphere.cpp | 424 ++++++++-------- Diagnostics/Radiation/rhd_shocktube.cpp | 238 ++++----- .../riemann_2d/problem_initialize.H | 14 +- .../problem_initialize_state_data.H | 158 +++--- .../test_convect/problem_diagnostics.H | 32 +- .../toy_convect/problem_diagnostics.H | 32 +- .../Rad2Tshock/python/LowrieEdwardsUnits.py | 8 +- .../Rad2Tshock/python/RadShock.py | 6 +- .../Rad2Tshock/python/initincgs.py | 6 +- .../RadShestakovBolstad/python/initunits.py | 6 +- .../RadSuOlsonMG/python/initunits.py | 6 +- .../problem_initialize_state_data.H | 2 +- .../problem_initialize_state_data.H | 52 +- Exec/unit_tests/particles_test/Prob.cpp | 132 ++--- Source/driver/sum_utils.cpp | 8 +- Source/driver/timestep.cpp | 4 +- Util/ConvertCheckpoint/Embiggen.cpp | 34 +- Util/code_checker/tab_exterminator.sh | 6 +- Util/model_parser/model_parser.H | 22 +- 26 files changed, 1363 insertions(+), 1363 deletions(-) diff --git a/CITATION.md b/CITATION.md index b38fe9a5cd..160568d222 100644 --- a/CITATION.md +++ b/CITATION.md @@ -29,23 +29,23 @@ provides some details on the algorithmic implementations): ``` @ARTICLE{2010ApJ...715.1221A, author = {{Almgren}, A.~S. and {Beckner}, V.~E. and {Bell}, - J.~B. and {Day}, M.~S. and {Howell}, L.~H. and - {Joggerst}, C.~C. and {Lijewski}, M.~J. and - {Nonaka}, A. and {Singer}, M. and {Zingale}, M.}, - title = "{CASTRO: A New Compressible Astrophysical - Solver. I. Hydrodynamics and Self-gravity}", + J.~B. and {Day}, M.~S. and {Howell}, L.~H. and + {Joggerst}, C.~C. and {Lijewski}, M.~J. and + {Nonaka}, A. and {Singer}, M. and {Zingale}, M.}, + title = "{CASTRO: A New Compressible Astrophysical + Solver. I. Hydrodynamics and Self-gravity}", journal = {\apj}, archivePrefix = "arXiv", eprint = {1005.0114}, primaryClass = "astro-ph.IM", keywords = {equation of state, gravitation, hydrodynamics, methods: - numerical, nuclear reactions, nucleosynthesis, - abundances}, - year = 2010, - month = jun, + numerical, nuclear reactions, nucleosynthesis, + abundances}, + year = 2010, + month = jun, volume = 715, - pages = {1221-1238}, - doi = {10.1088/0004-637X/715/2/1221}, + pages = {1221-1238}, + doi = {10.1088/0004-637X/715/2/1221}, adsurl = {http://adsabs.harvard.edu/abs/2010ApJ...715.1221A}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -64,21 +64,21 @@ cite the following: ``` @ARTICLE{2011ApJS..196...20Z, author = {{Zhang}, W. and {Howell}, L. and {Almgren}, A. and - {Burrows}, A. and {Bell}, J.}, - title = "{CASTRO: A New Compressible Astrophysical - Solver. II. Gray Radiation Hydrodynamics}", + {Burrows}, A. and {Bell}, J.}, + title = "{CASTRO: A New Compressible Astrophysical + Solver. II. Gray Radiation Hydrodynamics}", journal = {\apjs}, archivePrefix = "arXiv", eprint = {1105.2466}, primaryClass = "astro-ph.IM", keywords = {diffusion, hydrodynamics, methods: numerical, radiative - transfer}, - year = 2011, - month = oct, + transfer}, + year = 2011, + month = oct, volume = 196, - eid = {20}, - pages = {20}, - doi = {10.1088/0067-0049/196/2/20}, + eid = {20}, + pages = {20}, + doi = {10.1088/0067-0049/196/2/20}, adsurl = {http://adsabs.harvard.edu/abs/2011ApJS..196...20Z}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -87,21 +87,21 @@ cite the following: ``` @ARTICLE{2013ApJS..204....7Z, author = {{Zhang}, W. and {Howell}, L. and {Almgren}, A. and - {Burrows}, A. and {Dolence}, J. and {Bell}, J.}, - title = "{CASTRO: A New Compressible Astrophysical - Solver. III. Multigroup Radiation Hydrodynamics}", + {Burrows}, A. and {Dolence}, J. and {Bell}, J.}, + title = "{CASTRO: A New Compressible Astrophysical + Solver. III. Multigroup Radiation Hydrodynamics}", journal = {\apjs}, archivePrefix = "arXiv", eprint = {1207.3845}, primaryClass = "astro-ph.IM", keywords = {diffusion, hydrodynamics, methods: numerical, radiative - transfer }, - year = 2013, - month = jan, + transfer }, + year = 2013, + month = jan, volume = 204, - eid = {7}, - pages = {7}, - doi = {10.1088/0067-0049/204/1/7}, + eid = {7}, + pages = {7}, + doi = {10.1088/0067-0049/204/1/7}, adsurl = {http://adsabs.harvard.edu/abs/2013ApJS..204....7Z}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } @@ -116,22 +116,22 @@ hydro and gravity: ``` @ARTICLE{2016ApJ...819...94K, author = {{Katz}, M.~P. and {Zingale}, M. and {Calder}, A.~C. and - {Swesty}, F.~D. and {Almgren}, A.~S. and {Zhang}, - W.}, - title = "{White Dwarf Mergers on Adaptive Meshes. I. Methodology - and Code Verification}", + {Swesty}, F.~D. and {Almgren}, A.~S. and {Zhang}, + W.}, + title = "{White Dwarf Mergers on Adaptive Meshes. I. Methodology + and Code Verification}", journal = {\apj}, archivePrefix = "arXiv", eprint = {1512.06099}, primaryClass = "astro-ph.HE", keywords = {hydrodynamics, methods: numerical, supernovae: general, - white dwarfs}, - year = 2016, - month = mar, + white dwarfs}, + year = 2016, + month = mar, volume = 819, - eid = {94}, - pages = {94}, - doi = {10.3847/0004-637X/819/2/94}, + eid = {94}, + pages = {94}, + doi = {10.3847/0004-637X/819/2/94}, adsurl = {http://adsabs.harvard.edu/abs/2016ApJ...819...94K}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } diff --git a/Diagnostics/DustCollapse/main.cpp b/Diagnostics/DustCollapse/main.cpp index 49ed554e26..7e2a0763de 100644 --- a/Diagnostics/DustCollapse/main.cpp +++ b/Diagnostics/DustCollapse/main.cpp @@ -3,13 +3,13 @@ // and output the position of the interface as a function of time. // // For 2d: -// The initial dense sphere is assumed to be centered a r = 0 (x = 0). -// We take as default that it is centered vertically at y = 0, but -// this can be overridden with --yctr. +// The initial dense sphere is assumed to be centered a r = 0 (x = 0). +// We take as default that it is centered vertically at y = 0, but +// this can be overridden with --yctr. // // For 3d: -// The initial dense sphere is assumed to be centered a x = y = z = 0, -// but this can be overridden with --{x,y,z}ctr. +// The initial dense sphere is assumed to be centered a x = y = z = 0, +// but this can be overridden with --{x,y,z}ctr. // #include #include @@ -30,290 +30,290 @@ Vector GetCenter (const string pltfile); int main(int argc, char* argv[]) { - amrex::Initialize(argc, argv, false); + amrex::Initialize(argc, argv, false); - // timer for profiling - BL_PROFILE_VAR("main()", pmain); + // timer for profiling + BL_PROFILE_VAR("main()", pmain); - // Input arguments - if (argc < 2) - Abort("ERROR: Missing plotfiles"); + // Input arguments + if (argc < 2) + Abort("ERROR: Missing plotfiles"); - Print() << "\nUsing a density threshold of half the maximum analytic density" << std::endl; + Print() << "\nUsing a density threshold of half the maximum analytic density" << std::endl; - auto farg = 1; - bool profile = false; + auto farg = 1; + bool profile = false; #if (AMREX_SPACEDIM >= 2) - auto j = 1; - while (j < argc) { - if ( !strcmp(argv[j], "--profile") ) - { - profile = true;; - farg++; - } - // Go to the next parameter name - ++j; - } + auto j = 1; + while (j < argc) { + if ( !strcmp(argv[j], "--profile") ) + { + profile = true;; + farg++; + } + // Go to the next parameter name + ++j; + } #endif - // loop over plotfiles - for (auto f = farg; f < argc; f++) { + // loop over plotfiles + for (auto f = farg; f < argc; f++) { - string pltfile = argv[f]; + string pltfile = argv[f]; auto center = GetCenter(pltfile); - Real xctr = center[0]; - Real yctr = center[1]; - Real zctr = center[2]; + Real xctr = center[0]; + Real yctr = center[1]; + Real zctr = center[2]; - // Start dataservices - DataServices::SetBatchMode(); + // Start dataservices + DataServices::SetBatchMode(); - // Define the type of file - Amrvis::FileType fileType(Amrvis::NEWPLT); - DataServices dataServices (pltfile, fileType); + // Define the type of file + Amrvis::FileType fileType(Amrvis::NEWPLT); + DataServices dataServices (pltfile, fileType); - if (!dataServices.AmrDataOk()) - DataServices::Dispatch(DataServices::ExitRequest, NULL); + if (!dataServices.AmrDataOk()) + DataServices::Dispatch(DataServices::ExitRequest, NULL); - // get data from plot file - AmrData& data = dataServices.AmrDataRef(); + // get data from plot file + AmrData& data = dataServices.AmrDataRef(); - int finestLevel = data.FinestLevel(); + int finestLevel = data.FinestLevel(); - // get variable names - const Vector& varNames = data.PlotVarNames(); + // get variable names + const Vector& varNames = data.PlotVarNames(); - // get the index bounds and dx. - Box domain = data.ProbDomain()[finestLevel]; - auto dx = data.CellSize(finestLevel); - const Vector& problo = data.ProbLo(); - const Vector& probhi = data.ProbHi(); + // get the index bounds and dx. + Box domain = data.ProbDomain()[finestLevel]; + auto dx = data.CellSize(finestLevel); + const Vector& problo = data.ProbLo(); + const Vector& probhi = data.ProbHi(); - Vector rr = data.RefRatio(); + Vector rr = data.RefRatio(); - auto rmin = problo[0]; + auto rmin = problo[0]; - const auto dx_fine = *(std::min_element(dx.begin(), dx.end())); + const auto dx_fine = *(std::min_element(dx.begin(), dx.end())); - // compute the maximum number of zones, as if we were completely refined + // compute the maximum number of zones, as if we were completely refined #if (AMREX_SPACEDIM == 1) - int nbins = domain.length(0); + int nbins = domain.length(0); #elif (AMREX_SPACEDIM == 2) - auto x_maxdist = fmax(fabs(problo[0]), fabs(probhi[0])); - auto y_maxdist = fmax(fabs(problo[1] - yctr), fabs(probhi[1] - yctr)); + auto x_maxdist = fmax(fabs(problo[0]), fabs(probhi[0])); + auto y_maxdist = fmax(fabs(problo[1] - yctr), fabs(probhi[1] - yctr)); - auto max_dist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist); + auto max_dist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist); - int nbins = int(max_dist / dx_fine); + int nbins = int(max_dist / dx_fine); #else - auto x_maxdist = fmax(fabs(problo[0] - xctr), fabs(probhi[0] - xctr)); - auto y_maxdist = fmax(fabs(problo[1] - yctr), fabs(probhi[1] - yctr)); - auto z_maxdist = fmax(fabs(problo[2] - zctr), fabs(probhi[2] - zctr)); + auto x_maxdist = fmax(fabs(problo[0] - xctr), fabs(probhi[0] - xctr)); + auto y_maxdist = fmax(fabs(problo[1] - yctr), fabs(probhi[1] - yctr)); + auto z_maxdist = fmax(fabs(problo[2] - zctr), fabs(probhi[2] - zctr)); - auto max_dist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist + - z_maxdist*z_maxdist); + auto max_dist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist + + z_maxdist*z_maxdist); - int nbins = int(max_dist / dx_fine); + int nbins = int(max_dist / dx_fine); #endif - // radial coordinate - Vector r(nbins); - for (auto i = 0; i < nbins; i++) - r[i] = (i + 0.5) * dx_fine + rmin; + // radial coordinate + Vector r(nbins); + for (auto i = 0; i < nbins; i++) + r[i] = (i + 0.5) * dx_fine + rmin; - // find variable indices - auto dens_comp = data.StateNumber("density"); + // find variable indices + auto dens_comp = data.StateNumber("density"); - if (dens_comp < 0 ) - Abort("ERROR: density variable not found"); + if (dens_comp < 0 ) + Abort("ERROR: density variable not found"); - // allocate storage for data - Vector dens(nbins, 0.); - Vector volcount(nbins, 0); + // allocate storage for data + Vector dens(nbins, 0.); + Vector volcount(nbins, 0); - // r1 is the factor between the current level grid spacing and the - // FINEST level - auto r1 = 1.0; + // r1 is the factor between the current level grid spacing and the + // FINEST level + auto r1 = 1.0; - Vector fill_comps(data.NComp()); - for (auto i = 0; i < data.NComp(); i++) - fill_comps[i] = i; + Vector fill_comps(data.NComp()); + for (auto i = 0; i < data.NComp(); i++) + fill_comps[i] = i; - // imask will be set to false if we've already output the data. - // Note, imask is defined in terms of the finest level. As we loop - // over levels, we will compare to the finest level index space to - // determine if we've already output here - int mask_size = domain.length().max(); - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + // imask will be set to false if we've already output the data. + // Note, imask is defined in terms of the finest level. As we loop + // over levels, we will compare to the finest level index space to + // determine if we've already output here + int mask_size = domain.length().max(); + Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); - // counter - int cnt = 0; + // counter + int cnt = 0; - // loop over the data, starting at the finest grid, and if we haven't - // already stored data in that grid location (according to imask), - // store it. - for (int l = finestLevel; l >= 0; l--) { + // loop over the data, starting at the finest grid, and if we haven't + // already stored data in that grid location (according to imask), + // store it. + for (int l = finestLevel; l >= 0; l--) { - Vector level_dx = data.DxLevel()[l]; + Vector level_dx = data.DxLevel()[l]; - const BoxArray& ba = data.boxArray(l); - const DistributionMapping& dm = data.DistributionMap(l); + const BoxArray& ba = data.boxArray(l); + const DistributionMapping& dm = data.DistributionMap(l); - MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); - data.FillVar(lev_data_mf, l, varNames, fill_comps); + MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); + data.FillVar(lev_data_mf, l, varNames, fill_comps); #ifdef _OPENMP #pragma omp parallel #endif - for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); + for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); #if (AMREX_SPACEDIM == 1) - fdustcollapse1d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens.dataPtr(), - imask.dataPtr(), mask_size, r1, dens_comp, &cnt); + fdustcollapse1d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), + nbins, dens.dataPtr(), + imask.dataPtr(), mask_size, r1, dens_comp, &cnt); #elif (AMREX_SPACEDIM == 2) - fdustcollapse2d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens.dataPtr(), volcount.dataPtr(), - imask.dataPtr(), mask_size, r1, - ZFILL(level_dx), dx_fine, yctr, dens_comp); + fdustcollapse2d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), + nbins, dens.dataPtr(), volcount.dataPtr(), + imask.dataPtr(), mask_size, r1, + ZFILL(level_dx), dx_fine, yctr, dens_comp); #else - fdustcollapse3d(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens.dataPtr(), volcount.dataPtr(), - imask.dataPtr(), mask_size, r1, - level_dx.dataPtr(), dx_fine, xctr, yctr, zctr, dens_comp); + fdustcollapse3d(bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), + nbins, dens.dataPtr(), volcount.dataPtr(), + imask.dataPtr(), mask_size, r1, + level_dx.dataPtr(), dx_fine, xctr, yctr, zctr, dens_comp); #endif - } + } - // adjust r1 for the next lowest level - if (l != 0) r1 *= rr[l-1]; - } + // adjust r1 for the next lowest level + if (l != 0) r1 *= rr[l-1]; + } #if (AMREX_SPACEDIM == 1) - // sort the data based on the coordinates - auto isv = sort_indexes(r); + // sort the data based on the coordinates + auto isv = sort_indexes(r); #else - // normalize - for (int i = 0; i < nbins; i++) - if (volcount[i] != 0.) - dens[i] /= volcount[i]; + // normalize + for (int i = 0; i < nbins; i++) + if (volcount[i] != 0.) + dens[i] /= volcount[i]; #endif - // These are calculated analytically given initial density 1.e9 and the - // analytic expression for the radius as a function of time t = 0.00 - Real max_dens = 1.e9; - - if (fabs(data.Time()) <= 1.e-8) - max_dens = 1.e9; - else if (fabs(data.Time() - 0.01) <= 1.e-8) - max_dens = 1.043345e9; - else if (fabs(data.Time() - 0.02) <= 1.e-8) - max_dens = 1.192524e9; - else if (fabs(data.Time() - 0.03) <= 1.e-8) - max_dens = 1.527201e9; - else if (fabs(data.Time() - 0.04) <= 1.e-8) - max_dens = 2.312884e9; - else if (fabs(data.Time() - 0.05) <= 1.e-8) - max_dens = 4.779133e9; - else if (fabs(data.Time() - 0.06) <= 1.e-8) - max_dens = 24.472425e9; - else if (fabs(data.Time() - 0.065) <= 1.e-8) - max_dens = 423.447291e9; - else { - Print() << "Dont know the maximum density at this time: " << data.Time() < Vector sort_indexes(const Vector &v) { - // initialize original index locations - Vector idx(v.size()); - iota(idx.begin(), idx.end(), 0); + // initialize original index locations + Vector idx(v.size()); + iota(idx.begin(), idx.end(), 0); - // sort indexes based on comparing values in v - sort(idx.begin(), idx.end(), - [&v](size_t i1, size_t i2) { - return v[i1] < v[i2]; - }); + // sort indexes based on comparing values in v + sort(idx.begin(), idx.end(), + [&v](size_t i1, size_t i2) { + return v[i1] < v[i2]; + }); - return idx; + return idx; } /// @@ -321,40 +321,40 @@ Vector sort_indexes(const Vector &v) { /// string /// string GetVarFromJobInfo (const string pltfile, const string varname) { - string filename = pltfile + "/job_info"; - std::regex re("(?:[ \\t]*)" + varname + "\\s*:\\s*(.*)\\s*\\n"); - - std::smatch m; - - std::ifstream jobfile(filename); - if (jobfile.is_open()) { - std::stringstream buf; - buf << jobfile.rdbuf(); - string file_contents = buf.str(); - - if (std::regex_search(file_contents, m, re)) { - return m[1]; - } else { - Print() << "Unable to find " << varname << " in job_info file!" << std::endl; - } - } else { - Print() << "Could not open job_info file!" << std::endl; - } - - return ""; + string filename = pltfile + "/job_info"; + std::regex re("(?:[ \\t]*)" + varname + "\\s*:\\s*(.*)\\s*\\n"); + + std::smatch m; + + std::ifstream jobfile(filename); + if (jobfile.is_open()) { + std::stringstream buf; + buf << jobfile.rdbuf(); + string file_contents = buf.str(); + + if (std::regex_search(file_contents, m, re)) { + return m[1]; + } else { + Print() << "Unable to find " << varname << " in job_info file!" << std::endl; + } + } else { + Print() << "Could not open job_info file!" << std::endl; + } + + return ""; } // Get the center from the job info file and return as a Real Vector Vector GetCenter (const string pltfile) { - auto center_str = GetVarFromJobInfo(pltfile, "center"); + auto center_str = GetVarFromJobInfo(pltfile, "center"); - // split string - std::istringstream iss {center_str}; - Vector center; + // split string + std::istringstream iss {center_str}; + Vector center; - std::string s; - while (std::getline(iss, s, ',')) - center.push_back(stod(s)); + std::string s; + while (std::getline(iss, s, ',')) + center.push_back(stod(s)); - return center; + return center; } diff --git a/Diagnostics/Radiation/Radiation_utils.H b/Diagnostics/Radiation/Radiation_utils.H index 42b60210ae..af017f86f7 100644 --- a/Diagnostics/Radiation/Radiation_utils.H +++ b/Diagnostics/Radiation/Radiation_utils.H @@ -24,17 +24,17 @@ Vector GetCenter (const string pltfile); template Vector sort_indexes(const Vector &v) { - // initialize original index locations - Vector idx(v.size()); - iota(idx.begin(), idx.end(), 0); + // initialize original index locations + Vector idx(v.size()); + iota(idx.begin(), idx.end(), 0); - // sort indexes based on comparing values in v - sort(idx.begin(), idx.end(), - [&v](size_t i1, size_t i2) { - return v[i1] < v[i2]; - }); + // sort indexes based on comparing values in v + sort(idx.begin(), idx.end(), + [&v](size_t i1, size_t i2) { + return v[i1] < v[i2]; + }); - return idx; + return idx; } // @@ -44,43 +44,43 @@ void GetInputArgs ( const int argc, char** argv, string& pltfile, string& slcfile, int& dir) { - int i = 1; // skip program name - - while ( i < argc) - { - - if ( !strcmp(argv[i], "-p") || !strcmp(argv[i],"--pltfile") ) - { - pltfile = argv[++i]; - } - else if ( !strcmp(argv[i], "-s") || !strcmp(argv[i],"--slicefile") ) - { - slcfile = argv[++i]; - } - else if ( !strcmp(argv[i], "-d") || !strcmp(argv[i],"--direction") ) - { - dir = std::atoi(argv[++i]); - } - else - { - std::cout << "\n\nOption " << argv[i] << " not recognized" << std::endl; - PrintHelp (); - exit ( EXIT_FAILURE ); - } - - // Go to the next parameter name - ++i; - } - - if (pltfile.empty() && slcfile.empty()) - { - PrintHelp(); - Abort("Missing input file"); - } - - Print() << "\nplotfile = \"" << pltfile << "\"" << std::endl; - Print() << "slicefile = \"" << slcfile << "\"" << std::endl; - Print() << std::endl; + int i = 1; // skip program name + + while ( i < argc) + { + + if ( !strcmp(argv[i], "-p") || !strcmp(argv[i],"--pltfile") ) + { + pltfile = argv[++i]; + } + else if ( !strcmp(argv[i], "-s") || !strcmp(argv[i],"--slicefile") ) + { + slcfile = argv[++i]; + } + else if ( !strcmp(argv[i], "-d") || !strcmp(argv[i],"--direction") ) + { + dir = std::atoi(argv[++i]); + } + else + { + std::cout << "\n\nOption " << argv[i] << " not recognized" << std::endl; + PrintHelp (); + exit ( EXIT_FAILURE ); + } + + // Go to the next parameter name + ++i; + } + + if (pltfile.empty() && slcfile.empty()) + { + PrintHelp(); + Abort("Missing input file"); + } + + Print() << "\nplotfile = \"" << pltfile << "\"" << std::endl; + Print() << "slicefile = \"" << slcfile << "\"" << std::endl; + Print() << std::endl; } // @@ -88,12 +88,12 @@ void GetInputArgs ( const int argc, char** argv, // Vector GetComponents(AmrData& data, const Vector varNames) { - Vector varComps; + Vector varComps; - for (auto it=varNames.begin(); it!=varNames.end(); ++it) - varComps.push_back(data.StateNumber(*it)); + for (auto it=varNames.begin(); it!=varNames.end(); ++it) + varComps.push_back(data.StateNumber(*it)); - return varComps; + return varComps; } @@ -105,36 +105,36 @@ void WriteSlicefile(const int nbins, const Vector r, Vector > vars, const std::string slcfile) { - // now open the slicefile and write out the data - std::ofstream slicefile; - slicefile.open(slcfile); - slicefile.setf(std::ios::scientific); - slicefile.precision(12); - const auto w = 24; + // now open the slicefile and write out the data + std::ofstream slicefile; + slicefile.open(slcfile); + slicefile.setf(std::ios::scientific); + slicefile.precision(12); + const auto w = 24; - // write the header - slicefile << std::setw(w) << "x"; - for (auto it=varNames.begin(); it!=varNames.end(); ++it) - slicefile << std::setw(w) << *it; + // write the header + slicefile << std::setw(w) << "x"; + for (auto it=varNames.begin(); it!=varNames.end(); ++it) + slicefile << std::setw(w) << *it; - slicefile << std::endl; + slicefile << std::endl; - // write the data in columns - const auto SMALL = 1.e-20; - for (auto i = 0; i < nbins; i++) { + // write the data in columns + const auto SMALL = 1.e-20; + for (auto i = 0; i < nbins; i++) { - slicefile << std::setw(w) << r[i]; + slicefile << std::setw(w) << r[i]; - for (auto it=vars.begin(); it!=vars.end(); ++it) { - if(fabs((*it)[i]) < SMALL) (*it)[i] = 0.0; - slicefile << std::setw(w) << (*it)[i]; - } + for (auto it=vars.begin(); it!=vars.end(); ++it) { + if(fabs((*it)[i]) < SMALL) (*it)[i] = 0.0; + slicefile << std::setw(w) << (*it)[i]; + } - slicefile << std::endl; + slicefile << std::endl; - } + } - slicefile.close(); + slicefile.close(); } @@ -143,13 +143,13 @@ void WriteSlicefile(const int nbins, const Vector r, // void PrintHelp () { - Print() << "\nusage: executable_name args" - << "\nargs [-p|--pltfile] plotfile : plot file directory (required)" - << "\n [-s|--slicefile] slice file : slice file (required)" - << "\n [--xctr] xctr : central x coord (non-cartesian only)" - << "\n [--yctr] yctr : central y coord (non-cartesian only)" - << "\n [--zctr] zctr : central z coord (non-cartesian only)" - << "\n\n" << std::endl; + Print() << "\nusage: executable_name args" + << "\nargs [-p|--pltfile] plotfile : plot file directory (required)" + << "\n [-s|--slicefile] slice file : slice file (required)" + << "\n [--xctr] xctr : central x coord (non-cartesian only)" + << "\n [--yctr] yctr : central y coord (non-cartesian only)" + << "\n [--zctr] zctr : central z coord (non-cartesian only)" + << "\n\n" << std::endl; } @@ -158,42 +158,42 @@ void PrintHelp () /// string /// string GetVarFromJobInfo (const string pltfile, const string varname) { - string filename = pltfile + "/job_info"; - std::regex re("(?:[ \\t]*)" + varname + "\\s*:\\s*(.*)\\s*\\n"); - - std::smatch m; - - std::ifstream jobfile(filename); - if (jobfile.is_open()) { - std::stringstream buf; - buf << jobfile.rdbuf(); - string file_contents = buf.str(); - - if (std::regex_search(file_contents, m, re)) { - return m[1]; - } else { - Print() << "Unable to find " << varname << " in job_info file!" << std::endl; - } - } else { - Print() << "Could not open job_info file!" << std::endl; - } - - return ""; + string filename = pltfile + "/job_info"; + std::regex re("(?:[ \\t]*)" + varname + "\\s*:\\s*(.*)\\s*\\n"); + + std::smatch m; + + std::ifstream jobfile(filename); + if (jobfile.is_open()) { + std::stringstream buf; + buf << jobfile.rdbuf(); + string file_contents = buf.str(); + + if (std::regex_search(file_contents, m, re)) { + return m[1]; + } else { + Print() << "Unable to find " << varname << " in job_info file!" << std::endl; + } + } else { + Print() << "Could not open job_info file!" << std::endl; + } + + return ""; } // Get the center from the job info file and return as a Real Vector Vector GetCenter (const string pltfile) { - auto center_str = GetVarFromJobInfo(pltfile, "center"); + auto center_str = GetVarFromJobInfo(pltfile, "center"); - // split string - std::istringstream iss {center_str}; - Vector center; + // split string + std::istringstream iss {center_str}; + Vector center; - std::string s; - while (std::getline(iss, s, ',')) - center.push_back(stod(s)); + std::string s; + while (std::getline(iss, s, ',')) + center.push_back(stod(s)); - return center; + return center; } #endif diff --git a/Diagnostics/Radiation/gaussian_pulse.cpp b/Diagnostics/Radiation/gaussian_pulse.cpp index 22e537d296..b04020a567 100644 --- a/Diagnostics/Radiation/gaussian_pulse.cpp +++ b/Diagnostics/Radiation/gaussian_pulse.cpp @@ -13,141 +13,141 @@ std::string inputs_name = ""; int main(int argc, char* argv[]) { - amrex::Initialize(argc, argv, false); + amrex::Initialize(argc, argv, false); - // timer for profiling - BL_PROFILE_VAR("main()", pmain); + // timer for profiling + BL_PROFILE_VAR("main()", pmain); - // abort if dim != 2 - if (AMREX_SPACEDIM != 2) - Abort("ERROR: must compile gaussian pulse diagnostic with DIM=2"); + // abort if dim != 2 + if (AMREX_SPACEDIM != 2) + Abort("ERROR: must compile gaussian pulse diagnostic with DIM=2"); - // Input arguments - string pltfile, slcfile; - int dir; + // Input arguments + string pltfile, slcfile; + int dir; - GetInputArgs(argc, argv, pltfile, slcfile dir); + GetInputArgs(argc, argv, pltfile, slcfile dir); auto center = GetCenter(pltfile); - double xctr = center[0]; - double yctr = center[1]; + double xctr = center[0]; + double yctr = center[1]; - Print() << "xctr = " << xctr << std::endl; - Print() << "yctr = " << yctr << std::endl; - Print() << std::endl; + Print() << "xctr = " << xctr << std::endl; + Print() << "yctr = " << yctr << std::endl; + Print() << std::endl; - // Start dataservices - DataServices::SetBatchMode(); + // Start dataservices + DataServices::SetBatchMode(); - // Define the type of file - Amrvis::FileType fileType(Amrvis::NEWPLT); - DataServices dataServices (pltfile, fileType); + // Define the type of file + Amrvis::FileType fileType(Amrvis::NEWPLT); + DataServices dataServices (pltfile, fileType); - if (!dataServices.AmrDataOk()) - DataServices::Dispatch(DataServices::ExitRequest, NULL); + if (!dataServices.AmrDataOk()) + DataServices::Dispatch(DataServices::ExitRequest, NULL); - // get data from plot file - AmrData& data = dataServices.AmrDataRef(); + // get data from plot file + AmrData& data = dataServices.AmrDataRef(); - int finestLevel = data.FinestLevel(); + int finestLevel = data.FinestLevel(); - // get variable names - const Vector& varNames = data.PlotVarNames(); + // get variable names + const Vector& varNames = data.PlotVarNames(); - // get the index bounds and dx. - Box domain = data.ProbDomain()[finestLevel]; - Vector dx = data.CellSize(finestLevel); - const Vector& problo = data.ProbLo(); - const Vector& probhi = data.ProbHi(); + // get the index bounds and dx. + Box domain = data.ProbDomain()[finestLevel]; + Vector dx = data.CellSize(finestLevel); + const Vector& problo = data.ProbLo(); + const Vector& probhi = data.ProbHi(); - Vector rr = data.RefRatio(); + Vector rr = data.RefRatio(); - // compute the size of the radially-binned array -- we'll do it to - // the furtherest corner of the domain - double x_maxdist = max(fabs(probhi[0] - xctr), fabs(problo[0] - xctr)); - double y_maxdist = max(fabs(probhi[1] - yctr), fabs(problo[1] - yctr)); - double maxdist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist); + // compute the size of the radially-binned array -- we'll do it to + // the furtherest corner of the domain + double x_maxdist = max(fabs(probhi[0] - xctr), fabs(problo[0] - xctr)); + double y_maxdist = max(fabs(probhi[1] - yctr), fabs(problo[1] - yctr)); + double maxdist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist); - double dx_fine = *(std::min_element(dx.begin(), dx.end())); + double dx_fine = *(std::min_element(dx.begin(), dx.end())); - int nbins = int(maxdist / dx_fine); + int nbins = int(maxdist / dx_fine); - // radial coordinate - Vector r(nbins); - for (auto i = 0; i < nbins; i++) - r[i] = (i + 0.5) * dx_fine; + // radial coordinate + Vector r(nbins); + for (auto i = 0; i < nbins; i++) + r[i] = (i + 0.5) * dx_fine; - // find variable indices - Vector slcvarNames = {"rad"}; - const auto nvars = slcvarNames.size(); + // find variable indices + Vector slcvarNames = {"rad"}; + const auto nvars = slcvarNames.size(); - auto varComps = GetComponents(data, slcvarNames); - auto rad_comp = varComps[0]; + auto varComps = GetComponents(data, slcvarNames); + auto rad_comp = varComps[0]; - // allocate storage for data - Vector rad_bin(nbins, 0.); - Vector ncount(nbins, 0); + // allocate storage for data + Vector rad_bin(nbins, 0.); + Vector ncount(nbins, 0); - // r1 is the factor between the current level grid spacing and the - // FINEST level - auto r1 = 1.0; + // r1 is the factor between the current level grid spacing and the + // FINEST level + auto r1 = 1.0; - // fill a multifab with the data - Vector fill_comps(data.NComp()); - for (auto i = 0; i < data.NComp(); i++) - fill_comps[i] = i; + // fill a multifab with the data + Vector fill_comps(data.NComp()); + for (auto i = 0; i < data.NComp(); i++) + fill_comps[i] = i; - // imask will be set to false if we've already output the data. - // Note, imask is defined in terms of the finest level. As we loop - // over levels, we will compare to the finest level index space to - // determine if we've already output here - int mask_size = domain.length().max(); - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + // imask will be set to false if we've already output the data. + // Note, imask is defined in terms of the finest level. As we loop + // over levels, we will compare to the finest level index space to + // determine if we've already output here + int mask_size = domain.length().max(); + Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); - // loop over the data, starting at the finest grid, and if we haven't - // already stored data in that grid location (according to imask), - // store it. - for (int l = finestLevel; l >= 0; l--) { + // loop over the data, starting at the finest grid, and if we haven't + // already stored data in that grid location (according to imask), + // store it. + for (int l = finestLevel; l >= 0; l--) { - Vector level_dx = data.DxLevel()[l]; + Vector level_dx = data.DxLevel()[l]; - const BoxArray& ba = data.boxArray(l); - const DistributionMapping& dm = data.DistributionMap(l); + const BoxArray& ba = data.boxArray(l); + const DistributionMapping& dm = data.DistributionMap(l); - MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); - data.FillVar(lev_data_mf, l, varNames, fill_comps); + MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); + data.FillVar(lev_data_mf, l, varNames, fill_comps); #ifdef _OPENMP #pragma omp parallel #endif - for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); + for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); - fgaussian_pulse(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, rad_bin.dataPtr(), ncount.dataPtr(), - imask.dataPtr(), mask_size, r1, - rad_comp, ZFILL(dx), dx_fine, xctr, yctr); + fgaussian_pulse(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), + nbins, rad_bin.dataPtr(), ncount.dataPtr(), + imask.dataPtr(), mask_size, r1, + rad_comp, ZFILL(dx), dx_fine, xctr, yctr); - } + } - // adjust r1 for the next lowest level - if (l != 0) r1 *= rr[l-1]; - } + // adjust r1 for the next lowest level + if (l != 0) r1 *= rr[l-1]; + } - // normalize - for (int i = 0; i < nbins; i++) - if (ncount[i] != 0) - rad_bin[i] /= ncount[i]; + // normalize + for (int i = 0; i < nbins; i++) + if (ncount[i] != 0) + rad_bin[i] /= ncount[i]; - // write data to slicefile - Vector > vars(nvars); - vars[0] = rad_bin; + // write data to slicefile + Vector > vars(nvars); + vars[0] = rad_bin; - WriteSlicefile(nbins, r, slcvarNames, vars, slcfile); + WriteSlicefile(nbins, r, slcvarNames, vars, slcfile); - // destroy timer for profiling - BL_PROFILE_VAR_STOP(pmain); + // destroy timer for profiling + BL_PROFILE_VAR_STOP(pmain); - amrex::Finalize(); + amrex::Finalize(); } diff --git a/Diagnostics/Radiation/lgt_frnt1d.cpp b/Diagnostics/Radiation/lgt_frnt1d.cpp index c84ae99043..55f61a5ffe 100644 --- a/Diagnostics/Radiation/lgt_frnt1d.cpp +++ b/Diagnostics/Radiation/lgt_frnt1d.cpp @@ -14,132 +14,132 @@ std::string inputs_name = ""; int main(int argc, char* argv[]) { - amrex::Initialize(argc, argv, false); + amrex::Initialize(argc, argv, false); - // timer for profiling - BL_PROFILE_VAR("main()", pmain); + // timer for profiling + BL_PROFILE_VAR("main()", pmain); - // Input arguments - string pltfile, slcfile; - int dir; + // Input arguments + string pltfile, slcfile; + int dir; - GetInputArgs(argc, argv, pltfile, slcfile, dir); + GetInputArgs(argc, argv, pltfile, slcfile, dir); - // Start dataservices - DataServices::SetBatchMode(); + // Start dataservices + DataServices::SetBatchMode(); - // Define the type of file - Amrvis::FileType fileType(Amrvis::NEWPLT); - DataServices dataServices (pltfile, fileType); + // Define the type of file + Amrvis::FileType fileType(Amrvis::NEWPLT); + DataServices dataServices (pltfile, fileType); - if (!dataServices.AmrDataOk()) - DataServices::Dispatch(DataServices::ExitRequest, NULL); + if (!dataServices.AmrDataOk()) + DataServices::Dispatch(DataServices::ExitRequest, NULL); - // get data from plot file - AmrData& data = dataServices.AmrDataRef(); + // get data from plot file + AmrData& data = dataServices.AmrDataRef(); - int finestLevel = data.FinestLevel(); + int finestLevel = data.FinestLevel(); - // get variable names - const Vector& varNames = data.PlotVarNames(); + // get variable names + const Vector& varNames = data.PlotVarNames(); - // get the index bounds and dx. - Box domain = data.ProbDomain()[finestLevel]; - Vector dx = data.CellSize(finestLevel); - const Vector& problo = data.ProbLo(); - const Vector& probhi = data.ProbHi(); + // get the index bounds and dx. + Box domain = data.ProbDomain()[finestLevel]; + Vector dx = data.CellSize(finestLevel); + const Vector& problo = data.ProbLo(); + const Vector& probhi = data.ProbHi(); - Vector rr = data.RefRatio(); + Vector rr = data.RefRatio(); - // compute the size of the radially-binned array -- we'll do it to - // the furtherest corner of the domain - double maxdist = fabs(probhi[0] - problo[0]); - double dx_fine = *(std::min_element(dx.begin(), dx.end())); - int nbins = int(maxdist / dx_fine); + // compute the size of the radially-binned array -- we'll do it to + // the furtherest corner of the domain + double maxdist = fabs(probhi[0] - problo[0]); + double dx_fine = *(std::min_element(dx.begin(), dx.end())); + int nbins = int(maxdist / dx_fine); - // radial coordinate - Vector r(nbins); + // radial coordinate + Vector r(nbins); - for (auto i = 0; i < nbins; i++) - r[i] = (i + 0.5) * dx_fine; + for (auto i = 0; i < nbins; i++) + r[i] = (i + 0.5) * dx_fine; - // find variable indices - Vector compVarNames = {"density", "xmom", "pressure", "rad"}; + // find variable indices + Vector compVarNames = {"density", "xmom", "pressure", "rad"}; - auto varComps = GetComponents(data, compVarNames); - auto dens_comp = varComps[0]; - auto xmom_comp = varComps[1]; - auto pres_comp = varComps[2]; - auto rad_comp = varComps[3]; + auto varComps = GetComponents(data, compVarNames); + auto dens_comp = varComps[0]; + auto xmom_comp = varComps[1]; + auto pres_comp = varComps[2]; + auto rad_comp = varComps[3]; - // allocate storage for data - Vector dens_bin(nbins, 0.); - Vector vel_bin(nbins, 0.); - Vector pres_bin(nbins, 0.); - Vector rad_bin(nbins, 0.); - Vector ncount(nbins, 0); + // allocate storage for data + Vector dens_bin(nbins, 0.); + Vector vel_bin(nbins, 0.); + Vector pres_bin(nbins, 0.); + Vector rad_bin(nbins, 0.); + Vector ncount(nbins, 0); - // r1 is the factor between the current level grid spacing and the - // FINEST level - auto r1 = 1.0; + // r1 is the factor between the current level grid spacing and the + // FINEST level + auto r1 = 1.0; - // fill a multifab with the data - Vector fill_comps(data.NComp()); - for (auto i = 0; i < data.NComp(); i++) - fill_comps[i] = i; + // fill a multifab with the data + Vector fill_comps(data.NComp()); + for (auto i = 0; i < data.NComp(); i++) + fill_comps[i] = i; - // imask will be set to false if we've already output the data. - // Note, imask is defined in terms of the finest level. As we loop - // over levels, we will compare to the finest level index space to - // determine if we've already output here - int mask_size = domain.length().max(); - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + // imask will be set to false if we've already output the data. + // Note, imask is defined in terms of the finest level. As we loop + // over levels, we will compare to the finest level index space to + // determine if we've already output here + int mask_size = domain.length().max(); + Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); - // loop over the data, starting at the finest grid, and if we haven't - // already stored data in that grid location (according to imask), - // store it. - for (int l = finestLevel; l >= 0; l--) { + // loop over the data, starting at the finest grid, and if we haven't + // already stored data in that grid location (according to imask), + // store it. + for (int l = finestLevel; l >= 0; l--) { - Vector level_dx = data.DxLevel()[l]; + Vector level_dx = data.DxLevel()[l]; - const BoxArray& ba = data.boxArray(l); - const DistributionMapping& dm = data.DistributionMap(l); + const BoxArray& ba = data.boxArray(l); + const DistributionMapping& dm = data.DistributionMap(l); - MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); - data.FillVar(lev_data_mf, l, varNames, fill_comps); + MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); + data.FillVar(lev_data_mf, l, varNames, fill_comps); #ifdef _OPENMP #pragma omp parallel #endif - for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); + for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); - flgt_frnt1d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens_bin.dataPtr(), vel_bin.dataPtr(), - pres_bin.dataPtr(), rad_bin.dataPtr(), - imask.dataPtr(), mask_size, r1, - dens_comp, xmom_comp, pres_comp, rad_comp, - ZFILL(dx), dx_fine); - } + flgt_frnt1d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), + nbins, dens_bin.dataPtr(), vel_bin.dataPtr(), + pres_bin.dataPtr(), rad_bin.dataPtr(), + imask.dataPtr(), mask_size, r1, + dens_comp, xmom_comp, pres_comp, rad_comp, + ZFILL(dx), dx_fine); + } - // adjust r1 for the next lowest level - if (l != 0) r1 *= rr[l-1]; - } + // adjust r1 for the next lowest level + if (l != 0) r1 *= rr[l-1]; + } - //normalize - for (int i = 0; i < nbins; i++) - if (ncount[i] != 0) - rad_bin[i] /= ncount[i]; + //normalize + for (int i = 0; i < nbins; i++) + if (ncount[i] != 0) + rad_bin[i] /= ncount[i]; - // write data to slicefile - Vector > vars = {dens_bin, vel_bin, pres_bin, rad_bin}; - Vector slcvarNames = {"density", "velocity", "pressure", "rad"}; + // write data to slicefile + Vector > vars = {dens_bin, vel_bin, pres_bin, rad_bin}; + Vector slcvarNames = {"density", "velocity", "pressure", "rad"}; - WriteSlicefile(nbins, r, slcvarNames, vars, slcfile); + WriteSlicefile(nbins, r, slcvarNames, vars, slcfile); - // destroy timer for profiling - BL_PROFILE_VAR_STOP(pmain); + // destroy timer for profiling + BL_PROFILE_VAR_STOP(pmain); - amrex::Finalize(); + amrex::Finalize(); } diff --git a/Diagnostics/Radiation/rad_shock.cpp b/Diagnostics/Radiation/rad_shock.cpp index 37ad5ceb0c..43176c4e27 100644 --- a/Diagnostics/Radiation/rad_shock.cpp +++ b/Diagnostics/Radiation/rad_shock.cpp @@ -20,209 +20,209 @@ std::string inputs_name = ""; int main(int argc, char* argv[]) { - amrex::Initialize(argc, argv, false); + amrex::Initialize(argc, argv, false); - // timer for profiling - BL_PROFILE_VAR("main()", pmain); + // timer for profiling + BL_PROFILE_VAR("main()", pmain); - // Input arguments - string pltfile; - string slcfile; - int idir = 1; + // Input arguments + string pltfile; + string slcfile; + int idir = 1; - GetInputArgs (argc, argv, pltfile, slcfile, idir); + GetInputArgs (argc, argv, pltfile, slcfile, idir); - Print() << "idir = " << idir << std::endl; + Print() << "idir = " << idir << std::endl; - // check that idir <= DIM - if (idir > AMREX_SPACEDIM) - Abort("ERROR: idir must be <= DIM"); + // check that idir <= DIM + if (idir > AMREX_SPACEDIM) + Abort("ERROR: idir must be <= DIM"); - // Start dataservices - DataServices::SetBatchMode(); + // Start dataservices + DataServices::SetBatchMode(); - // Define the type of file - Amrvis::FileType fileType(Amrvis::NEWPLT); - DataServices dataServices (pltfile, fileType); + // Define the type of file + Amrvis::FileType fileType(Amrvis::NEWPLT); + DataServices dataServices (pltfile, fileType); - if (!dataServices.AmrDataOk()) - DataServices::Dispatch(DataServices::ExitRequest, NULL); + if (!dataServices.AmrDataOk()) + DataServices::Dispatch(DataServices::ExitRequest, NULL); - // get data from plot file - AmrData& data = dataServices.AmrDataRef(); + // get data from plot file + AmrData& data = dataServices.AmrDataRef(); - int finestLevel = data.FinestLevel(); + int finestLevel = data.FinestLevel(); - // get variable names - const Vector& varNames = data.PlotVarNames(); + // get variable names + const Vector& varNames = data.PlotVarNames(); - // get the index bounds and dx. - Box domain = data.ProbDomain()[finestLevel]; - Vector dx = data.CellSize(finestLevel); - Vector problo = data.ProbLo(); - Vector probhi = data.ProbHi(); + // get the index bounds and dx. + Box domain = data.ProbDomain()[finestLevel]; + Vector dx = data.CellSize(finestLevel); + Vector problo = data.ProbLo(); + Vector probhi = data.ProbHi(); - Vector rr = data.RefRatio(); + Vector rr = data.RefRatio(); - // compute the size of the radially-binned array - int nbins = domain.length()[idir-1]; + // compute the size of the radially-binned array + int nbins = domain.length()[idir-1]; - // radial coordinate - Vector r(nbins); + // radial coordinate + Vector r(nbins); - for (auto i = 0; i < nbins; i++) - r[i] = (i + 0.5) * dx[idir-1] + problo[idir-1]; + for (auto i = 0; i < nbins; i++) + r[i] = (i + 0.5) * dx[idir-1] + problo[idir-1]; - // find variable indices + // find variable indices #if (AMREX_SPACEDIM == 1) - Vector compVarNames = {"density", "eint_E", "Temp", - "pressure", "rad", "x_velocity"}; + Vector compVarNames = {"density", "eint_E", "Temp", + "pressure", "rad", "x_velocity"}; #elif (AMREX_SPACEDIM == 2) - Vector compVarNames = {"density", - "eint_E", "Temp", "pressure", "rad", - "x_velocity", "y_velocity"}; + Vector compVarNames = {"density", + "eint_E", "Temp", "pressure", "rad", + "x_velocity", "y_velocity"}; #else - Vector compVarNames = {"density", "eint_E", "Temp", "pressure", - "rad", "x_velocity", "y_velocity","z_velocity"}; + Vector compVarNames = {"density", "eint_E", "Temp", "pressure", + "rad", "x_velocity", "y_velocity","z_velocity"}; #endif - auto varComps = GetComponents(data, compVarNames); + auto varComps = GetComponents(data, compVarNames); - // r1 is the factor between the current level grid spacing and the - // FINEST level - auto r1 = 1.0; + // r1 is the factor between the current level grid spacing and the + // FINEST level + auto r1 = 1.0; - // number of variables in plotfile - const int nvars = data.NComp(); + // number of variables in plotfile + const int nvars = data.NComp(); - Vector fill_comps(nvars); - for (auto i = 0; i < nvars; i++) - fill_comps[i] = i; + Vector fill_comps(nvars); + for (auto i = 0; i < nvars; i++) + fill_comps[i] = i; - // allocate storage for data - Vector vars_bin(nbins * (nvars + 1), 0.); + // allocate storage for data + Vector vars_bin(nbins * (nvars + 1), 0.); - // imask will be set to false if we've already output the data. - // Note, imask is defined in terms of the finest level. As we loop - // over levels, we will compare to the finest level index space to - // determine if we've already output here - int mask_size = domain.length().max(); - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + // imask will be set to false if we've already output the data. + // Note, imask is defined in terms of the finest level. As we loop + // over levels, we will compare to the finest level index space to + // determine if we've already output here + int mask_size = domain.length().max(); + Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); - // counter - auto cnt = 0; + // counter + auto cnt = 0; - // loop over the data, starting at the finest grid, and if we haven't - // already stored data in that grid location (according to imask), - // store it. - for (int l = finestLevel; l >= 0; l--) { + // loop over the data, starting at the finest grid, and if we haven't + // already stored data in that grid location (according to imask), + // store it. + for (int l = finestLevel; l >= 0; l--) { - // refratio is the factor between the COARSEST level grid spacing and + // refratio is the factor between the COARSEST level grid spacing and // the current level - int refratio = 1; - for (auto lev = 0; lev < l; lev++) refratio *= rr[lev]; + int refratio = 1; + for (auto lev = 0; lev < l; lev++) refratio *= rr[lev]; - Vector level_dx = data.DxLevel()[l]; + Vector level_dx = data.DxLevel()[l]; - const BoxArray& ba = data.boxArray(l); - const DistributionMapping& dm = data.DistributionMap(l); + const BoxArray& ba = data.boxArray(l); + const DistributionMapping& dm = data.DistributionMap(l); - MultiFab lev_data_mf(ba, dm, nvars, data.NGrow()); - data.FillVar(lev_data_mf, l, varNames, fill_comps); + MultiFab lev_data_mf(ba, dm, nvars, data.NGrow()); + data.FillVar(lev_data_mf, l, varNames, fill_comps); #ifdef _OPENMP #pragma omp parallel #endif - for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); + for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); - fradshock(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - problo.dataPtr(), probhi.dataPtr(), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, vars_bin.dataPtr(), - imask.dataPtr(), mask_size, r1, refratio, - dx.dataPtr(), idir, &cnt); + fradshock(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + problo.dataPtr(), probhi.dataPtr(), + BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), + nbins, vars_bin.dataPtr(), + imask.dataPtr(), mask_size, r1, refratio, + dx.dataPtr(), idir, &cnt); - } + } - // adjust r1 for the next lowest level - if (l != 0) r1 *= rr[l-1]; - } + // adjust r1 for the next lowest level + if (l != 0) r1 *= rr[l-1]; + } - // sort the data based on the coordinates - Vector coords(cnt); - for (auto i = 0; i < cnt; i++) - coords[i] = vars_bin[i]; + // sort the data based on the coordinates + Vector coords(cnt); + for (auto i = 0; i < cnt; i++) + coords[i] = vars_bin[i]; - auto isv = sort_indexes(coords); + auto isv = sort_indexes(coords); - // find variable indices + // find variable indices #if (AMREX_SPACEDIM == 1) - compVarNames = {"density", "x_velocity", "pressure", - "eint_E", "Temp", "rad", "rad"}; + compVarNames = {"density", "x_velocity", "pressure", + "eint_E", "Temp", "rad", "rad"}; #elif (AMREX_SPACEDIM == 2) - compVarNames = {"density", "x_velocity", "y_velocity", "pressure", - "eint_E", "Temp", "rad", "rad"}; + compVarNames = {"density", "x_velocity", "y_velocity", "pressure", + "eint_E", "Temp", "rad", "rad"}; #else - compVarNames = {"density", "x_velocity", "y_velocity", "z_velocity", - "pressure","eint_E", "Temp", "rad", "rad"}; + compVarNames = {"density", "x_velocity", "y_velocity", "z_velocity", + "pressure","eint_E", "Temp", "rad", "rad"}; #endif - Vector comps = GetComponents(data, compVarNames); + Vector comps = GetComponents(data, compVarNames); #if (AMREX_SPACEDIM == 1) - Vector slcvarNames = {"density", "x-velocity", - "pressure", "int. energy", "temperature", - "rad energy", "rad temp"}; + Vector slcvarNames = {"density", "x-velocity", + "pressure", "int. energy", "temperature", + "rad energy", "rad temp"}; #elif (AMREX_SPACEDIM == 2) - Vector slcvarNames = {"density", "x-velocity", "y-velocity", - "pressure", "int. energy", "temperature", - "rad energy", "rad temp"}; + Vector slcvarNames = {"density", "x-velocity", "y-velocity", + "pressure", "int. energy", "temperature", + "rad energy", "rad temp"}; #else - Vector slcvarNames = {"density", "x-velocity", "y-velocity", - "z-velocity", "pressure", "int. energy", "temperature", "rad energy", "rad temp"}; + Vector slcvarNames = {"density", "x-velocity", "y-velocity", + "z-velocity", "pressure", "int. energy", "temperature", "rad energy", "rad temp"}; #endif - // write to file - std::ofstream slicefile; - slicefile.open(slcfile); - slicefile.setf(std::ios::scientific); - slicefile.precision(12); - const auto w = 24; + // write to file + std::ofstream slicefile; + slicefile.open(slcfile); + slicefile.setf(std::ios::scientific); + slicefile.precision(12); + const auto w = 24; - // write the header - if (idir == 1) - slicefile << std::setw(w) << "x"; - else if (idir == 2) - slicefile << std::setw(w) << "y"; - else - slicefile << std::setw(w) << "z"; + // write the header + if (idir == 1) + slicefile << std::setw(w) << "x"; + else if (idir == 2) + slicefile << std::setw(w) << "y"; + else + slicefile << std::setw(w) << "z"; - for (auto it=slcvarNames.begin(); it!=slcvarNames.end(); ++it) - slicefile << std::setw(w) << *it; + for (auto it=slcvarNames.begin(); it!=slcvarNames.end(); ++it) + slicefile << std::setw(w) << *it; - slicefile << std::endl; + slicefile << std::endl; - // NOTE: I could not find the constant arad anywhere, so I'm - // setting it to 1 here. - const Real arad = 1.0; + // NOTE: I could not find the constant arad anywhere, so I'm + // setting it to 1 here. + const Real arad = 1.0; - // write the data in columns - for (auto i = 0; i < cnt; i++) { + // write the data in columns + for (auto i = 0; i < cnt; i++) { - slicefile << std::setw(w) << vars_bin[isv[i]]; + slicefile << std::setw(w) << vars_bin[isv[i]]; - for (auto it=comps.begin(); it!=(comps.end()-1); ++it) - slicefile << std::setw(w) << vars_bin[isv[i]+((*it)+1)*nbins]; + for (auto it=comps.begin(); it!=(comps.end()-1); ++it) + slicefile << std::setw(w) << vars_bin[isv[i]+((*it)+1)*nbins]; - auto it = comps.end()-1; - slicefile << std::setw(w) << pow(vars_bin[isv[i]+((*it)+1)*nbins] / arad, 0.25); + auto it = comps.end()-1; + slicefile << std::setw(w) << pow(vars_bin[isv[i]+((*it)+1)*nbins] / arad, 0.25); - slicefile << std::endl; - } + slicefile << std::endl; + } - slicefile.close(); + slicefile.close(); - // destroy timer for profiling - BL_PROFILE_VAR_STOP(pmain); + // destroy timer for profiling + BL_PROFILE_VAR_STOP(pmain); - amrex::Finalize(); + amrex::Finalize(); } diff --git a/Diagnostics/Radiation/rad_source.cpp b/Diagnostics/Radiation/rad_source.cpp index eb890cd198..fae6885035 100644 --- a/Diagnostics/Radiation/rad_source.cpp +++ b/Diagnostics/Radiation/rad_source.cpp @@ -20,87 +20,87 @@ std::string inputs_name = ""; int main(int argc, char* argv[]) { - amrex::Initialize(argc, argv, false); + amrex::Initialize(argc, argv, false); - // timer for profiling - BL_PROFILE_VAR("main()", pmain); + // timer for profiling + BL_PROFILE_VAR("main()", pmain); - // check DIM = 1 - if (AMREX_SPACEDIM != 1) - Abort("ERROR: rad_source diagnostic can only be compiled with DIM=1"); + // check DIM = 1 + if (AMREX_SPACEDIM != 1) + Abort("ERROR: rad_source diagnostic can only be compiled with DIM=1"); - // Input arguments - if (argc < 2) - Abort("ERROR: Missing plotfiles"); + // Input arguments + if (argc < 2) + Abort("ERROR: Missing plotfiles"); - // all the arguments are assumed to be plotfiles - // loop over all the plotfiles - for (auto i = 1; i < argc; i++) { + // all the arguments are assumed to be plotfiles + // loop over all the plotfiles + for (auto i = 1; i < argc; i++) { - string pltfile = argv[i]; + string pltfile = argv[i]; - // Start dataservices - DataServices::SetBatchMode(); + // Start dataservices + DataServices::SetBatchMode(); - // Define the type of file - Amrvis::FileType fileType(Amrvis::NEWPLT); - DataServices dataServices (pltfile, fileType); + // Define the type of file + Amrvis::FileType fileType(Amrvis::NEWPLT); + DataServices dataServices (pltfile, fileType); - if (!dataServices.AmrDataOk()) - DataServices::Dispatch(DataServices::ExitRequest, NULL); + if (!dataServices.AmrDataOk()) + DataServices::Dispatch(DataServices::ExitRequest, NULL); - // get data from plot file - AmrData& data = dataServices.AmrDataRef(); + // get data from plot file + AmrData& data = dataServices.AmrDataRef(); - // get variable names - const Vector& varNames = data.PlotVarNames(); + // get variable names + const Vector& varNames = data.PlotVarNames(); - // find variable indices - Vector compVarNames = {"rho_e", "rad"}; + // find variable indices + Vector compVarNames = {"rho_e", "rad"}; - auto varComps = GetComponents(data, compVarNames); - auto rhoe_comp = varComps[0]; - auto rad_comp = varComps[1]; + auto varComps = GetComponents(data, compVarNames); + auto rhoe_comp = varComps[0]; + auto rad_comp = varComps[1]; - Vector fill_comps(data.NComp()); - for (auto j = 0; j < data.NComp(); j++) - fill_comps[j] = j; + Vector fill_comps(data.NComp()); + for (auto j = 0; j < data.NComp(); j++) + fill_comps[j] = j; - // only work with the finest level's data - auto lev = data.FinestLevel();; + // only work with the finest level's data + auto lev = data.FinestLevel();; - const BoxArray& ba = data.boxArray(lev); - const DistributionMapping& dm = data.DistributionMap(lev); + const BoxArray& ba = data.boxArray(lev); + const DistributionMapping& dm = data.DistributionMap(lev); - MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); - data.FillVar(lev_data_mf, lev, varNames, fill_comps); + MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); + data.FillVar(lev_data_mf, lev, varNames, fill_comps); - // we only care about a single zone, so take the first box - MFIter mfi(lev_data_mf); - const Box& bx = mfi.tilebox(); + // we only care about a single zone, so take the first box + MFIter mfi(lev_data_mf); + const Box& bx = mfi.tilebox(); - Real rhoe, rad; + Real rhoe, rad; - fradsource(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - &rhoe, &rad, rhoe_comp, rad_comp); + fradsource(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), + &rhoe, &rad, rhoe_comp, rad_comp); - const auto w = 20; + const auto w = 20; - std::cout.setf(std::ios::scientific); - std::cout.precision(12); + std::cout.setf(std::ios::scientific); + std::cout.precision(12); - if (i == 1) - std::cout << std::setw(w) << "time" << std::setw(w) << "rho e" - << std::setw(w) << "rad" << std::endl; + if (i == 1) + std::cout << std::setw(w) << "time" << std::setw(w) << "rho e" + << std::setw(w) << "rad" << std::endl; - std::cout << std::setw(w) << data.Time() << std::setw(w) << rhoe << std::setw(w) - << rad << std::endl; + std::cout << std::setw(w) << data.Time() << std::setw(w) << rhoe << std::setw(w) + << rad << std::endl; - } + } - // destroy timer for profiling - BL_PROFILE_VAR_STOP(pmain); + // destroy timer for profiling + BL_PROFILE_VAR_STOP(pmain); - amrex::Finalize(); + amrex::Finalize(); } diff --git a/Diagnostics/Radiation/rad_sphere.cpp b/Diagnostics/Radiation/rad_sphere.cpp index 64bf6483fb..8eda551489 100644 --- a/Diagnostics/Radiation/rad_sphere.cpp +++ b/Diagnostics/Radiation/rad_sphere.cpp @@ -13,256 +13,256 @@ using namespace amrex; std::string inputs_name = ""; void Print_Help() { - Print() << "\nPrint out the radiation quantities at a specified distance from" - << "\nthe origin. This is written for the 1-d radiating sphere problem." - << "\n" - << "\n./fradsphere -p plotfile -r radius -g groupfile -v variable" - << "\n" - << "\nHere groupfile is the file containing the group structure information" - << "\nas output by Castro (usually group_structure.dat)." - << "\nvariable should be either rad or rad_analytic_ (to obtain the numerical" - << "\nand analytic data, respectively.)" - << "\n\n" << std::endl; + Print() << "\nPrint out the radiation quantities at a specified distance from" + << "\nthe origin. This is written for the 1-d radiating sphere problem." + << "\n" + << "\n./fradsphere -p plotfile -r radius -g groupfile -v variable" + << "\n" + << "\nHere groupfile is the file containing the group structure information" + << "\nas output by Castro (usually group_structure.dat)." + << "\nvariable should be either rad or rad_analytic_ (to obtain the numerical" + << "\nand analytic data, respectively.)" + << "\n\n" << std::endl; } int main(int argc, char* argv[]) { - amrex::Initialize(argc, argv, false); - - // timer for profiling - BL_PROFILE_VAR("main()", pmain); - - if (AMREX_SPACEDIM != 1) - Abort("ERROR: rad_sphere diagnostic only works for DIM=1"); - - // Input arguments - string pltfile, groupfile, variable; - Real radius = 0.; - int j = 1; // skip program name - - while ( j < argc) - { - - if ( !strcmp(argv[j], "-p") || !strcmp(argv[j],"--pltfile") ) - { - pltfile = argv[++j]; - } - else if ( !strcmp(argv[j], "-g") || !strcmp(argv[j],"--groupfile") ) - { - groupfile = argv[++j]; - } - else if ( !strcmp(argv[j], "-r") || !strcmp(argv[j],"--radius") ) - { - radius = std::atof(argv[++j]); - } - else if ( !strcmp(argv[j], "-v") || !strcmp(argv[j],"--variable") ) - { - variable = argv[++j]; - } - else - { - std::cout << "\n\nOption " << argv[j] << " not recognized" << std::endl; - Print_Help(); - exit ( EXIT_FAILURE ); - } - - // Go to the next parameter name - ++j; - } - - if (pltfile.empty() || groupfile.empty()) - { - Print_Help(); - Abort("Missing input file(s)"); - } - - Print() << "\nplotfile = \"" << pltfile << "\"" << std::endl; - Print() << "groupfile = \"" << groupfile << "\"" << std::endl; - Print() << "radius = " << radius << std::endl; - Print() << "variable = " << variable << std::endl; - Print() << std::endl; - - // Start dataservices - DataServices::SetBatchMode(); - - // Define the type of file - Amrvis::FileType fileType(Amrvis::NEWPLT); - DataServices dataServices (pltfile, fileType); - - if (!dataServices.AmrDataOk()) - DataServices::Dispatch(DataServices::ExitRequest, NULL); - - // get data from plot file - AmrData& data = dataServices.AmrDataRef(); - - int finestLevel = data.FinestLevel(); - - // get variable names - const Vector& varNames = data.PlotVarNames(); - - // get the index bounds and dx. - Box domain = data.ProbDomain()[finestLevel]; - Vector dx = data.CellSize(finestLevel); - const Vector& problo = data.ProbLo(); - const Vector& probhi = data.ProbHi(); - - Vector rr = data.RefRatio(); - - if (radius < problo[0] || radius > probhi[0]) - Abort("ERROR: specified observer radius outside of domain"); - - std::cout.setf(std::ios::scientific); - std::cout.precision(12); - std::cout << "rmin = " << problo[0] << std::endl; - std::cout << "rmax = " << probhi[0] << std::endl << std::endl; - - int nbins = domain.length(0); - - // find variable indices - Vector compVarNames = {variable + "0"}; - - auto varComps = GetComponents(data, compVarNames); - auto rad_comp = varComps[0]; - - const int nvars = data.NComp(); - - // allocate storage for data - Vector vars_bin(nbins * (nvars + 1), 0.); - - // r1 is the factor between the current level grid spacing and the - // FINEST level - auto r1 = 1.0; - - Vector fill_comps(data.NComp()); - for (auto i = 0; i < data.NComp(); i++) - fill_comps[i] = i; - - // imask will be set to false if we've already output the data. - // Note, imask is defined in terms of the finest level. As we loop - // over levels, we will compare to the finest level index space to - // determine if we've already output here - int mask_size = domain.length()[0]; - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); - - // counter - int cnt = 0; - - // loop over the data, starting at the finest grid, and if we haven't - // already stored data in that grid location (according to imask), - // store it. - for (int l = finestLevel; l >= 0; l--) { + amrex::Initialize(argc, argv, false); - Vector level_dx = data.CellSize(l); + // timer for profiling + BL_PROFILE_VAR("main()", pmain); - const BoxArray& ba = data.boxArray(l); - const DistributionMapping& dm = data.DistributionMap(l); + if (AMREX_SPACEDIM != 1) + Abort("ERROR: rad_sphere diagnostic only works for DIM=1"); - MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); - data.FillVar(lev_data_mf, l, varNames, fill_comps); + // Input arguments + string pltfile, groupfile, variable; + Real radius = 0.; + int j = 1; // skip program name - for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); + while ( j < argc) + { - fradsphere(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - ZFILL(problo), ZFILL(probhi), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, vars_bin.dataPtr(), - imask.dataPtr(), mask_size, r1, - ZFILL(dx), &cnt); - } + if ( !strcmp(argv[j], "-p") || !strcmp(argv[j],"--pltfile") ) + { + pltfile = argv[++j]; + } + else if ( !strcmp(argv[j], "-g") || !strcmp(argv[j],"--groupfile") ) + { + groupfile = argv[++j]; + } + else if ( !strcmp(argv[j], "-r") || !strcmp(argv[j],"--radius") ) + { + radius = std::atof(argv[++j]); + } + else if ( !strcmp(argv[j], "-v") || !strcmp(argv[j],"--variable") ) + { + variable = argv[++j]; + } + else + { + std::cout << "\n\nOption " << argv[j] << " not recognized" << std::endl; + Print_Help(); + exit ( EXIT_FAILURE ); + } + + // Go to the next parameter name + ++j; + } + + if (pltfile.empty() || groupfile.empty()) + { + Print_Help(); + Abort("Missing input file(s)"); + } + + Print() << "\nplotfile = \"" << pltfile << "\"" << std::endl; + Print() << "groupfile = \"" << groupfile << "\"" << std::endl; + Print() << "radius = " << radius << std::endl; + Print() << "variable = " << variable << std::endl; + Print() << std::endl; + + // Start dataservices + DataServices::SetBatchMode(); + + // Define the type of file + Amrvis::FileType fileType(Amrvis::NEWPLT); + DataServices dataServices (pltfile, fileType); + + if (!dataServices.AmrDataOk()) + DataServices::Dispatch(DataServices::ExitRequest, NULL); + + // get data from plot file + AmrData& data = dataServices.AmrDataRef(); + + int finestLevel = data.FinestLevel(); + + // get variable names + const Vector& varNames = data.PlotVarNames(); + + // get the index bounds and dx. + Box domain = data.ProbDomain()[finestLevel]; + Vector dx = data.CellSize(finestLevel); + const Vector& problo = data.ProbLo(); + const Vector& probhi = data.ProbHi(); + + Vector rr = data.RefRatio(); + + if (radius < problo[0] || radius > probhi[0]) + Abort("ERROR: specified observer radius outside of domain"); + + std::cout.setf(std::ios::scientific); + std::cout.precision(12); + std::cout << "rmin = " << problo[0] << std::endl; + std::cout << "rmax = " << probhi[0] << std::endl << std::endl; - // adjust r1 for the next lowest level - if (l != 0) r1 *= rr[l-1]; - } + int nbins = domain.length(0); - // sort the data based on the coordinates - Vector coords(cnt); - for (auto i = 0; i < cnt; i++) - coords[i] = vars_bin[i]; + // find variable indices + Vector compVarNames = {variable + "0"}; - auto isv = sort_indexes(coords); + auto varComps = GetComponents(data, compVarNames); + auto rad_comp = varComps[0]; + + const int nvars = data.NComp(); + + // allocate storage for data + Vector vars_bin(nbins * (nvars + 1), 0.); + + // r1 is the factor between the current level grid spacing and the + // FINEST level + auto r1 = 1.0; + + Vector fill_comps(data.NComp()); + for (auto i = 0; i < data.NComp(); i++) + fill_comps[i] = i; + + // imask will be set to false if we've already output the data. + // Note, imask is defined in terms of the finest level. As we loop + // over levels, we will compare to the finest level index space to + // determine if we've already output here + int mask_size = domain.length()[0]; + Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + + // counter + int cnt = 0; + + // loop over the data, starting at the finest grid, and if we haven't + // already stored data in that grid location (according to imask), + // store it. + for (int l = finestLevel; l >= 0; l--) { + + Vector level_dx = data.CellSize(l); + + const BoxArray& ba = data.boxArray(l); + const DistributionMapping& dm = data.DistributionMap(l); + + MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); + data.FillVar(lev_data_mf, l, varNames, fill_comps); + + for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + + fradsphere(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + ZFILL(problo), ZFILL(probhi), + BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), + nbins, vars_bin.dataPtr(), + imask.dataPtr(), mask_size, r1, + ZFILL(dx), &cnt); + } + + // adjust r1 for the next lowest level + if (l != 0) r1 *= rr[l-1]; + } + + // sort the data based on the coordinates + Vector coords(cnt); + for (auto i = 0; i < cnt; i++) + coords[i] = vars_bin[i]; + + auto isv = sort_indexes(coords); Print() << "coords_min = " << coords[0] << " coords_max = " << coords[cnt-1] << std::endl; - // open the group file and read in the group information - std::ifstream group_file; - group_file.open(groupfile); + // open the group file and read in the group information + std::ifstream group_file; + group_file.open(groupfile); - string header_line; + string header_line; - string line; + string line; - Vector nu_groups; - Vector dnu_groups; + Vector nu_groups; + Vector dnu_groups; - bool first_line = true; + bool first_line = true; - int ngroups = 0; + int ngroups = 0; - while(std::getline(group_file, line)) { - if (first_line) { + while(std::getline(group_file, line)) { + if (first_line) { - const std::regex re("=\\s*([0-9]*)"); - std::smatch m; - std::regex_search(line, m, re); + const std::regex re("=\\s*([0-9]*)"); + std::smatch m; + std::regex_search(line, m, re); - ngroups = stoi(m[1].str()); + ngroups = stoi(m[1].str()); - first_line = false; + first_line = false; - // skip next line - std::getline(group_file, line); + // skip next line + std::getline(group_file, line); - } else { + } else { - // read in the group centers and weights - Real nu, dnu; - std::istringstream iss(line); - iss >> nu >> dnu; + // read in the group centers and weights + Real nu, dnu; + std::istringstream iss(line); + iss >> nu >> dnu; - nu_groups.push_back(nu); - dnu_groups.push_back(dnu); - } - } + nu_groups.push_back(nu); + dnu_groups.push_back(dnu); + } + } - group_file.close(); + group_file.close(); - // find the index corresponding to the desired observer radius - auto idx_obs = -1; + // find the index corresponding to the desired observer radius + auto idx_obs = -1; - for (auto i = 0; i < cnt; i++) { - if (radius >= vars_bin[isv[i]] && radius < vars_bin[isv[i+1]]) { - idx_obs = i; - break; + for (auto i = 0; i < cnt; i++) { + if (radius >= vars_bin[isv[i]] && radius < vars_bin[isv[i+1]]) { + idx_obs = i; + break; } - } + } - if (idx_obs == -1) Abort("ERROR: radius not found in domain"); + if (idx_obs == -1) Abort("ERROR: radius not found in domain"); - // output all the radiation energies - const auto w = 28; + // output all the radiation energies + const auto w = 28; - std::ofstream slicefile; - slicefile.open("rad_sphere.out"); - slicefile.setf(std::ios::scientific); - slicefile.precision(12); + std::ofstream slicefile; + slicefile.open("rad_sphere.out"); + slicefile.setf(std::ios::scientific); + slicefile.precision(12); - slicefile << std::setw(15) << "# group name" - << std::setw(w) << "group center energy" - << std::setw(w) << "E_rad(nu)*dnu (erg/cm^3)" - << std::setw(w) << "E_rad(nu) (erg/cm^3/Hz)" << std::endl; + slicefile << std::setw(15) << "# group name" + << std::setw(w) << "group center energy" + << std::setw(w) << "E_rad(nu)*dnu (erg/cm^3)" + << std::setw(w) << "E_rad(nu) (erg/cm^3/Hz)" << std::endl; - for (auto i = 0; i < ngroups; i++) { - slicefile << std::setw(15) << varNames[rad_comp+i] - << std::setw(w) << nu_groups[i] - << std::setw(w) << vars_bin[isv[idx_obs] + (rad_comp+i+1) * nbins] - << std::setw(w) << vars_bin[isv[idx_obs] + (rad_comp+i+1) * nbins] / dnu_groups[i] << std::endl; - } + for (auto i = 0; i < ngroups; i++) { + slicefile << std::setw(15) << varNames[rad_comp+i] + << std::setw(w) << nu_groups[i] + << std::setw(w) << vars_bin[isv[idx_obs] + (rad_comp+i+1) * nbins] + << std::setw(w) << vars_bin[isv[idx_obs] + (rad_comp+i+1) * nbins] / dnu_groups[i] << std::endl; + } - slicefile.close(); + slicefile.close(); - // destroy timer for profiling - BL_PROFILE_VAR_STOP(pmain); + // destroy timer for profiling + BL_PROFILE_VAR_STOP(pmain); - amrex::Finalize(); + amrex::Finalize(); } diff --git a/Diagnostics/Radiation/rhd_shocktube.cpp b/Diagnostics/Radiation/rhd_shocktube.cpp index 269f55a189..9ed79f00ca 100644 --- a/Diagnostics/Radiation/rhd_shocktube.cpp +++ b/Diagnostics/Radiation/rhd_shocktube.cpp @@ -14,157 +14,157 @@ std::string inputs_name = ""; int main(int argc, char* argv[]) { - amrex::Initialize(argc, argv, false); + amrex::Initialize(argc, argv, false); - // timer for profiling - BL_PROFILE_VAR("main()", pmain); + // timer for profiling + BL_PROFILE_VAR("main()", pmain); - if (AMREX_SPACEDIM != 1) - Abort("ERROR: rhd_shocktube diagnostic only works for DIM=1"); - - // Input arguments - string pltfile, groupfile, slcfile; - int j = 1; // skip program name - - while ( j < argc) - { + if (AMREX_SPACEDIM != 1) + Abort("ERROR: rhd_shocktube diagnostic only works for DIM=1"); + + // Input arguments + string pltfile, groupfile, slcfile; + int j = 1; // skip program name + + while ( j < argc) + { - if ( !strcmp(argv[j], "-p") || !strcmp(argv[j],"--pltfile") ) - { - pltfile = argv[++j]; - } - else if ( !strcmp(argv[j], "-g") || !strcmp(argv[j],"--groupfile") ) - { - groupfile = argv[++j]; - } - else if ( !strcmp(argv[j], "-s") || !strcmp(argv[j],"--slicefile") ) - { - slcfile = argv[++j]; - } - else - { - std::cout << "\n\nOption " << argv[j] << " not recognized" << std::endl; - PrintHelp(); - exit ( EXIT_FAILURE ); - } + if ( !strcmp(argv[j], "-p") || !strcmp(argv[j],"--pltfile") ) + { + pltfile = argv[++j]; + } + else if ( !strcmp(argv[j], "-g") || !strcmp(argv[j],"--groupfile") ) + { + groupfile = argv[++j]; + } + else if ( !strcmp(argv[j], "-s") || !strcmp(argv[j],"--slicefile") ) + { + slcfile = argv[++j]; + } + else + { + std::cout << "\n\nOption " << argv[j] << " not recognized" << std::endl; + PrintHelp(); + exit ( EXIT_FAILURE ); + } - // Go to the next parameter name - ++j; - } + // Go to the next parameter name + ++j; + } - if (pltfile.empty() || groupfile.empty() || slcfile.empty()) - { - PrintHelp(); - Abort("Missing input file"); - } + if (pltfile.empty() || groupfile.empty() || slcfile.empty()) + { + PrintHelp(); + Abort("Missing input file"); + } - Print() << "\nplotfile = \"" << pltfile << "\"" << std::endl; - Print() << "groupfile = \"" << groupfile << "\"" << std::endl; - Print() << "slicefile = \"" << slcfile << "\"" << std::endl; - Print() << std::endl; + Print() << "\nplotfile = \"" << pltfile << "\"" << std::endl; + Print() << "groupfile = \"" << groupfile << "\"" << std::endl; + Print() << "slicefile = \"" << slcfile << "\"" << std::endl; + Print() << std::endl; - // open the group file and read in the number of groups - std::ifstream group_file; - group_file.open(groupfile); + // open the group file and read in the number of groups + std::ifstream group_file; + group_file.open(groupfile); - string header_line; - string line; + string header_line; + string line; - std::getline(group_file, line); - const std::regex re("=\\s*([0-9]*)"); - std::smatch m; - std::regex_search(line, m, re); + std::getline(group_file, line); + const std::regex re("=\\s*([0-9]*)"); + std::smatch m; + std::regex_search(line, m, re); - int ngroups = stoi(m[1].str()); + int ngroups = stoi(m[1].str()); - group_file.close(); + group_file.close(); - // Start dataservices - DataServices::SetBatchMode(); + // Start dataservices + DataServices::SetBatchMode(); - // Define the type of file - Amrvis::FileType fileType(Amrvis::NEWPLT); - DataServices dataServices (pltfile, fileType); + // Define the type of file + Amrvis::FileType fileType(Amrvis::NEWPLT); + DataServices dataServices (pltfile, fileType); - if (!dataServices.AmrDataOk()) - DataServices::Dispatch(DataServices::ExitRequest, NULL); + if (!dataServices.AmrDataOk()) + DataServices::Dispatch(DataServices::ExitRequest, NULL); - // get data from plot file - AmrData& data = dataServices.AmrDataRef(); + // get data from plot file + AmrData& data = dataServices.AmrDataRef(); - if (data.FinestLevel() != 0) - Abort("ERROR: rhd_shocktube only works for single level"); + if (data.FinestLevel() != 0) + Abort("ERROR: rhd_shocktube only works for single level"); - // get variable names - const Vector& varNames = data.PlotVarNames(); + // get variable names + const Vector& varNames = data.PlotVarNames(); - // get the index bounds and dx. - Box domain = data.ProbDomain()[0]; - Vector dx = data.CellSize(0); - const Vector& problo = data.ProbLo(); - const Vector& probhi = data.ProbHi(); + // get the index bounds and dx. + Box domain = data.ProbDomain()[0]; + Vector dx = data.CellSize(0); + const Vector& problo = data.ProbLo(); + const Vector& probhi = data.ProbHi(); - double dx_coarse = *(std::min_element(dx.begin(), dx.end())); + double dx_coarse = *(std::min_element(dx.begin(), dx.end())); - int nbins = domain.length()[0]; + int nbins = domain.length()[0]; - // x-coord - Vector x(nbins); + // x-coord + Vector x(nbins); - for (auto i = 0; i < nbins; i++) - x[i] = (i + 0.5) * dx_coarse; + for (auto i = 0; i < nbins; i++) + x[i] = (i + 0.5) * dx_coarse; - // find variable indices - Vector compVarNames = {"density", "x_velocity", "pressure", "rad0"}; - auto varComps = GetComponents(data, compVarNames); + // find variable indices + Vector compVarNames = {"density", "x_velocity", "pressure", "rad0"}; + auto varComps = GetComponents(data, compVarNames); - auto dens_comp = varComps[0]; - auto velx_comp = varComps[1]; - auto pres_comp = varComps[2]; - auto rad0_comp = varComps[3]; + auto dens_comp = varComps[0]; + auto velx_comp = varComps[1]; + auto pres_comp = varComps[2]; + auto rad0_comp = varComps[3]; - // allocate storage for data - Vector dens_bin(nbins, 0.); - Vector vel_bin(nbins, 0.); - Vector pres_bin(nbins, 0.); - Vector rad_bin(nbins, 0.); + // allocate storage for data + Vector dens_bin(nbins, 0.); + Vector vel_bin(nbins, 0.); + Vector pres_bin(nbins, 0.); + Vector rad_bin(nbins, 0.); - Vector fill_comps(data.NComp()); - for (auto i = 0; i < data.NComp(); i++) - fill_comps[i] = i; + Vector fill_comps(data.NComp()); + for (auto i = 0; i < data.NComp(); i++) + fill_comps[i] = i; - // extract the 1d data - auto l = 0; + // extract the 1d data + auto l = 0; - Vector level_dx = data.DxLevel()[l]; + Vector level_dx = data.DxLevel()[l]; - const BoxArray& ba = data.boxArray(l); - const DistributionMapping& dm = data.DistributionMap(l); + const BoxArray& ba = data.boxArray(l); + const DistributionMapping& dm = data.DistributionMap(l); - MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); - data.FillVar(lev_data_mf, l, varNames, fill_comps); + MultiFab lev_data_mf(ba, dm, data.NComp(), data.NGrow()); + data.FillVar(lev_data_mf, l, varNames, fill_comps); #ifdef _OPENMP #pragma omp parallel #endif - for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); - - frhdshocktube(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), - nbins, dens_bin.dataPtr(), vel_bin.dataPtr(), - pres_bin.dataPtr(), rad_bin.dataPtr(), - dens_comp, velx_comp, pres_comp, rad0_comp, - ngroups); - } - - // write slicefile - Vector > vars = {dens_bin, vel_bin, pres_bin, rad_bin}; - Vector slcvarNames = {"density", "velocity", "pressure", "rad"}; - WriteSlicefile(nbins, x, slcvarNames, vars, slcfile); - - // destroy timer for profiling - BL_PROFILE_VAR_STOP(pmain); - - amrex::Finalize(); + for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + + frhdshocktube(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), + nbins, dens_bin.dataPtr(), vel_bin.dataPtr(), + pres_bin.dataPtr(), rad_bin.dataPtr(), + dens_comp, velx_comp, pres_comp, rad0_comp, + ngroups); + } + + // write slicefile + Vector > vars = {dens_bin, vel_bin, pres_bin, rad_bin}; + Vector slcvarNames = {"density", "velocity", "pressure", "rad"}; + WriteSlicefile(nbins, x, slcvarNames, vars, slcfile); + + // destroy timer for profiling + BL_PROFILE_VAR_STOP(pmain); + + amrex::Finalize(); } diff --git a/Exec/hydro_tests/riemann_2d/problem_initialize.H b/Exec/hydro_tests/riemann_2d/problem_initialize.H index 80e32b36ba..26811f372b 100644 --- a/Exec/hydro_tests/riemann_2d/problem_initialize.H +++ b/Exec/hydro_tests/riemann_2d/problem_initialize.H @@ -4,15 +4,15 @@ AMREX_INLINE void problem_initialize() { - const Geometry& dgeom = DefaultGeometry(); + const Geometry& dgeom = DefaultGeometry(); - const Real* problo = dgeom.ProbLo(); - const Real* probhi = dgeom.ProbHi(); + const Real* problo = dgeom.ProbLo(); + const Real* probhi = dgeom.ProbHi(); - for (int n=0; n < AMREX_SPACEDIM; ++n) - { - problem::center[n] = 0.5_rt * (problo[n] + probhi[n]); - } + for (int n=0; n < AMREX_SPACEDIM; ++n) + { + problem::center[n] = 0.5_rt * (problo[n] + probhi[n]); + } diff --git a/Exec/hydro_tests/riemann_2d/problem_initialize_state_data.H b/Exec/hydro_tests/riemann_2d/problem_initialize_state_data.H index a1ae1f4e6a..9049bd9bdc 100644 --- a/Exec/hydro_tests/riemann_2d/problem_initialize_state_data.H +++ b/Exec/hydro_tests/riemann_2d/problem_initialize_state_data.H @@ -11,95 +11,95 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE //Declare boxes index i, j, k, the state Array4, and //the geometry information, which emcompases the void problem_initialize_state_data(int i, int j, int k, - Array4 const& state, - const GeometryData& geomdata) + Array4 const& state, + const GeometryData& geomdata) { - //The integer "coord_type" represents the coordinate system - int coord_type = geomdata.Coord(); - - //The pointer *dx represents the real cell size value - const Real* dx = geomdata.CellSize(); + //The integer "coord_type" represents the coordinate system + int coord_type = geomdata.Coord(); + + //The pointer *dx represents the real cell size value + const Real* dx = geomdata.CellSize(); - //The pointer *problo represents the bottom of each dimension - const Real* problo = geomdata.ProbLo(); + //The pointer *problo represents the bottom of each dimension + const Real* problo = geomdata.ProbLo(); #if AMREX_SPACEDIM != 2 - amrex::Abort("The Riemman Problem considered is 2D."); - + amrex::Abort("The Riemman Problem considered is 2D."); + #endif #if AMREX_SPACEDIM == 2 - Real xx = problo[0] + dx[0] * (static_cast(i) + 0.5_rt); - Real yy = problo[1] + dx[1] * (static_cast(j) + 0.5_rt); + Real xx = problo[0] + dx[0] * (static_cast(i) + 0.5_rt); + Real yy = problo[1] + dx[1] * (static_cast(j) + 0.5_rt); #endif - //Now we want to fracture the domain in 4 regions as described in - //Liska & Wendroff (2003). - // - //In order to define the regions we use the following convention - //-------- - // ---|---- - // -2-|--1- - //----+---- - //--3-|--4- - //--------- - //Now, after setting the z-momentum component to zero: - state(i, j, k, UMZ) = 0.0; - - if (xx < problem::center[0] && yy < problem::center[1]) - { - //State in domain 3 - state(i, j, k, URHO) = problem::rho_3; - - state(i, j, k, UMX) = state(i, j, k, URHO) * problem::ux_3; - state(i, j, k, UMY) = state(i, j, k, URHO) * problem::uy_3; - state(i, j, k, UEINT) = problem::p_3 / (eos_gamma - 1.0_rt); - state(i, j, k, UEDEN) = problem::p_3 / (eos_gamma - 1.0_rt) - + 0.5 * problem::rho_3 * problem::ux_3 * problem::ux_3 - + 0.5 * problem::rho_3 * problem::uy_3 * problem::uy_3; - - } else if (xx > problem::center[0] && yy < problem::center[1]) - { - //State in domain 4 - state(i, j, k, URHO) = problem::rho_4; - - state(i, j, k, UMX) = state(i, j, k, URHO) * problem::ux_4; - state(i, j, k, UMY) = state(i, j, k, URHO) * problem::uy_4; - state(i, j, k, UEINT) = problem::p_4 / (eos_gamma - 1.0_rt); - - state(i, j, k, UEDEN) = problem::p_4 / (eos_gamma - 1.0_rt) - + 0.5 * problem::rho_4 * problem::ux_4 * problem::ux_4 - + 0.5 * problem::rho_4 * problem::uy_4 * problem::uy_4; - - } else if (xx < problem::center[0] && yy > problem::center[1]) - { - //State in domain 2 - state(i, j, k, URHO) = problem::rho_2; - - state(i, j, k, UMX) = state(i, j, k, URHO) * problem::ux_2; - state(i, j, k, UMY) = state(i, j, k, URHO) * problem::uy_2; - state(i, j, k, UEINT) = problem::p_2 / (eos_gamma - 1.0_rt); - - state(i, j, k, UEDEN) = problem::p_2 / (eos_gamma - 1.0_rt) - + 0.5 * problem::rho_2 * problem::ux_2 * problem::ux_2 - + 0.5 * problem::rho_2 * problem::uy_2 * problem::uy_2; - } else if (xx > problem::center[0] & yy > problem::center[1]) - { - //State in domain 1 - state(i, j, k, URHO) = problem::rho_1; - - state(i, j, k, UMX) = state(i, j, k, URHO) * problem::ux_1; - state(i, j, k, UMY) = state(i, j, k, URHO) * problem::uy_1; - state(i, j, k, UEINT) = problem::p_1 / (eos_gamma - 1.0_rt); - - state(i, j, k, UEDEN) = problem::p_1 / (eos_gamma - 1.0_rt) - + 0.5 * problem::rho_1 * problem::ux_1 * problem::ux_1 - + 0.5 * problem::rho_1 * problem::uy_1 * problem::uy_1; - } - - //We will consider only one specie. - state(i, j, k,UFS) = state(i,j,k,URHO); + //Now we want to fracture the domain in 4 regions as described in + //Liska & Wendroff (2003). + // + //In order to define the regions we use the following convention + //-------- + // ---|---- + // -2-|--1- + //----+---- + //--3-|--4- + //--------- + //Now, after setting the z-momentum component to zero: + state(i, j, k, UMZ) = 0.0; + + if (xx < problem::center[0] && yy < problem::center[1]) + { + //State in domain 3 + state(i, j, k, URHO) = problem::rho_3; + + state(i, j, k, UMX) = state(i, j, k, URHO) * problem::ux_3; + state(i, j, k, UMY) = state(i, j, k, URHO) * problem::uy_3; + state(i, j, k, UEINT) = problem::p_3 / (eos_gamma - 1.0_rt); + state(i, j, k, UEDEN) = problem::p_3 / (eos_gamma - 1.0_rt) + + 0.5 * problem::rho_3 * problem::ux_3 * problem::ux_3 + + 0.5 * problem::rho_3 * problem::uy_3 * problem::uy_3; + + } else if (xx > problem::center[0] && yy < problem::center[1]) + { + //State in domain 4 + state(i, j, k, URHO) = problem::rho_4; + + state(i, j, k, UMX) = state(i, j, k, URHO) * problem::ux_4; + state(i, j, k, UMY) = state(i, j, k, URHO) * problem::uy_4; + state(i, j, k, UEINT) = problem::p_4 / (eos_gamma - 1.0_rt); + + state(i, j, k, UEDEN) = problem::p_4 / (eos_gamma - 1.0_rt) + + 0.5 * problem::rho_4 * problem::ux_4 * problem::ux_4 + + 0.5 * problem::rho_4 * problem::uy_4 * problem::uy_4; + + } else if (xx < problem::center[0] && yy > problem::center[1]) + { + //State in domain 2 + state(i, j, k, URHO) = problem::rho_2; + + state(i, j, k, UMX) = state(i, j, k, URHO) * problem::ux_2; + state(i, j, k, UMY) = state(i, j, k, URHO) * problem::uy_2; + state(i, j, k, UEINT) = problem::p_2 / (eos_gamma - 1.0_rt); + + state(i, j, k, UEDEN) = problem::p_2 / (eos_gamma - 1.0_rt) + + 0.5 * problem::rho_2 * problem::ux_2 * problem::ux_2 + + 0.5 * problem::rho_2 * problem::uy_2 * problem::uy_2; + } else if (xx > problem::center[0] & yy > problem::center[1]) + { + //State in domain 1 + state(i, j, k, URHO) = problem::rho_1; + + state(i, j, k, UMX) = state(i, j, k, URHO) * problem::ux_1; + state(i, j, k, UMY) = state(i, j, k, URHO) * problem::uy_1; + state(i, j, k, UEINT) = problem::p_1 / (eos_gamma - 1.0_rt); + + state(i, j, k, UEDEN) = problem::p_1 / (eos_gamma - 1.0_rt) + + 0.5 * problem::rho_1 * problem::ux_1 * problem::ux_1 + + 0.5 * problem::rho_1 * problem::uy_1 * problem::uy_1; + } + + //We will consider only one specie. + state(i, j, k,UFS) = state(i,j,k,URHO); } #endif diff --git a/Exec/hydro_tests/test_convect/problem_diagnostics.H b/Exec/hydro_tests/test_convect/problem_diagnostics.H index f75711f189..1285d2dd8c 100644 --- a/Exec/hydro_tests/test_convect/problem_diagnostics.H +++ b/Exec/hydro_tests/test_convect/problem_diagnostics.H @@ -19,22 +19,22 @@ Castro::problem_diagnostics () Castro& ca_lev = getLevel(lev); mass += ca_lev.volWgtSum("density", time); - - rho_E += ca_lev.volWgtSum("rho_E", time); + + rho_E += ca_lev.volWgtSum("rho_E", time); auto temp_mf = ca_lev.derive("Temp", time, 0); - Tmax_level = temp_mf->max(0, 0); - if (Tmax_level > Tmax) - { - Tmax = Tmax_level; - } + Tmax_level = temp_mf->max(0, 0); + if (Tmax_level > Tmax) + { + Tmax = Tmax_level; + } auto mach_mf = ca_lev.derive("MachNumber", time, 0); - MachMax_level = mach_mf->max(0, 0); - if (MachMax_level > MachMax) - { - MachMax = MachMax_level; - } + MachMax_level = mach_mf->max(0, 0); + if (MachMax_level > MachMax) + { + MachMax = MachMax_level; + } } @@ -45,13 +45,13 @@ Castro::problem_diagnostics () std::cout << "TIME= " << time << " MASS = " << mass << '\n'; std::cout << "TIME= " << time << " RHO*E = " << rho_E << '\n'; - // get output files - std::ostream& data_log1 = *Castro::problem_data_logs[0]; + // get output files + std::ostream& data_log1 = *Castro::problem_data_logs[0]; if (time == 0.0) { data_log1 << std::setw(14) << "# time "; data_log1 << std::setw(14) << " max(T) "; - data_log1 << std::setw(14) << " max(M) "; + data_log1 << std::setw(14) << " max(M) "; data_log1 << std::setw(14) << " rho_E "; data_log1 << std::endl; } @@ -63,7 +63,7 @@ Castro::problem_diagnostics () data_log1 << std::setw(14) << std::setprecision(6) << rho_E; data_log1 << std::endl; - std::cout<<'\n'; + std::cout<<'\n'; } } diff --git a/Exec/hydro_tests/toy_convect/problem_diagnostics.H b/Exec/hydro_tests/toy_convect/problem_diagnostics.H index f75711f189..1285d2dd8c 100644 --- a/Exec/hydro_tests/toy_convect/problem_diagnostics.H +++ b/Exec/hydro_tests/toy_convect/problem_diagnostics.H @@ -19,22 +19,22 @@ Castro::problem_diagnostics () Castro& ca_lev = getLevel(lev); mass += ca_lev.volWgtSum("density", time); - - rho_E += ca_lev.volWgtSum("rho_E", time); + + rho_E += ca_lev.volWgtSum("rho_E", time); auto temp_mf = ca_lev.derive("Temp", time, 0); - Tmax_level = temp_mf->max(0, 0); - if (Tmax_level > Tmax) - { - Tmax = Tmax_level; - } + Tmax_level = temp_mf->max(0, 0); + if (Tmax_level > Tmax) + { + Tmax = Tmax_level; + } auto mach_mf = ca_lev.derive("MachNumber", time, 0); - MachMax_level = mach_mf->max(0, 0); - if (MachMax_level > MachMax) - { - MachMax = MachMax_level; - } + MachMax_level = mach_mf->max(0, 0); + if (MachMax_level > MachMax) + { + MachMax = MachMax_level; + } } @@ -45,13 +45,13 @@ Castro::problem_diagnostics () std::cout << "TIME= " << time << " MASS = " << mass << '\n'; std::cout << "TIME= " << time << " RHO*E = " << rho_E << '\n'; - // get output files - std::ostream& data_log1 = *Castro::problem_data_logs[0]; + // get output files + std::ostream& data_log1 = *Castro::problem_data_logs[0]; if (time == 0.0) { data_log1 << std::setw(14) << "# time "; data_log1 << std::setw(14) << " max(T) "; - data_log1 << std::setw(14) << " max(M) "; + data_log1 << std::setw(14) << " max(M) "; data_log1 << std::setw(14) << " rho_E "; data_log1 << std::endl; } @@ -63,7 +63,7 @@ Castro::problem_diagnostics () data_log1 << std::setw(14) << std::setprecision(6) << rho_E; data_log1 << std::endl; - std::cout<<'\n'; + std::cout<<'\n'; } } diff --git a/Exec/radiation_tests/Rad2Tshock/python/LowrieEdwardsUnits.py b/Exec/radiation_tests/Rad2Tshock/python/LowrieEdwardsUnits.py index 544408ee94..758ede71cc 100755 --- a/Exec/radiation_tests/Rad2Tshock/python/LowrieEdwardsUnits.py +++ b/Exec/radiation_tests/Rad2Tshock/python/LowrieEdwardsUnits.py @@ -40,10 +40,10 @@ def LEunits(P0=1.0e-4, gamma=5./3., uL=1.0e5, uT=100.0, mu=1.0, printmessg=0): if __name__ == "__main__": try: - opts, args = getopt.getopt(sys.argv[1:], "L:T:", ["gamma=", "P0=", "mu="]) + opts, args = getopt.getopt(sys.argv[1:], "L:T:", ["gamma=", "P0=", "mu="]) except getopt.GetoptError: - print '???' - sys.exit(1) + print '???' + sys.exit(1) mu = 1.0 P0 = 1.0e-4 @@ -52,7 +52,7 @@ def LEunits(P0=1.0e-4, gamma=5./3., uL=1.0e5, uT=100.0, mu=1.0, printmessg=0): uT = 100.0 for o, a in opts: - if o == "-L": + if o == "-L": uL = float(a) elif o == "-T": uT = float(a) diff --git a/Exec/radiation_tests/Rad2Tshock/python/RadShock.py b/Exec/radiation_tests/Rad2Tshock/python/RadShock.py index 5da2e523a4..d0e02ee277 100755 --- a/Exec/radiation_tests/Rad2Tshock/python/RadShock.py +++ b/Exec/radiation_tests/Rad2Tshock/python/RadShock.py @@ -293,10 +293,10 @@ def dxTdM(xT, M): # Eq. 37 if __name__ == "__main__": try: - opts, args = getopt.getopt(sys.argv[1:], "", ["P0","gamma=","sigma=","kappa=", "M0="]) + opts, args = getopt.getopt(sys.argv[1:], "", ["P0","gamma=","sigma=","kappa=", "M0="]) except getopt.GetoptError: print '???' - sys.exit(1) + sys.exit(1) P0 = 1.0e-4 gamma = 5./3. @@ -305,7 +305,7 @@ def dxTdM(xT, M): # Eq. 37 M0 = 2.0 for o, a in opts: - if o == "--P0": + if o == "--P0": P0 = float(a) elif o == "--gamma": gamma = float(a) diff --git a/Exec/radiation_tests/Rad2Tshock/python/initincgs.py b/Exec/radiation_tests/Rad2Tshock/python/initincgs.py index a8a9dec302..60f6a80236 100755 --- a/Exec/radiation_tests/Rad2Tshock/python/initincgs.py +++ b/Exec/radiation_tests/Rad2Tshock/python/initincgs.py @@ -7,10 +7,10 @@ def main(): try: - opts, args = getopt.getopt(sys.argv[1:], "L:T:", ["P0","gamma=","sigma=","kappa=", "M0=", "mu="]) + opts, args = getopt.getopt(sys.argv[1:], "L:T:", ["P0","gamma=","sigma=","kappa=", "M0=", "mu="]) except getopt.GetoptError: print '???' - sys.exit(1) + sys.exit(1) P0 = 1.0e-4 gamma = 5./3. @@ -21,7 +21,7 @@ def main(): uT = 100.0 mu = 1.0 for o, a in opts: - if o == "--P0": + if o == "--P0": P0 = float(a) elif o == "--gamma": gamma = float(a) diff --git a/Exec/radiation_tests/RadShestakovBolstad/python/initunits.py b/Exec/radiation_tests/RadShestakovBolstad/python/initunits.py index b4e6841228..3570e7779a 100755 --- a/Exec/radiation_tests/RadShestakovBolstad/python/initunits.py +++ b/Exec/radiation_tests/RadShestakovBolstad/python/initunits.py @@ -30,10 +30,10 @@ def SBunits(T0, rho0, kap0, printmessg=0): if __name__ == "__main__": try: - opts, args = getopt.getopt(sys.argv[1:], "", ["rho0=", "kap0=", "T0="]) + opts, args = getopt.getopt(sys.argv[1:], "", ["rho0=", "kap0=", "T0="]) except getopt.GetoptError: - print '???' - sys.exit(1) + print '???' + sys.exit(1) rho0 = 1.8212111e-5 kap0 = 4.0628337e43 diff --git a/Exec/radiation_tests/RadSuOlsonMG/python/initunits.py b/Exec/radiation_tests/RadSuOlsonMG/python/initunits.py index d894593803..3f55f39c43 100755 --- a/Exec/radiation_tests/RadSuOlsonMG/python/initunits.py +++ b/Exec/radiation_tests/RadSuOlsonMG/python/initunits.py @@ -32,10 +32,10 @@ def SOunits(eps, printmessg=0): if __name__ == "__main__": try: - opts, args = getopt.getopt(sys.argv[1:], "", ["eps="]) + opts, args = getopt.getopt(sys.argv[1:], "", ["eps="]) except getopt.GetoptError: - print '???' - sys.exit(1) + print '???' + sys.exit(1) eps = 1.0 diff --git a/Exec/science/Detonation/problem_initialize_state_data.H b/Exec/science/Detonation/problem_initialize_state_data.H index 1372340168..3b08ce86a9 100644 --- a/Exec/science/Detonation/problem_initialize_state_data.H +++ b/Exec/science/Detonation/problem_initialize_state_data.H @@ -66,7 +66,7 @@ void problem_initialize_state_data (int i, int j, int k, auto nse_state = get_actual_nse_state(burn_state); for (int n = 0; n < NumSpec; ++n) { - state(i,j,k,UFS+n) = state(i,j,k,URHO) * nse_state.xn[n]; + state(i,j,k,UFS+n) = state(i,j,k,URHO) * nse_state.xn[n]; } state(i,j,k,UMUP) = nse_state.mu_p; diff --git a/Exec/science/massive_star/problem_initialize_state_data.H b/Exec/science/massive_star/problem_initialize_state_data.H index 68fd2acfd8..4276fa879a 100644 --- a/Exec/science/massive_star/problem_initialize_state_data.H +++ b/Exec/science/massive_star/problem_initialize_state_data.H @@ -130,32 +130,32 @@ void problem_initialize_state_data (int i, int j, int k, if (problem::interpolate_pres == 1) { - // we need to get T from P, rho, but consistent with NSE (if - // we are in NSE) - if (nse_check) { - // this will also change the composition, since the new T - // alters the NSE - nse_T_abar_from_p(eos_state.rho, eos_state.p, eos_state.aux[AuxZero::iye], - eos_state.T, eos_state.aux[AuxZero::iabar]); - state(i,j,k,UTEMP) = eos_state.T; - state(i,j,k,UFX+AuxZero::iabar) = eos_state.aux[AuxZero::iabar]; - - // now call the EOS with the new T and abar to get e - eos(eos_input_rt, eos_state); - - // finally, get the updated B/A - nse_table_t nse_state; - nse_state.T = state(i,j,k,UTEMP); - nse_state.rho = state(i,j,k,URHO); - nse_state.Ye = state(i,j,k,UFX+AuxZero::iye); - nse_interp(nse_state); - - state(i,j,k,UFX+AuxZero::ibea) = nse_state.bea; - - } else { - eos(eos_input_rp, eos_state); - state(i,j,k,UTEMP) = eos_state.T; - } + // we need to get T from P, rho, but consistent with NSE (if + // we are in NSE) + if (nse_check) { + // this will also change the composition, since the new T + // alters the NSE + nse_T_abar_from_p(eos_state.rho, eos_state.p, eos_state.aux[AuxZero::iye], + eos_state.T, eos_state.aux[AuxZero::iabar]); + state(i,j,k,UTEMP) = eos_state.T; + state(i,j,k,UFX+AuxZero::iabar) = eos_state.aux[AuxZero::iabar]; + + // now call the EOS with the new T and abar to get e + eos(eos_input_rt, eos_state); + + // finally, get the updated B/A + nse_table_t nse_state; + nse_state.T = state(i,j,k,UTEMP); + nse_state.rho = state(i,j,k,URHO); + nse_state.Ye = state(i,j,k,UFX+AuxZero::iye); + nse_interp(nse_state); + + state(i,j,k,UFX+AuxZero::ibea) = nse_state.bea; + + } else { + eos(eos_input_rp, eos_state); + state(i,j,k,UTEMP) = eos_state.T; + } } else { eos(eos_input_rt, eos_state); } diff --git a/Exec/unit_tests/particles_test/Prob.cpp b/Exec/unit_tests/particles_test/Prob.cpp index 04f18ae213..73c2d08ac5 100644 --- a/Exec/unit_tests/particles_test/Prob.cpp +++ b/Exec/unit_tests/particles_test/Prob.cpp @@ -10,111 +10,111 @@ using ParticleType = Particle<1, 1>; #ifdef DO_PROBLEM_POST_SIMULATION void Castro::problem_post_simulation(Vector >& amr_level) { - int nlevels = amr_level.size(); + int nlevels = amr_level.size(); - Real err = -1.e30; + Real err = -1.e30; - auto lev = 0; - auto ipass = 0; + auto lev = 0; + auto ipass = 0; - Castro& castro = dynamic_cast(*amr_level[lev]); + Castro& castro = dynamic_cast(*amr_level[lev]); - const Real strttime = amrex::second(); - const Real* dx = castro.geom.CellSize(); - const Real* plo = castro.geom.ProbLo(); - const Real* phi = castro.geom.ProbHi(); + const Real strttime = amrex::second(); + const Real* dx = castro.geom.CellSize(); + const Real* plo = castro.geom.ProbLo(); + const Real* phi = castro.geom.ProbHi(); - // read initial positions from Ascii file + // read initial positions from Ascii file - const int MyProc = ParallelDescriptor::MyProc(); + const int MyProc = ParallelDescriptor::MyProc(); - Vector nparticles; + Vector nparticles; - std::string file = "particle_file"; + std::string file = "particle_file"; - VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); + VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); - std::ifstream ifs; + std::ifstream ifs; - ifs.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); + ifs.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); - ifs.open(file.c_str(), std::ios::in); + ifs.open(file.c_str(), std::ios::in); - if (!ifs.good()) - amrex::FileOpenFailed(file); + if (!ifs.good()) + amrex::FileOpenFailed(file); - int num_particles = 0; + int num_particles = 0; - ifs >> num_particles >> std::ws; + ifs >> num_particles >> std::ws; - ParticleLocData pld; - ParticleType p, p_rep; + ParticleLocData pld; + ParticleType p, p_rep; - for (int i = 0; i < num_particles; i++) - { - AMREX_D_TERM(ifs >> p.m_pos[0]; , - ifs >> p.m_pos[1]; , - ifs >> p.m_pos[2]; ); + for (int i = 0; i < num_particles; i++) + { + AMREX_D_TERM(ifs >> p.m_pos[0]; , + ifs >> p.m_pos[1]; , + ifs >> p.m_pos[2]; ); - auto extradata = 0; + auto extradata = 0; - for (int n = 0; n < extradata; n++) - { - ifs >> p.m_rdata[AMREX_SPACEDIM+n]; - } + for (int n = 0; n < extradata; n++) + { + ifs >> p.m_rdata[AMREX_SPACEDIM+n]; + } - p.id() = ParticleType::NextID(); - p.cpu() = MyProc; + p.id() = ParticleType::NextID(); + p.cpu() = MyProc; - nparticles.push_back(p); - } + nparticles.push_back(p); + } - // compute the change in position of the particles wrt to their initial positions + // compute the change in position of the particles wrt to their initial positions - auto& pmap = TracerPC->GetParticles(lev); + auto& pmap = TracerPC->GetParticles(lev); double total_err = 0; - for (auto& kv : pmap) { - int grid = kv.first.first; - auto& pbox = kv.second.GetArrayOfStructs(); - const int n = pbox.size(); + for (auto& kv : pmap) { + int grid = kv.first.first; + auto& pbox = kv.second.GetArrayOfStructs(); + const int n = pbox.size(); - for (int i = 0; i < n; i++) - { - auto& p = pbox[i]; + for (int i = 0; i < n; i++) + { + auto& p = pbox[i]; - bool match = false; + bool match = false; - auto it = nparticles.begin(); + auto it = nparticles.begin(); - // find the original particle and calculate change in position - for (; !match && it != nparticles.end(); ++it) { - if (it->id() == p.id()) - match = true; - } + // find the original particle and calculate change in position + for (; !match && it != nparticles.end(); ++it) { + if (it->id() == p.id()) + match = true; + } - if (!match) Print() << "haven't found a match :(" << std::endl; + if (!match) Print() << "haven't found a match :(" << std::endl; - // calculate change in position - auto deltax = it->m_pos[0] - p.m_pos[0]; - auto deltay = it->m_pos[1] - p.m_pos[1]; + // calculate change in position + auto deltax = it->m_pos[0] - p.m_pos[0]; + auto deltay = it->m_pos[1] - p.m_pos[1]; - auto delta = sqrt(deltax*deltax + deltay*deltay); + auto delta = sqrt(deltax*deltax + deltay*deltay); // quick nan check here if (delta == delta) total_err += delta; - // Print() << "delta r = " << delta << std::endl; + // Print() << "delta r = " << delta << std::endl; - } - } + } + } - const std::string stars(78,'*'); - amrex::Print() << stars << "\n" - << " particles problem post_simulation() \n" - << " Average change in position: " << total_err / num_particles << "\n" - << stars << "\n\n"; + const std::string stars(78,'*'); + amrex::Print() << stars << "\n" + << " particles problem post_simulation() \n" + << " Average change in position: " << total_err / num_particles << "\n" + << stars << "\n\n"; } #endif diff --git a/Source/driver/sum_utils.cpp b/Source/driver/sum_utils.cpp index 433c01393f..f0fa10e35e 100644 --- a/Source/driver/sum_utils.cpp +++ b/Source/driver/sum_utils.cpp @@ -314,10 +314,10 @@ Castro::locSquaredSum (const std::string& name, #ifdef GRAVITY void Castro::gwstrain (Real time, - Real& h_plus_1, Real& h_cross_1, - Real& h_plus_2, Real& h_cross_2, - Real& h_plus_3, Real& h_cross_3, - bool local) { + Real& h_plus_1, Real& h_cross_1, + Real& h_plus_2, Real& h_cross_2, + Real& h_plus_3, Real& h_cross_3, + bool local) { amrex::ignore_unused(time); diff --git a/Source/driver/timestep.cpp b/Source/driver/timestep.cpp index 01af85ce3c..42ac0325b2 100644 --- a/Source/driver/timestep.cpp +++ b/Source/driver/timestep.cpp @@ -438,8 +438,8 @@ Castro::estdt_burning (int is_new) #endif #ifdef NSE_NET - burn_state.mu_p = S(i,j,k,UMUP); - burn_state.mu_n = S(i,j,k,UMUN); + burn_state.mu_p = S(i,j,k,UMUP); + burn_state.mu_n = S(i,j,k,UMUN); #endif if (!in_nse(burn_state)) { diff --git a/Util/ConvertCheckpoint/Embiggen.cpp b/Util/ConvertCheckpoint/Embiggen.cpp index d5a6136eb9..f8c2d43b97 100644 --- a/Util/ConvertCheckpoint/Embiggen.cpp +++ b/Util/ConvertCheckpoint/Embiggen.cpp @@ -403,7 +403,7 @@ static void ReadCheckpointFile(const std::string& fileName) { FullPathName = fileName; if( ! fileName.empty() && fileName[fileName.length()-1] != '/') { FullPathName += '/'; - } + } FullPathName += mf_name; VisMF::Read(*(falRef.state[i].old_data), FullPathName); } @@ -458,7 +458,7 @@ static void ReadCheckpointFile(const std::string& fileName) { int ncomp = falRef_orig.state[i].new_data->nComp(); int ngrow = falRef_orig.state[i].new_data->nGrow(); - DistributionMapping dmap {falRef.grids}; + DistributionMapping dmap {falRef.grids}; falRef.state[i].new_data = new MultiFab(falRef.grids, dmap, ncomp, ngrow); falRef.state[i].new_data->setVal(0.); @@ -502,17 +502,17 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string if (oldCastroHeaderFile.good()) { - std::ofstream newCastroHeaderFile; + std::ofstream newCastroHeaderFile; - std::string newCastroHeaderName = outFileName + "/CastroHeader"; - newCastroHeaderFile.open(newCastroHeaderName.c_str(), std::ios::binary); + std::string newCastroHeaderName = outFileName + "/CastroHeader"; + newCastroHeaderFile.open(newCastroHeaderName.c_str(), std::ios::binary); - if (newCastroHeaderFile.good()) { - newCastroHeaderFile << oldCastroHeaderFile.rdbuf(); - newCastroHeaderFile.close(); - } + if (newCastroHeaderFile.good()) { + newCastroHeaderFile << oldCastroHeaderFile.rdbuf(); + newCastroHeaderFile.close(); + } - oldCastroHeaderFile.close(); + oldCastroHeaderFile.close(); } @@ -533,15 +533,15 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string if(ParallelDescriptor::IOProcessor()) { // Only the IOProcessor() writes to the header file. HeaderFile.open(HeaderFileName.c_str(), - std::ios::out|std::ios::trunc|std::ios::binary); + std::ios::out|std::ios::trunc|std::ios::binary); if( ! HeaderFile.good()) { amrex::FileOpenFailed(HeaderFileName); - } + } old_prec = HeaderFile.precision(15); - int max_level(fakeAmr.finest_level); + int max_level(fakeAmr.finest_level); HeaderFile << CheckPointVersion << '\n' << AMREX_SPACEDIM << '\n' << fakeAmr.cumtime << '\n' @@ -627,9 +627,9 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string if(ParallelDescriptor::IOProcessor()) { // The relative name gets written to the Header file. std::string mf_name_old = name; - mf_name_old += OldSuffix; + mf_name_old += OldSuffix; std::string mf_name_new = name; - mf_name_new += NewSuffix; + mf_name_new += NewSuffix; os << falRef.state[i].domain << '\n'; @@ -663,7 +663,7 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string BL_ASSERT(dump_old); BL_ASSERT(falRef.state[i].old_data); std::string mf_fullpath_old = fullpathname; - mf_fullpath_old += OldSuffix; + mf_fullpath_old += OldSuffix; VisMF::Write(*(falRef.state[i].old_data),mf_fullpath_old,how); } // ++++++++++++ @@ -689,7 +689,7 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string if( ! HeaderFile.good()) { amrex::Error("Amr::checkpoint() failed"); - } + } } FArrayBox::setFormat(thePrevFormat); diff --git a/Util/code_checker/tab_exterminator.sh b/Util/code_checker/tab_exterminator.sh index 80ec329600..61f954a234 100755 --- a/Util/code_checker/tab_exterminator.sh +++ b/Util/code_checker/tab_exterminator.sh @@ -1,6 +1,6 @@ for f in $(grep -rIPl "\t" ../../Source) do - echo "Converting tabs to spaces in $f" - expand -t 8 $f > tmp_file - mv tmp_file $f + echo "Converting tabs to spaces in $f" + expand -t 8 $f > tmp_file + mv tmp_file $f done diff --git a/Util/model_parser/model_parser.H b/Util/model_parser/model_parser.H index fa7dee3f1f..db9f4e522b 100644 --- a/Util/model_parser/model_parser.H +++ b/Util/model_parser/model_parser.H @@ -44,22 +44,22 @@ namespace model_string { inline std::string& ltrim(std::string& s) { - auto it = std::find_if(s.begin(), s.end(), + auto it = std::find_if(s.begin(), s.end(), [](int c) { return !std::isspace(c); }); - s.erase(s.begin(), it); - return s; + s.erase(s.begin(), it); + return s; } inline std::string& rtrim(std::string& s) { - auto it = std::find_if(s.rbegin(), s.rend(), + auto it = std::find_if(s.rbegin(), s.rend(), [](int c) { return !std::isspace(c); }); - s.erase(it.base(), s.end()); - return s; + s.erase(it.base(), s.end()); + return s; } } @@ -125,11 +125,11 @@ interpolate(const Real r, const int var_index, const int model_index=0) { // safety check to make sure interp lies within the bounding points. We don't // do this at the lower boundary, which usually corresponds to the center of the star. if (r >= model::profile(model_index).r(0)) { - Real minvar = std::min(model::profile(model_index).state(id+1, var_index), - model::profile(model_index).state(id, var_index)); - Real maxvar = std::max(model::profile(model_index).state(id+1, var_index), - model::profile(model_index).state(id, var_index)); - interp = std::clamp(interp, minvar, maxvar); + Real minvar = std::min(model::profile(model_index).state(id+1, var_index), + model::profile(model_index).state(id, var_index)); + Real maxvar = std::max(model::profile(model_index).state(id+1, var_index), + model::profile(model_index).state(id, var_index)); + interp = std::clamp(interp, minvar, maxvar); } return interp;