diff --git a/CMakeLists.txt b/CMakeLists.txt index d112397..b75a4dc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.22) project( euler2d DESCRIPTION "euler2d_kokkos is mini application using C++/Kokkos library, used for teaching purpose" - LANGUAGES CXX C) + LANGUAGES CXX C Fortran) # # default local cmake macro repository diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bb7136a..f329a08 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -17,6 +17,8 @@ add_executable(euler2d SimpleTimer.h SimpleTimer.cpp CudaTimer.h + cnpy/cnpy.h + cnpy/cnpy.cpp main.cpp) target_include_directories(euler2d @@ -34,3 +36,5 @@ configure_file(test_implode.ini test_implode.ini COPYONLY) configure_file(test_implode_big.ini test_implode_big.ini COPYONLY) configure_file(test_four_quadrant.ini test_four_quadrant.ini COPYONLY) configure_file(test_discontinuity.ini test_discontinuity.ini COPYONLY) + +add_subdirectory(sedov_blast) diff --git a/src/ComputeRadialProfileFunctor.h b/src/ComputeRadialProfileFunctor.h new file mode 100644 index 0000000..65d1bf2 --- /dev/null +++ b/src/ComputeRadialProfileFunctor.h @@ -0,0 +1,169 @@ +/** + * \file ComputeRadialProfile.h + */ +#ifndef EULER2D_COMPUTE_RADIAL_PROFILE_FUNCTOR_H_ +#define EULER2D_COMPUTE_RADIAL_PROFILE_FUNCTOR_H_ + +#include "kokkos_shared.h" +#include "HydroBaseFunctor.h" +#include "HydroParams.h" +#include "real_type.h" + +// other utilities +#include "cnpy/cnpy_io.h" + +namespace euler2d +{ + +/*************************************************/ +/*************************************************/ +/*************************************************/ +/** + * Compute radial profile by averaging all direction. + * + * This functor is mainly aimed at checking numerical and analytical solution of the radial Sedov + * blast. + * + * \tparam device_t is the Kokkos device use for computation (CPU, GPU, ...) + */ +template +class ComputeRadialProfileFunctor : HydroBaseFunctor +{ + +public: + using DataArray_t = DataArray; + + //! our kokkos execution space + using exec_space = typename device_t::execution_space; + +private: + //! heavy data - conservative variables + DataArray_t m_Udata; + + //! number of bins in radial direction + const int m_nbins; + + //! max radial distance + const real_t m_max_radial_distance; + + //! center of the box + Kokkos::Array m_box_center; + + //! radial distance histogram + Kokkos::View> m_radial_distance_histo; + + //! averaged density radial profile + Kokkos::View> m_density_profile; + +public: + ComputeRadialProfileFunctor(HydroParams params, + DataArray_t Udata, + real_t max_radial_distance, + Kokkos::Array box_center) + : HydroBaseFunctor(params) + , m_Udata(Udata) + , m_nbins(params.blast_nbins) + , m_max_radial_distance(max_radial_distance) + , m_box_center(box_center) + , m_radial_distance_histo("radial_distance_histo", m_nbins) + , m_density_profile("density_profile", m_nbins) + { + Kokkos::deep_copy(m_radial_distance_histo, 0); + Kokkos::deep_copy(m_density_profile, ZERO_F); + }; + + // ==================================================================== + // ==================================================================== + //! static method which does it all: create and execute functor using range policy + //! + static void + apply(HydroParams params, DataArray_t Udata) + { + + const real_t xmin = params.xmin; + const real_t ymin = params.ymin; + const real_t xmax = params.xmax; + const real_t ymax = params.ymax; + const real_t dx = params.dx; + const real_t dy = params.dy; + + // get center of the box coordinates + Kokkos::Array box_center; + + box_center[IX] = (xmin + xmax) / 2; + box_center[IY] = (ymin + ymax) / 2; + + // compute maximum radial distance + const auto Dx = (xmax - xmin) / 2; + const auto Dy = (ymax - ymin) / 2; + const auto max_radial_distance = sqrt(Dx * Dx + Dy * Dy); + + ComputeRadialProfileFunctor functor(params, Udata, max_radial_distance, box_center); + + Kokkos::parallel_for( + "euler2d::ComputeRadialProfileFunctor", + Kokkos::MDRangePolicy>({ 0, 0 }, { params.isize, params.jsize }), + functor); + + // once the radial profile is computed, gather all the MPI pieces + const auto nbins = params.blast_nbins; + auto total_radial_distance_histo = + Kokkos::View("total_radial_distance_histo", nbins); + + const auto total_density_profile = + Kokkos::View("total_density_profile", nbins); + + // then we compute the profile and output in numpy format + auto rad_dist = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, functor.m_radial_distance_histo); + auto dens_prof = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, functor.m_density_profile); + + auto distances = Kokkos::View("distances", nbins); + const auto dr = max_radial_distance / nbins; + + for (int i = 0; i < nbins; ++i) + { + distances(i) = (i + 0.5) * dr; + dens_prof(i) /= rad_dist(i); + } + + // now we can output results + save_cnpy(distances, "sedov_blast_radial_distances"); + save_cnpy(dens_prof, "sedov_blast_density_profile"); + + } // apply + + // ==================================================================== + // ==================================================================== + KOKKOS_INLINE_FUNCTION + void + operator()(const int & i, const int & j) const + { + + const int ghostWidth = params.ghostWidth; + const real_t xmin = params.xmin; + const real_t ymin = params.ymin; + const real_t dx = params.dx; + const real_t dy = params.dy; + const real_t x = xmin + dx / 2 + (i - ghostWidth) * dx; + const real_t y = ymin + dy / 2 + (j - ghostWidth) * dy; + + const auto distance = sqrt((x - m_box_center[IX]) * (x - m_box_center[IX]) + + (y - m_box_center[IY]) * (y - m_box_center[IY])); + + // compute bin number from distance to center of the box + const auto bin = (int)(distance / m_max_radial_distance * m_nbins); + + // update density profile + m_radial_distance_histo(bin) += 1; + + m_density_profile(bin) += m_Udata(i, j, ID); + + } // operator() + +}; // ComputeRadialProfileFunctor + +} // namespace euler2d + +#endif // EULER2D_COMPUTE_RADIAL_PROFILE_FUNCTOR_H_ diff --git a/src/HydroParams.cpp b/src/HydroParams.cpp index 2208af7..5a9f38e 100644 --- a/src/HydroParams.cpp +++ b/src/HydroParams.cpp @@ -112,6 +112,8 @@ HydroParams::setup(ConfigMap & configMap) blast_density_out = configMap.getFloat("blast", "density_out", 1.2); blast_pressure_in = configMap.getFloat("blast", "pressure_in", 10.0); blast_pressure_out = configMap.getFloat("blast", "pressure_out", 0.1); + blast_total_energy_inside = configMap.getFloat("blast", "total_energy_inside", 0.0); + blast_nbins = configMap.getInteger("blast", "nbins", 100); implementationVersion = configMap.getFloat("OTHER", "implementationVersion", 0); if (implementationVersion != 0 and implementationVersion != 1 and implementationVersion != 2) diff --git a/src/HydroParams.h b/src/HydroParams.h index 1bb6580..10bff0f 100644 --- a/src/HydroParams.h +++ b/src/HydroParams.h @@ -183,6 +183,8 @@ struct HydroParams real_t blast_density_out; real_t blast_pressure_in; real_t blast_pressure_out; + real_t blast_total_energy_inside; + int blast_nbins; // other parameters int implementationVersion = 0; /*!< triggers which implementation to use (currently 3 versions)*/ @@ -224,6 +226,7 @@ struct HydroParams , blast_density_out(1.2) , blast_pressure_in(10.0) , blast_pressure_out(0.1) + , blast_total_energy_inside(0) , implementationVersion(0) {} diff --git a/src/HydroRun.h b/src/HydroRun.h index 7e77a3e..393c31c 100644 --- a/src/HydroRun.h +++ b/src/HydroRun.h @@ -118,7 +118,7 @@ class HydroRun * Paraview/Visit to read these h5 files as a time series. */ void - write_xdmf_time_series(); + write_xdmf_time_series(int last_time_step); #endif }; // class HydroRun @@ -696,7 +696,8 @@ HydroRun::godunov_unsplit_impl(DataArray_t data_in, * point to data file : xdmf2d.h5 */ template - void HydroRun::write_xdmf_time_series() + void + HydroRun::write_xdmf_time_series(int last_time_step) { FILE * xdmf = 0; const int & nx = params.nx; @@ -722,8 +723,15 @@ HydroRun::godunov_unsplit_impl(DataArray_t data_in, " \n"); // for each time step write a item - for (int iStep = 0; iStep <= params.nStepmax; iStep += params.nOutput) + const int nbOutput = (params.nStepmax + params.nOutput - 1) / params.nOutput + 1; + for (int iOut = 0; iOut < nbOutput; ++iOut) { + int iStep = iOut * params.nOutput; + + if (iOut == (nbOutput - 1)) + { + iStep = last_time_step; + } std::ostringstream outNum; outNum.width(7); @@ -786,7 +794,7 @@ HydroRun::godunov_unsplit_impl(DataArray_t data_in, fprintf(xdmf, "\n"); fclose(xdmf); - } // HydroRun::write_xdmf_xml + } // HydroRun::write_xdmf_xml #endif // USE_HDF5 } // namespace euler2d diff --git a/src/HydroRunFunctors.h b/src/HydroRunFunctors.h index b8e78e4..22bfe0f 100644 --- a/src/HydroRunFunctors.h +++ b/src/HydroRunFunctors.h @@ -1384,32 +1384,63 @@ class InitImplodeFunctor : public HydroBaseFunctor /*************************************************/ /*************************************************/ /*************************************************/ -template +/** + * Initialize blast test case. + * + * Two types of initialization: + * + * - if parameter total_energy_inside is positive, we initialize in total energy + * - if parameter total_energy_inside is negative or null, we initialize using pressure in / + * pressure out + * + * If you want to do the well know sedov blast test, you need to initialize using total energy, not + * pressure. + */ +template class InitBlastFunctor : public HydroBaseFunctor { public: + struct TagRegularInit + {}; + struct TagCorrectTotalEnergy + {}; + using DataArray_t = DataArray; using exec_space = typename device_t::execution_space; InitBlastFunctor(HydroParams params, DataArray_t Udata) : HydroBaseFunctor(params) - , Udata(Udata){}; + , Udata(Udata) {}; // static method which does it all: create and execute functor static void apply(HydroParams params, DataArray_t Udata) { + real_t volume_inside = ZERO_F; + Kokkos::Sum reducer(volume_inside); + InitBlastFunctor functor(params, Udata); - Kokkos::parallel_for( - "InitBlast", - Kokkos::MDRangePolicy>({ 0, 0 }, { params.isize, params.jsize }), - functor); + Kokkos::parallel_reduce("InitBlast - regular", + Kokkos::MDRangePolicy, TagRegularInit>( + { 0, 0 }, { params.isize, params.jsize }), + functor, + reducer); + + if (params.blast_total_energy_inside > 0) + { + functor.m_volume_inside = volume_inside; + Kokkos::parallel_for( + "InitBlast - correct total energy", + Kokkos::MDRangePolicy, TagCorrectTotalEnergy>( + { 0, 0 }, { params.isize, params.jsize }), + functor); + } } KOKKOS_INLINE_FUNCTION void - operator()(const int & i, const int & j) const + operator()(TagRegularInit const &, const int & i, const int & j, real_t & volume) const { const int ghostWidth = params.ghostWidth; @@ -1443,6 +1474,8 @@ class InitBlastFunctor : public HydroBaseFunctor Udata(i, j, IP) = blast_pressure_in / (gamma0 - 1.0); Udata(i, j, IU) = 0.0; Udata(i, j, IV) = 0.0; + + volume += dx * dy; } else { @@ -1452,9 +1485,44 @@ class InitBlastFunctor : public HydroBaseFunctor Udata(i, j, IV) = 0.0; } - } // end operator () + } // end operator () - TagRegularInit + + KOKKOS_INLINE_FUNCTION + void + operator()(TagCorrectTotalEnergy const &, const int & i, const int & j) const + { + + const int ghostWidth = params.ghostWidth; + + const real_t xmin = params.xmin; + const real_t ymin = params.ymin; + const real_t dx = params.dx; + const real_t dy = params.dy; + + const real_t gamma0 = params.settings.gamma0; + + // blast problem parameters + const real_t blast_radius = params.blast_radius; + const real_t radius2 = blast_radius * blast_radius; + const real_t blast_center_x = params.blast_center_x; + const real_t blast_center_y = params.blast_center_y; + const real_t blast_total_energy_inside = params.blast_total_energy_inside; + + real_t x = xmin + dx / 2 + (i - ghostWidth) * dx; + real_t y = ymin + dy / 2 + (j - ghostWidth) * dy; + + real_t d2 = + (x - blast_center_x) * (x - blast_center_x) + (y - blast_center_y) * (y - blast_center_y); + + if (d2 < radius2) + { + Udata(i, j, IP) = blast_total_energy_inside / m_volume_inside; + } + + } // end operator () - TagCorrectTotalEnergy DataArray_t Udata; + real_t m_volume_inside = 0.0; }; // InitBlastFunctor diff --git a/src/cnpy/cnpy.cpp b/src/cnpy/cnpy.cpp new file mode 100644 index 0000000..9bae156 --- /dev/null +++ b/src/cnpy/cnpy.cpp @@ -0,0 +1,341 @@ +//Copyright (C) 2011 Carl Rogers +//Released under MIT License +//license available in LICENSE file, or at http://www.opensource.org/licenses/mit-license.php + +#include"cnpy.h" +#include +#include +#include +#include +#include +#include +#include +#include + +char cnpy::BigEndianTest() { + int x = 1; + return (((char *)&x)[0]) ? '<' : '>'; +} + +char cnpy::map_type(const std::type_info& t) +{ + if(t == typeid(float) ) return 'f'; + if(t == typeid(double) ) return 'f'; + if(t == typeid(long double) ) return 'f'; + + if(t == typeid(int) ) return 'i'; + if(t == typeid(char) ) return 'i'; + if(t == typeid(short) ) return 'i'; + if(t == typeid(long) ) return 'i'; + if(t == typeid(long long) ) return 'i'; + + if(t == typeid(unsigned char) ) return 'u'; + if(t == typeid(unsigned short) ) return 'u'; + if(t == typeid(unsigned long) ) return 'u'; + if(t == typeid(unsigned long long) ) return 'u'; + if(t == typeid(unsigned int) ) return 'u'; + + if(t == typeid(bool) ) return 'b'; + + if(t == typeid(std::complex) ) return 'c'; + if(t == typeid(std::complex) ) return 'c'; + if(t == typeid(std::complex) ) return 'c'; + + else return '?'; +} + +template<> std::vector& cnpy::operator+=(std::vector& lhs, const std::string rhs) { + lhs.insert(lhs.end(),rhs.begin(),rhs.end()); + return lhs; +} + +template<> std::vector& cnpy::operator+=(std::vector& lhs, const char* rhs) { + //write in little endian + size_t len = strlen(rhs); + lhs.reserve(len); + for(size_t byte = 0; byte < len; byte++) { + lhs.push_back(rhs[byte]); + } + return lhs; +} + +void cnpy::parse_npy_header(unsigned char* buffer,size_t& word_size, std::vector& shape, bool& fortran_order) { + //std::string magic_string(buffer,6); + uint8_t major_version = *reinterpret_cast(buffer+6); + uint8_t minor_version = *reinterpret_cast(buffer+7); + uint16_t header_len = *reinterpret_cast(buffer+8); + std::string header(reinterpret_cast(buffer+9),header_len); + + size_t loc1, loc2; + + //fortran order + loc1 = header.find("fortran_order")+16; + fortran_order = (header.substr(loc1,4) == "True" ? true : false); + + //shape + loc1 = header.find("("); + loc2 = header.find(")"); + + std::regex num_regex("[0-9][0-9]*"); + std::smatch sm; + shape.clear(); + + std::string str_shape = header.substr(loc1+1,loc2-loc1-1); + while(std::regex_search(str_shape, sm, num_regex)) { + shape.push_back(std::stoi(sm[0].str())); + str_shape = sm.suffix().str(); + } + + //endian, word size, data type + //byte order code | stands for not applicable. + //not sure when this applies except for byte array + loc1 = header.find("descr")+9; + bool littleEndian = (header[loc1] == '<' || header[loc1] == '|' ? true : false); + assert(littleEndian); + + //char type = header[loc1+1]; + //assert(type == map_type(T)); + + std::string str_ws = header.substr(loc1+2); + loc2 = str_ws.find("'"); + word_size = atoi(str_ws.substr(0,loc2).c_str()); +} + +void cnpy::parse_npy_header(FILE* fp, size_t& word_size, std::vector& shape, bool& fortran_order) { + char buffer[256]; + size_t res = fread(buffer,sizeof(char),11,fp); + if(res != 11) + throw std::runtime_error("parse_npy_header: failed fread"); + std::string header = fgets(buffer,256,fp); + assert(header[header.size()-1] == '\n'); + + size_t loc1, loc2; + + //fortran order + loc1 = header.find("fortran_order"); + if (loc1 == std::string::npos) + throw std::runtime_error("parse_npy_header: failed to find header keyword: 'fortran_order'"); + loc1 += 16; + fortran_order = (header.substr(loc1,4) == "True" ? true : false); + + //shape + loc1 = header.find("("); + loc2 = header.find(")"); + if (loc1 == std::string::npos || loc2 == std::string::npos) + throw std::runtime_error("parse_npy_header: failed to find header keyword: '(' or ')'"); + + std::regex num_regex("[0-9][0-9]*"); + std::smatch sm; + shape.clear(); + + std::string str_shape = header.substr(loc1+1,loc2-loc1-1); + while(std::regex_search(str_shape, sm, num_regex)) { + shape.push_back(std::stoi(sm[0].str())); + str_shape = sm.suffix().str(); + } + + //endian, word size, data type + //byte order code | stands for not applicable. + //not sure when this applies except for byte array + loc1 = header.find("descr"); + if (loc1 == std::string::npos) + throw std::runtime_error("parse_npy_header: failed to find header keyword: 'descr'"); + loc1 += 9; + bool littleEndian = (header[loc1] == '<' || header[loc1] == '|' ? true : false); + assert(littleEndian); + + //char type = header[loc1+1]; + //assert(type == map_type(T)); + + std::string str_ws = header.substr(loc1+2); + loc2 = str_ws.find("'"); + word_size = atoi(str_ws.substr(0,loc2).c_str()); +} + +void cnpy::parse_zip_footer(FILE* fp, uint16_t& nrecs, size_t& global_header_size, size_t& global_header_offset) +{ + std::vector footer(22); + fseek(fp,-22,SEEK_END); + size_t res = fread(&footer[0],sizeof(char),22,fp); + if(res != 22) + throw std::runtime_error("parse_zip_footer: failed fread"); + + uint16_t disk_no, disk_start, nrecs_on_disk, comment_len; + disk_no = *(uint16_t*) &footer[4]; + disk_start = *(uint16_t*) &footer[6]; + nrecs_on_disk = *(uint16_t*) &footer[8]; + nrecs = *(uint16_t*) &footer[10]; + global_header_size = *(uint32_t*) &footer[12]; + global_header_offset = *(uint32_t*) &footer[16]; + comment_len = *(uint16_t*) &footer[20]; + + assert(disk_no == 0); + assert(disk_start == 0); + assert(nrecs_on_disk == nrecs); + assert(comment_len == 0); +} + +cnpy::NpyArray load_the_npy_file(FILE* fp) +{ + std::vector shape; + size_t word_size; + bool fortran_order; + cnpy::parse_npy_header(fp,word_size,shape,fortran_order); + + cnpy::NpyArray arr(shape, word_size, fortran_order); + size_t nread = fread(arr.data(),1,arr.num_bytes(),fp); + if(nread != arr.num_bytes()) + throw std::runtime_error("load_the_npy_file: failed fread"); + return arr; +} + +cnpy::NpyArray load_the_npz_array(FILE* fp, uint32_t compr_bytes, uint32_t uncompr_bytes) { + + std::vector buffer_compr(compr_bytes); + std::vector buffer_uncompr(uncompr_bytes); + size_t nread = fread(&buffer_compr[0],1,compr_bytes,fp); + if(nread != compr_bytes) + throw std::runtime_error("load_the_npy_file: failed fread"); + + int err; + z_stream d_stream; + + d_stream.zalloc = Z_NULL; + d_stream.zfree = Z_NULL; + d_stream.opaque = Z_NULL; + d_stream.avail_in = 0; + d_stream.next_in = Z_NULL; + err = inflateInit2(&d_stream, -MAX_WBITS); + + d_stream.avail_in = compr_bytes; + d_stream.next_in = &buffer_compr[0]; + d_stream.avail_out = uncompr_bytes; + d_stream.next_out = &buffer_uncompr[0]; + + err = inflate(&d_stream, Z_FINISH); + err = inflateEnd(&d_stream); + + std::vector shape; + size_t word_size; + bool fortran_order; + cnpy::parse_npy_header(&buffer_uncompr[0],word_size,shape,fortran_order); + + cnpy::NpyArray array(shape, word_size, fortran_order); + + size_t offset = uncompr_bytes - array.num_bytes(); + memcpy(array.data(),&buffer_uncompr[0]+offset,array.num_bytes()); + + return array; +} + +cnpy::npz_t cnpy::npz_load(std::string fname) { + FILE* fp = fopen(fname.c_str(),"rb"); + + if(!fp) { + throw std::runtime_error("npz_load: Error! Unable to open file "+fname+"!"); + } + + cnpy::npz_t arrays; + + while(1) { + std::vector local_header(30); + size_t headerres = fread(&local_header[0],sizeof(char),30,fp); + if(headerres != 30) + throw std::runtime_error("npz_load: failed fread"); + + //if we've reached the global header, stop reading + if(local_header[2] != 0x03 || local_header[3] != 0x04) break; + + //read in the variable name + uint16_t name_len = *(uint16_t*) &local_header[26]; + std::string varname(name_len,' '); + size_t vname_res = fread(&varname[0],sizeof(char),name_len,fp); + if(vname_res != name_len) + throw std::runtime_error("npz_load: failed fread"); + + //erase the lagging .npy + varname.erase(varname.end()-4,varname.end()); + + //read in the extra field + uint16_t extra_field_len = *(uint16_t*) &local_header[28]; + if(extra_field_len > 0) { + std::vector buff(extra_field_len); + size_t efield_res = fread(&buff[0],sizeof(char),extra_field_len,fp); + if(efield_res != extra_field_len) + throw std::runtime_error("npz_load: failed fread"); + } + + uint16_t compr_method = *reinterpret_cast(&local_header[0]+8); + uint32_t compr_bytes = *reinterpret_cast(&local_header[0]+18); + uint32_t uncompr_bytes = *reinterpret_cast(&local_header[0]+22); + + if(compr_method == 0) {arrays[varname] = load_the_npy_file(fp);} + else {arrays[varname] = load_the_npz_array(fp,compr_bytes,uncompr_bytes);} + } + + fclose(fp); + return arrays; +} + +cnpy::NpyArray cnpy::npz_load(std::string fname, std::string varname) { + FILE* fp = fopen(fname.c_str(),"rb"); + + if(!fp) throw std::runtime_error("npz_load: Unable to open file "+fname); + + while(1) { + std::vector local_header(30); + size_t header_res = fread(&local_header[0],sizeof(char),30,fp); + if(header_res != 30) + throw std::runtime_error("npz_load: failed fread"); + + //if we've reached the global header, stop reading + if(local_header[2] != 0x03 || local_header[3] != 0x04) break; + + //read in the variable name + uint16_t name_len = *(uint16_t*) &local_header[26]; + std::string vname(name_len,' '); + size_t vname_res = fread(&vname[0],sizeof(char),name_len,fp); + if(vname_res != name_len) + throw std::runtime_error("npz_load: failed fread"); + vname.erase(vname.end()-4,vname.end()); //erase the lagging .npy + + //read in the extra field + uint16_t extra_field_len = *(uint16_t*) &local_header[28]; + fseek(fp,extra_field_len,SEEK_CUR); //skip past the extra field + + uint16_t compr_method = *reinterpret_cast(&local_header[0]+8); + uint32_t compr_bytes = *reinterpret_cast(&local_header[0]+18); + uint32_t uncompr_bytes = *reinterpret_cast(&local_header[0]+22); + + if(vname == varname) { + NpyArray array = (compr_method == 0) ? load_the_npy_file(fp) : load_the_npz_array(fp,compr_bytes,uncompr_bytes); + fclose(fp); + return array; + } + else { + //skip past the data + uint32_t size = *(uint32_t*) &local_header[22]; + fseek(fp,size,SEEK_CUR); + } + } + + fclose(fp); + + //if we get here, we haven't found the variable in the file + throw std::runtime_error("npz_load: Variable name "+varname+" not found in "+fname); +} + +cnpy::NpyArray cnpy::npy_load(std::string fname) { + + FILE* fp = fopen(fname.c_str(), "rb"); + + if(!fp) throw std::runtime_error("npy_load: Unable to open file "+fname); + + NpyArray arr = load_the_npy_file(fp); + + fclose(fp); + return arr; +} + + + diff --git a/src/cnpy/cnpy.h b/src/cnpy/cnpy.h new file mode 100644 index 0000000..18d0486 --- /dev/null +++ b/src/cnpy/cnpy.h @@ -0,0 +1,319 @@ +// Copyright (C) 2011 Carl Rogers +// Released under MIT License +// license available in LICENSE file, or at +// http://www.opensource.org/licenses/mit-license.php + +#ifndef LIBCNPY_H_ +#define LIBCNPY_H_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace cnpy { + +struct NpyArray { + NpyArray(const std::vector &_shape, size_t _word_size, + bool _fortran_order) + : shape(_shape), word_size(_word_size), fortran_order(_fortran_order) { + num_vals = 1; + for (size_t i = 0; i < shape.size(); i++) + num_vals *= shape[i]; + data_holder = std::shared_ptr>( + new std::vector(num_vals * word_size)); + } + + NpyArray() : shape(0), word_size(0), fortran_order(0), num_vals(0) {} + + template T *data() { + return reinterpret_cast(&(*data_holder)[0]); + } + + template const T *data() const { + return reinterpret_cast(&(*data_holder)[0]); + } + + template std::vector as_vec() const { + const T *p = data(); + return std::vector(p, p + num_vals); + } + + size_t num_bytes() const { return data_holder->size(); } + + std::shared_ptr> data_holder; + std::vector shape; + size_t word_size; + bool fortran_order; + size_t num_vals; +}; + +using npz_t = std::map; + +char BigEndianTest(); +char map_type(const std::type_info &t); +template +std::vector create_npy_header(const std::vector &shape); +void parse_npy_header(FILE *fp, size_t &word_size, std::vector &shape, + bool &fortran_order); +void parse_npy_header(unsigned char *buffer, size_t &word_size, + std::vector &shape, bool &fortran_order); +void parse_zip_footer(FILE *fp, uint16_t &nrecs, size_t &global_header_size, + size_t &global_header_offset); +npz_t npz_load(std::string fname); +NpyArray npz_load(std::string fname, std::string varname); +NpyArray npy_load(std::string fname); + +template +std::vector &operator+=(std::vector &lhs, const T rhs) { + // write in little endian + for (size_t byte = 0; byte < sizeof(T); byte++) { + char val = *((char *)&rhs + byte); + lhs.push_back(val); + } + return lhs; +} + +template <> +std::vector &operator+=(std::vector &lhs, const std::string rhs); +template <> +std::vector &operator+=(std::vector &lhs, const char *rhs); + +template +void npy_save(std::string fname, const T *data, const std::vector shape, + std::string mode = "w") { + FILE *fp = NULL; + std::vector + true_data_shape; // if appending, the shape of existing + new data + + if (mode == "a") + fp = fopen(fname.c_str(), "r+b"); + + if (fp) { + // file exists. we need to append to it. read the header, modify the array + // size + size_t word_size; + bool fortran_order; + parse_npy_header(fp, word_size, true_data_shape, fortran_order); + assert(!fortran_order); + + if (word_size != sizeof(T)) { + std::cout << "libnpy error: " << fname << " has word size " << word_size + << " but npy_save appending data sized " << sizeof(T) << "\n"; + assert(word_size == sizeof(T)); + } + if (true_data_shape.size() != shape.size()) { + std::cout << "libnpy error: npy_save attempting to append misdimensioned " + "data to " + << fname << "\n"; + assert(true_data_shape.size() != shape.size()); + } + + for (size_t i = 1; i < shape.size(); i++) { + if (shape[i] != true_data_shape[i]) { + std::cout + << "libnpy error: npy_save attempting to append misshaped data to " + << fname << "\n"; + assert(shape[i] == true_data_shape[i]); + } + } + true_data_shape[0] += shape[0]; + } else { + fp = fopen(fname.c_str(), "wb"); + true_data_shape = shape; + } + + std::vector header = create_npy_header(true_data_shape); + size_t nels = + std::accumulate(shape.begin(), shape.end(), 1, std::multiplies()); + + fseek(fp, 0, SEEK_SET); + fwrite(&header[0], sizeof(char), header.size(), fp); + fseek(fp, 0, SEEK_END); + fwrite(data, sizeof(T), nels, fp); + fclose(fp); +} + +template +void npy_save(std::string fname, const T *data, const unsigned int *shape, + const unsigned int ndims, std::string mode = "w") { + std::vector shapev; + for (unsigned int i = 0; i < ndims; ++i) + shapev.push_back(shape[i]); + + npy_save(fname, data, shapev, mode); +} + +template +void npz_save(std::string zipname, std::string fname, const T *data, + const std::vector &shape, std::string mode = "w") { + // first, append a .npy to the fname + fname += ".npy"; + + // now, on with the show + FILE *fp = NULL; + uint16_t nrecs = 0; + size_t global_header_offset = 0; + std::vector global_header; + + if (mode == "a") + fp = fopen(zipname.c_str(), "r+b"); + + if (fp) { + // zip file exists. we need to add a new npy file to it. + // first read the footer. this gives us the offset and size of the global + // header then read and store the global header. below, we will write the the + // new data at the start of the global header then append the global header + // and footer below it + size_t global_header_size; + parse_zip_footer(fp, nrecs, global_header_size, global_header_offset); + fseek(fp, global_header_offset, SEEK_SET); + global_header.resize(global_header_size); + size_t res = fread(&global_header[0], sizeof(char), global_header_size, fp); + if (res != global_header_size) { + throw std::runtime_error( + "npz_save: header read error while adding to existing zip"); + } + fseek(fp, global_header_offset, SEEK_SET); + } else { + fp = fopen(zipname.c_str(), "wb"); + } + + std::vector npy_header = create_npy_header(shape); + + size_t nels = + std::accumulate(shape.begin(), shape.end(), 1, std::multiplies()); + size_t nbytes = nels * sizeof(T) + npy_header.size(); + + // get the CRC of the data to be added + uint32_t crc = crc32(0L, (uint8_t *)&npy_header[0], npy_header.size()); + crc = crc32(crc, (uint8_t *)data, nels * sizeof(T)); + + // build the local header + std::vector local_header; + local_header += "PK"; // first part of sig + local_header += (uint16_t)0x0403; // second part of sig + local_header += (uint16_t)20; // min version to extract + local_header += (uint16_t)0; // general purpose bit flag + local_header += (uint16_t)0; // compression method + local_header += (uint16_t)0; // file last mod time + local_header += (uint16_t)0; // file last mod date + local_header += (uint32_t)crc; // crc + local_header += (uint32_t)nbytes; // compressed size + local_header += (uint32_t)nbytes; // uncompressed size + local_header += (uint16_t)fname.size(); // fname length + local_header += (uint16_t)0; // extra field length + local_header += fname; + + // build global header + global_header += "PK"; // first part of sig + global_header += (uint16_t)0x0201; // second part of sig + global_header += (uint16_t)20; // version made by + global_header.insert(global_header.end(), local_header.begin() + 4, + local_header.begin() + 30); + global_header += (uint16_t)0; // file comment length + global_header += (uint16_t)0; // disk number where file starts + global_header += (uint16_t)0; // internal file attributes + global_header += (uint32_t)0; // external file attributes + global_header += (uint32_t) + global_header_offset; // relative offset of local file header, since it + // begins where the global header used to begin + global_header += fname; + + // build footer + std::vector footer; + footer += "PK"; // first part of sig + footer += (uint16_t)0x0605; // second part of sig + footer += (uint16_t)0; // number of this disk + footer += (uint16_t)0; // disk where footer starts + footer += (uint16_t)(nrecs + 1); // number of records on this disk + footer += (uint16_t)(nrecs + 1); // total number of records + footer += (uint32_t)global_header.size(); // nbytes of global headers + footer += (uint32_t)( + global_header_offset + nbytes + + local_header.size()); // offset of start of global headers, since global + // header now starts after newly written array + footer += (uint16_t)0; // zip file comment length + + // write everything + fwrite(&local_header[0], sizeof(char), local_header.size(), fp); + fwrite(&npy_header[0], sizeof(char), npy_header.size(), fp); + fwrite(data, sizeof(T), nels, fp); + fwrite(&global_header[0], sizeof(char), global_header.size(), fp); + fwrite(&footer[0], sizeof(char), footer.size(), fp); + fclose(fp); +} + +template +void npz_save(std::string zipname, std::string fname, const T *data, + const unsigned int *shape, const unsigned int ndims, + std::string mode = "w") { + std::vector shapev; + for (unsigned int i = 0; i < ndims; ++i) + shapev.push_back(shape[i]); + + npz_save(zipname, fname, data, shapev, mode); +} + +template +void npy_save(std::string fname, const std::vector data, + std::string mode = "w") { + std::vector shape; + shape.push_back(data.size()); + npy_save(fname, &data[0], shape, mode); +} + +template +void npz_save(std::string zipname, std::string fname, const std::vector data, + std::string mode = "w") { + std::vector shape; + shape.push_back(data.size()); + npz_save(zipname, fname, &data[0], shape, mode); +} + +template +std::vector create_npy_header(const std::vector &shape) { + + std::vector dict; + dict += "{'descr': '"; + dict += BigEndianTest(); + dict += map_type(typeid(T)); + dict += std::to_string(sizeof(T)); + dict += "', 'fortran_order': False, 'shape': ("; + dict += std::to_string(shape[0]); + for (size_t i = 1; i < shape.size(); i++) { + dict += ", "; + dict += std::to_string(shape[i]); + } + if (shape.size() == 1) + dict += ","; + dict += "), }"; + // pad with spaces so that preamble+dict is modulo 16 bytes. preamble is 10 + // bytes. dict needs to end with \n + int remainder = 16 - (10 + dict.size()) % 16; + dict.insert(dict.end(), remainder, ' '); + dict.back() = '\n'; + + std::vector header; + header += (char)0x93; + header += "NUMPY"; + header += (char)0x01; // major version of numpy format + header += (char)0x00; // minor version of numpy format + header += (uint16_t)dict.size(); + header.insert(header.end(), dict.begin(), dict.end()); + + return header; +} + +} // namespace cnpy + +#endif diff --git a/src/cnpy/cnpy_io.h b/src/cnpy/cnpy_io.h new file mode 100644 index 0000000..8301229 --- /dev/null +++ b/src/cnpy/cnpy_io.h @@ -0,0 +1,63 @@ +/** + * \file cnpy_io.h + * \brief converting numpy files to/from a Kokkos::View + */ +#ifndef EULER2D_CNPY_IO_H_ +#define EULER2D_CNPY_IO_H_ + +#include "../kokkos_shared.h" + +#include +#include +#include +#include // for std::setfill + +#include "./cnpy.h" + +namespace euler2d +{ + +// ============================================================================= +// ============================================================================= +/** + * Save a multidimensional Kokkos::View into a file using numpy format. + * + * example of use in python: + * data = np.load('data.npy') + * + * \tparam View is a Kokkos::View class, we check that the view is accessible from host + * + * \param[in] v a View instance containing a data array to be saved + */ +template +void +save_cnpy(View v, std::string name = "") +{ + constexpr bool view_accessible_from_host = + Kokkos::SpaceAccessibility::accessible; + + if constexpr (view_accessible_from_host) + { + + // save in numpy format + std::vector shape; + if constexpr (View::rank >= 1) + shape.push_back(v.extent(0)); + if constexpr (View::rank >= 2) + shape.push_back(v.extent(1)); + if constexpr (View::rank >= 3) + shape.push_back(v.extent(2)); + + std::ostringstream ss; + if (name.size() > 0) + ss << name << ".npy"; + else + ss << v.label() << ".npy"; + cnpy::npy_save(ss.str(), v.data(), shape, "w"); + } +} // save_cnpy + +} // namespace euler2d + +#endif // EULER2D_CNPY_IO_H_ diff --git a/src/main.cpp b/src/main.cpp index 791fb48..30c0fad 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -11,6 +11,7 @@ #include "kokkos_shared.h" #include "HydroBaseFunctor.h" +#include "ComputeRadialProfileFunctor.h" #include "HydroParams.h" // read parameter file #include "HydroRun.h" // memory allocation for hydro arrays #include "real_type.h" // choose between single and double precision @@ -20,7 +21,7 @@ int main(int argc, char * argv[]) { using DefaultDevice = - Kokkos::Device; + Kokkos::Device; using device = DefaultDevice; // using device = Kokkos::OpenMP::device_type; // using device = Kokkos::CudaSpace::device_type; @@ -121,6 +122,13 @@ main(int argc, char * argv[]) // compute new dt dt_timer.start(); dt = hydro->compute_dt(nStep % 2); + + // correct dt if necessary + if (t + dt > params.tEnd) + { + dt = params.tEnd - t; + } + dt_timer.stop(); // perform one step integration @@ -135,7 +143,7 @@ main(int argc, char * argv[]) // save last time step if (params.enableOutput) { - if (params.nOutput > 0 and nStep % params.nOutput == 0) + if (params.nOutput > 0) { std::cout << "Output results at time t=" << t << " step " << nStep << " dt=" << dt << std::endl; @@ -146,7 +154,7 @@ main(int argc, char * argv[]) hydro->saveData(hydro->U2, nStep, "U"); io_timer.stop(); } // end output - } // end enable output + } // end enable output // end of computation total_timer.stop(); @@ -154,10 +162,18 @@ main(int argc, char * argv[]) // write XDMF wrapper #ifdef USE_HDF5 + printf("Last time step is : %d\n", nStep); + if (params.nOutput > 0 and params.ioHDF5) - hydro->write_xdmf_time_series(); + hydro->write_xdmf_time_series(nStep); #endif + // post-processing for Sedov blast + if (params.problemType == euler2d::PROBLEM_BLAST and params.blast_total_energy_inside > 0) + { + euler2d::ComputeRadialProfileFunctor::apply(params, hydro->U); + } + // print monitoring information { int isize = params.isize; diff --git a/src/sedov_blast/CMakeLists.txt b/src/sedov_blast/CMakeLists.txt new file mode 100644 index 0000000..4a6ad2b --- /dev/null +++ b/src/sedov_blast/CMakeLists.txt @@ -0,0 +1,4 @@ +add_subdirectory(sedov_fortran) + +configure_file(test_sedov_blast_2d.ini test_sedov_blast_2d.ini COPYONLY) +configure_file(sedov_plot.py sedov_plot.py COPYONLY) diff --git a/src/sedov_blast/README.md b/src/sedov_blast/README.md new file mode 100644 index 0000000..2706e1f --- /dev/null +++ b/src/sedov_blast/README.md @@ -0,0 +1,15 @@ +# Comparing numerical and analytical solution of the Sedov blast run + +1. run kanop to get the radial profile of density: +```shell +../euler2d ./test_sedov_blast_2d.ini +``` +this will create 2 output files: +- `sedov_blast_radial_distances.npy` : radial distance data +- `sedov_blast_density_profile.npy` : density data + +2. run python script `sedov_plot.py`: +```shell +./sedov_plot.py --ini ./test_sedov_blast_2d.ini +``` +this will compute the analytical solution (using fortran code from F.X. Timmes) and plot the results. diff --git a/src/sedov_blast/sedov_fortran/CMakeLists.txt b/src/sedov_blast/sedov_fortran/CMakeLists.txt new file mode 100644 index 0000000..df28a2a --- /dev/null +++ b/src/sedov_blast/sedov_fortran/CMakeLists.txt @@ -0,0 +1,5 @@ +add_executable(sedov3_qp + sedov3_qp.f90) + +# TODO: refactor sedov3_qp.f90 to avoid all the warnings +remove_definitions(-Wall) diff --git a/src/sedov_blast/sedov_fortran/README b/src/sedov_blast/sedov_fortran/README new file mode 100644 index 0000000..f70305b --- /dev/null +++ b/src/sedov_blast/sedov_fortran/README @@ -0,0 +1,10 @@ +compile, with say + +% ifort -o sedov3_qp.exe sedov3_qp.f90 + +run the executable + +% ./sedov3_qp.exe + standard + xgeom = 3.000000E+00 gamma= 1.400000E+00 omega = 0.000000E+00 alpha = 8.510719E-01 j1 = 2.962690E-02 j2 = 2.116508E-02 + r2 = 1.000000E+00 rho2 = 6.000000E+00 u2 = 3.333333E-01 e2 = 5.555556E-02 p2 = 1.333333E-01 cs2 = 1.763834E-01 diff --git a/src/sedov_blast/sedov_fortran/README_kanop.md b/src/sedov_blast/sedov_fortran/README_kanop.md new file mode 100644 index 0000000..c7370f6 --- /dev/null +++ b/src/sedov_blast/sedov_fortran/README_kanop.md @@ -0,0 +1,5 @@ +# Sedov blast analytical solution + +the fortran code here is a slightly modified version of [sedov.tbz](https://cococubed.com/codes/sedov/sedov.tbz) developped by F.X. Timmes + +We just modified the user interface so we can configure at which time the solution is computed. diff --git a/src/sedov_blast/sedov_fortran/sedov3_qp.dek b/src/sedov_blast/sedov_fortran/sedov3_qp.dek new file mode 100644 index 0000000..d57bada --- /dev/null +++ b/src/sedov_blast/sedov_fortran/sedov3_qp.dek @@ -0,0 +1,9 @@ +! common block communication + logical :: ltransitional,lstandard,lvacuum,lomega2,lomega3 + real*16 :: gamma,gamm1,gamp1,gpogm,xgeom,xg2,rwant,r2, & + a0,a1,a2,a3,a4,a5,a_val,b_val,c_val,d_val,e_val, & + omega,vv,xlam_want,vwant,rvv + common /slap/ gamma,gamm1,gamp1,gpogm,xgeom,xg2,rwant,r2, & + a0,a1,a2,a3,a4,a5,a_val,b_val,c_val,d_val,e_val, & + omega,vv,xlam_want,vwant,rvv, & + ltransitional,lstandard,lvacuum,lomega2,lomega3 diff --git a/src/sedov_blast/sedov_fortran/sedov3_qp.f90 b/src/sedov_blast/sedov_fortran/sedov3_qp.f90 new file mode 100644 index 0000000..7845df7 --- /dev/null +++ b/src/sedov_blast/sedov_fortran/sedov3_qp.f90 @@ -0,0 +1,172 @@ +program sedov3 + implicit none + + ! exercises the sedov solver + + ! declare + character*80 :: outfile,string + integer :: i,nstep,iargc + integer, parameter :: nmax = 1000 + real*16 :: time,zpos(nmax), & + eblast,rho0,omega,vel0,ener0,pres0,cs0,gamma, & + xgeom,rshock,rho2,u2,e2,p2,cs2,rvacuum, & + den(nmax),ener(nmax),pres(nmax),vel(nmax), & + cs(nmax), & + zlo,zhi,zstep,value + + + ! popular formats +01 format(1x,t5,a,t8,a,t32,a,t56,a,t80,a,t104,a,t128,a,t132,a) +02 format(1x,i4,1p8e12.4) +03 format(1x,i4,1p8e24.16) + + + + ! if your compiler/linker handles command line arguments + ! get the number of spatial points, blast energy, geometry type, + ! density exponent, and output file name + + i = iargc() + + if (i .eq. 0) then + nstep = 120 + eblast = 0.851072q0 + xgeom = 3.0q0 + omega = 0.0q0 + gamma = 1.4q0 + time = 0.1q0 + outfile = 'spherical_standard_omega0p00.dat' + + else if (i .eq. 7) then + + call getarg(1,string) + nstep = int(value(string)) + + call getarg(2,string) + eblast = value(string) + + call getarg(3,string) + xgeom = value(string) + + call getarg(4,string) + omega = value(string) + + call getarg(5,string) + gamma = value(string) + + call getarg(6,string) + time = value(string) + + call getarg(7,outfile) + + else + stop 'pass in 7 parameters: nstep eblast xgeom omega gamma tend outfile' + end if + + ! input parameters in cgs + rho0 = 1.0q0 + vel0 = 0.0q0 + ener0 = 0.0q0 + pres0 = 0.0q0 + cs0 = 0.0q0 + + ! or explicitly set stuff + ! standard cases + ! spherical constant density should reach r=1 at t=1 + ! nstep = 120 + ! eblast = 0.851072q0 + ! xgeom = 3.0q0 + ! omega = 0.0q0 + ! outfile = 'spherical_standard_omega0p00.dat' + + ! cylindrical constant density should reach r=1 at t=1 + ! nstep = 120 + ! eblast = 0.311357q0 + ! xgeom = 2.0q0 + ! omega = 0.0q0 + ! outfile = 'cylindrical_standard_omega0p00.dat' + + ! planar constant density should reach x=1 at t=1 + ! nstep = 120 + ! eblast = 0.0673185q0 + ! xgeom = 1.0q0 + ! omega = 0.0q0 + ! outfile = 'planar_standard_omega0p00.dat' + + + ! singular cases + ! spherical should reach r=1 at t=1 + ! nstep = 120 + ! eblast = 4.90875q0 + ! xgeom = 3.0q0 + ! omega = 7.0q0/3.0q0 + ! outfile = 'spherical_singular_omega2p33.dat' + + ! cylindrical should reach r=0.75 at t=1 + ! nstep = 120 + ! eblast = 2.45749q0 + ! xgeom = 2.0q0 + ! omega = 5.0q0/3.0q0 + ! outfile = 'cylindrical_singular_omega1p66.dat' + + + ! vacuum cases + ! spherical should reach r=1 at t=1 + ! nstep = 120 + ! eblast = 5.45670q0 + ! xgeom = 3.0q0 + ! omega = 2.4q0 + ! outfile = 'spherical_vacuum_omega2p40.dat' + + ! cylindrical should reach r=0.75 at t=1 + ! nstep = 120 + ! eblast = 2.67315q0 + ! xgeom = 2.0q0 + ! omega = 1.7q0 + ! outfile = 'cylindrical_vacuum_omega1p70.dat' + + + + ! number of grid points, spatial domain, spatial step size + ! to match rage output, use the mid-cell points + zlo = 0.0q0 + zhi = 1.2q0 + zstep = (zhi - zlo)/float(nstep-1) + do i=1,nstep + ! zpos(i) = zlo + 0.5q0*zstep + float(i-1)*zstep + zpos(i) = zlo + float(i-1)*zstep + enddo + + + ! get the solution for all spatial points at once + + call sed_1d(time,nstep,zpos, & + eblast,omega,xgeom, & + rho0,vel0,ener0,pres0,cs0,gamma, & + rshock,rho2,u2,e2,p2,cs2,rvacuum, & + den,ener,pres,vel,cs) + + + + ! output file + open(unit=2,file=outfile,status='unknown') + write(2,02) nstep,time + write(2,01) 'i','x','den','ener','pres','vel','cs' + do i=1,nstep + + write(2,03) i,zpos(i),den(i),ener(i),pres(i),vel(i),cs(i) + + ! write(2,03) i,zpos(i),den(i)/rho2,ener(i)/e2, & + ! pres(i)/p2,vel(i)/u2,cs(i)/cs2 + + enddo + close(unit=2) + + write(6,*) + + ! close up shop +end program sedov3 + + + +include 'sedov_library_qp.f90' diff --git a/src/sedov_blast/sedov_fortran/sedov_library_qp.f90 b/src/sedov_blast/sedov_fortran/sedov_library_qp.f90 new file mode 100644 index 0000000..c555c5e --- /dev/null +++ b/src/sedov_blast/sedov_fortran/sedov_library_qp.f90 @@ -0,0 +1,1221 @@ + subroutine sed_1d(time,nstep,xpos, & + eblast,omega_in,xgeom_in, & + rho0,vel0,ener0,pres0,cs0,gam0, & + rshock,rho2,u2,e2,p2,cs2,rvacuum, & + den,ener,pres,vel,cs) + implicit none + + +! this routine produces 1d solutions for a sedov blast wave propagating +! through a density gradient rho = rho**(-omega) +! in planar, cylindrical or spherical geometry +! for the standard, transitional and vaccum cases: + +! standard : a nonzero solution extends from the shock to the origin, where the pressure is finite. +! transitional : a nonzero solution extends from the shock to the origin, where the pressure vanishes. +! vacuum : a nonzero solution extends from the shock to a boundary point, where the density vanishes. + + +! input: +! time = temporal point where solution is desired seconds +! xpos(i) = spatial points where solution is desired cm +! eblast = energy of blast erg +! rho0 = ambient density g/cm**3 rho = rho0 * r**(-omega_in) +! omegain = density power law exponent rho = rho0 * r**(-omega_in) +! vel0 = ambient material speed cm/s +! pres0 = ambient pressure erg/cm**3 +! cs0 = ambient sound speed cm/s +! gam0 = gamma law equation of state +! xgeom_in = geometry factor, 3=spherical, 2=cylindircal, 1=planar + + +! for efficiency (i.e., doing the energy integrals only once), +! this routine returns the solution for an array of spatial points +! at the desired time point. + +! output: +! den(i) = density g/cm**3 +! ener(i) = specific internal energy erg/g +! pres(i) = presssure erg/cm**3 +! vel(i) = velocity cm/s +! cs(i) = sound speed cm/s + + +! this routine is based upon two papers: +! "evaluation of the sedov-von neumann-taylor blast wave solution" +! jim kamm, la-ur-00-6055 +! "the sedov self-similiar point blast solutions in nonuniform media" +! david book, shock waves, 4, 1, 1994 + +! although the ordinary differential equations are analytic, +! the sedov expressions appear to become singular for various +! combinations of parameters and at the lower limits of the integration +! range. all these singularies are removable and done so by this routine. + +! real*16 because the real*8 implementations run out of precision, for example, +! "near" the origin in the standard case or the transition region in the vacuum case. + + +! declare the pass + integer :: nstep + real*16 :: time,xpos(*), & + eblast,rho0,omega_in,vel0,ener0,pres0,cs0, & + gam0,xgeom_in,rshock,rho2,u2,e2,p2,cs2,rvacuum, & + den(*),ener(*),pres(*),vel(*),cs(*) + +! local variables + external :: midpnt,midpowl,midpowl2,sed_v_find,sed_r_find, & + efun01,efun02 + integer :: i + real*16 :: efun01,efun02,eval1,eval2 + real*16 :: v0,v2,vstar,vmin,alpha,us, & + zeroin,sed_v_find,sed_r_find, & + vat,l_fun,dlamdv,f_fun,g_fun,h_fun, & + denom2,denom3,rho1 + + +! eps controls the integration accuracy, don't get too greedy or the number +! of function evaluations required kills. +! eps2 controls the root find accuracy +! osmall controls the size of transition regions + + integer, parameter :: iprint = 1 + real*16, parameter :: eps = 1.0q-10, & + eps2 = 1.0q-28, & + osmall = 1.0q-4, & + pi = 3.1415926535897932384626433832795029q0 + + +! common block communication + include 'sedov3_qp.dek' + + +! common block communication with the integration stepper + real*16 :: gam_int + common /cmidp/ gam_int + + +! popular formats + 87 format(1x,1p10e14.6) + 88 format(1x,8(a7,1pe14.6,' ')) + + +! initialize the solution + do i=1,nstep + den(i) = 0.0q0 + vel(i) = 0.0q0 + pres(i) = 0.0q0 + ener(i) = 0.0q0 + cs(i) = 0.0q0 + end do + + +! return on unphysical cases +! infinite mass + if (omega_in .ge. xgeom_in) return + + + +! transfer the pass to common block and create some frequent combinations + gamma = gam0 + gamm1 = gamma - 1.0q0 + gamp1 = gamma + 1.0q0 + gpogm = gamp1 / gamm1 + xgeom = xgeom_in + omega = omega_in + xg2 = xgeom + 2.0q0 - omega + denom2 = 2.0q0*gamm1 + xgeom - gamma*omega + denom3 = xgeom * (2.0q0 - gamma) - omega + + +! post shock location v2 and location of transitional point vstar +! kamm equation 18 and 19 + + v2 = 4.0q0 / (xg2 * gamp1) + vstar = 2.0q0 / (gamm1*xgeom + 2.0q0) + + +! set logicals that determines the type of solution + + lstandard = .false. + ltransitional = .false. + lvacuum = .false. + + if (abs(v2 - vstar) .le. osmall) then + ltransitional = .true. + if (iprint .eq. 1) write(6,*) 'transitional' + else if (v2 .lt. vstar - osmall) then + lstandard = .true. + if (iprint .eq. 1) write(6,*) 'standard' + else if (v2 .gt. vstar + osmall) then + lvacuum = .true. + if (iprint .eq. 1) write(6,*) 'vacuum' + end if + +! two apparent singularies, book's notation for omega2 and omega3 +! resetting denom2 and denom3 is harmless as these two singularities are explicitely addressed + + lomega2 = .false. + lomega3 = .false. + + if (abs(denom2) .le. osmall) then + if (iprint .eq. 1) write(6,'(a,1pe12.4)') 'omega = omega2 ',denom2 + lomega2 = .true. + denom2 = 1.0q-8 + + else if (abs(denom3) .le. osmall) then + if (iprint .eq. 1) write(6,'(a,1pe12.4)') 'omega = omega3 ',denom3 + lomega3 = .true. + denom3 = 1.0q-8 + end if + + +! various exponents, kamm equations 42-47 +! in terms of book's notation: +! a0=beta6 a1=beta1 a2=-beta2 a3=beta3 a4=beta4 a5=-beta5 + + a0 = 2.0q0/xg2 + a2 = -gamm1/denom2 + a1 = xg2*gamma/(2.0q0 + xgeom*gamm1) * & + (((2.0q0*(xgeom*(2.0q0-gamma) - omega))/(gamma*xg2*xg2))-a2) + a3 = (xgeom - omega) / denom2 + a4 = xg2 * (xgeom - omega) * a1 /denom3 + a5 = (omega*gamp1 - 2.0q0*xgeom)/denom3 + + +! frequent combinations, kamm equations 33-37 + a_val = 0.25q0 * xg2 * gamp1 + b_val = gpogm + c_val = 0.5q0 * xg2 * gamma + d_val = (xg2 * gamp1)/(xg2*gamp1 - 2.0q0*(2.0q0 + xgeom*gamm1)) + e_val = 0.5q0 * (2.0q0 + xgeom * gamm1) + + + +! evaluate the energy integrals +! the transitional case can be done by hand; save some cpu cycles +! kamm equations 80, 81, and 85 + + if (ltransitional) then + + eval2 = gamp1/(xgeom*(gamm1*xgeom + 2.0q0)**2) + eval1 = 2.0q0/gamm1 * eval2 + alpha = gpogm * 2**(xgeom)/(xgeom*(gamm1*xgeom + 2.0q0)**2) + if (int(xgeom) .ne. 1) alpha = pi * alpha + + + +! for the standard or vacuum cases +! v0 = post-shock origin v0 and vv = vacuum boundary vv +! set the radius corespondin to vv to zero for now +! kamm equations 18, and 28. + + else + v0 = 2.0q0 / (xg2 * gamma) + vv = 2.0q0 / xg2 + rvv = 0.0q0 + if (lstandard) vmin = v0 + if (lvacuum) vmin = vv + + + +! the first energy integral +! in the standard case the term (c_val*v - 1) might be singular at v=vmin + + if (lstandard) then + gam_int = a3 - a2*xg2 - 1.0q0 + if (gam_int .ge. 0) then + call qromo(efun01,vmin,v2,eps,eval1,midpnt) + else + gam_int = abs(gam_int) + call qromo(efun01,vmin,v2,eps,eval1,midpowl) + end if + +! in the vacuum case the term (1 - c_val/gamma*v) might be singular at v=vmin + + else if (lvacuum) then + gam_int = a5 + if (gam_int .ge. 0) then + call qromo(efun01,vmin,v2,eps,eval1,midpnt) + else + gam_int = abs(gam_int) + call qromo(efun01,vmin,v2,eps,eval1,midpowl2) + end if + end if + + + +! the second energy integral +! in the standard case the term (c_val*v - 1) might be singular at v=vmin + + if (lstandard) then + gam_int = a3 - a2*xg2 - 2.0q0 + if (gam_int .ge. 0) then + call qromo(efun02,vmin,v2,eps,eval2,midpnt) + else + gam_int = abs(gam_int) + call qromo(efun02,vmin,v2,eps,eval2,midpowl) + end if + +! in the vacuum case the term (1 - c_val/gamma*v) might be singular at v=vmin + + else if (lvacuum) then + gam_int = a5 + if (gam_int .ge. 0) then + call qromo(efun02,vmin,v2,eps,eval2,midpnt) + else + gam_int = abs(gam_int) + call qromo(efun02,vmin,v2,eps,eval2,midpowl2) + end if + end if + + +! kamm, bolstad & timmes equation 57 for alpha + if (int(xgeom) .eq. 1) then + +! bug reported Irina Sagert march 2019 +! alpha = 0.5q0*eval1 + eval2/gamm1 + + alpha = eval1 + 2.0q0*eval2/gamm1 + + else + alpha = (xgeom - 1.0q0) * pi * (eval1 + 2.0q0 * eval2/gamm1) + end if + end if + + +! write what we have for the energy integrals + if (iprint .eq. 1) & + write(6,88) 'xgeom =',xgeom,'gamma=',gamma, & + 'omega =',omega,'alpha =',alpha, & + 'j1 =',eval1,'j2 =',eval2 + + + + +! immediate post-shock values +! kamm page 14 or equations 14, 16, 5, 13 +! r2 = shock position, u2 = shock speed, rho1 = pre-shock density, +! u2 = post-shock material speed, rho2 = post-shock density, +! p2 = post-shock pressure, e2 = post-shoock specific internal energy, +! and cs2 = post-shock sound speed + + r2 = (eblast/(alpha*rho0))**(1.0q0/xg2) * time**(2.0q0/xg2) + rshock = r2 + us = (2.0q0/xg2) * r2 / time + rho1 = rho0 * r2**(-omega) + u2 = 2.0q0 * us / gamp1 + rho2 = gpogm * rho1 + p2 = 2.0q0 * rho1 * us**2 / gamp1 + e2 = p2/(gamm1*rho2) + cs2 = sqrt(gamma*p2/rho2) + + +! find the radius corresponding to vv + if (lvacuum) then + vwant = vv + rvv = zeroin(0.0q0,r2,sed_r_find,eps2) + rvacuum = rvv + end if + +! write a summary + if (lstandard .and. iprint .eq. 1) & + write(6,88) 'r2 =',r2,'rho2 =',rho2, & + 'u2 =',u2,'e2 =',e2, & + 'p2 =',p2,'cs2 =',cs2 + + if (lvacuum .and. iprint .eq. 1) & + write(6,88) & + 'rv =',rvv, & + 'r2 =',r2,'rho2 =',rho2, & + 'u2 =',u2,'e2 =',e2, & + 'p2 =',p2,'cs2 =',cs2 + + + + +! now start the loop over spatial positions + do i=1,nstep + rwant = xpos(i) + + +! if we are upstream from the shock front + if (rwant .gt. r2) then + den(i) = rho0 * rwant**(-omega) + vel(i) = vel0 + pres(i) = pres0 + ener(i) = ener0 + cs(i) = cs0 + + +! if we are between the origin and the shock front +! find the correct similarity value for this radius in the standard or vacuum cases + else + if (lstandard) then + vat = zeroin(0.90q0*v0,v2,sed_v_find,eps2) + else if (lvacuum) then + vat = zeroin(v2,1.2q0*vv,sed_v_find,eps2) + end if + +! the physical solution + call sedov_funcs(vat,l_fun,dlamdv,f_fun,g_fun,h_fun) + den(i) = rho2 * g_fun + vel(i) = u2 * f_fun + pres(i) = p2 * h_fun +! den(i) = g_fun +! vel(i) = f_fun +! pres(i) = h_fun + ener(i) = 0.0q0 + cs(i) = 0.0q0 + if (den(i) .ne. 0.0) then + ener(i) = pres(i) / (gamm1 * den(i)) + cs(i) = sqrt(gamma * pres(i)/den(i)) + end if + end if + +! end of loop over positions + enddo + + return + end subroutine sed_1d + + + + + + + real*16 function efun01(v) + implicit none + save + +! evaluates the first energy integrand, kamm equations 67 and 10. +! the (c_val*v - 1) term might be singular at v=vmin in the standard case. +! the (1 - c_val/gamma * v) term might be singular at v=vmin in the vacuum case. +! due care should be taken for these removable singularities by the integrator. + +! declare the pass + real*16 :: v + +! common block communication + include 'sedov3_qp.dek' + +! local variables + real*16 :: l_fun,dlamdv,f_fun,g_fun,h_fun + +! go + call sedov_funcs(v,l_fun,dlamdv,f_fun,g_fun,h_fun) + efun01 = dlamdv * l_fun**(xgeom + 1.0q0) * gpogm * g_fun * v**2 + + return + end function efun01 + + + + + + + real*16 function efun02(v) + implicit none + save + +! evaluates the second energy integrand, kamm equations 68 and 11. +! the (c_val*v - 1) term might be singular at v=vmin in the standard case. +! the (1 - c_val/gamma * v) term might be singular at v=vmin in the vacuum case. +! due care should be taken for these removable singularities by the integrator. + +! declare the pass + real*16 :: v + + +! common block communication + include 'sedov3_qp.dek' + +! local variables + real*16 :: l_fun,dlamdv,f_fun,g_fun,h_fun,z + +! go + call sedov_funcs(v,l_fun,dlamdv,f_fun,g_fun,h_fun) + z = 8.0q0/( (xgeom + 2.0q0 - omega)**2 * gamp1) + efun02 = dlamdv * l_fun**(xgeom - 1.0q0 ) * h_fun * z + + return + end function efun02 + + + + + + + + real*16 function sed_v_find(v) + implicit none + save + +! given corresponding physical distances, find the similarity variable v +! kamm equation 38 as a root find + +! declare the pass + real*16 :: v + + +! common block communication + include 'sedov3_qp.dek' + +! local variables + real*16 :: l_fun,dlamdv,f_fun,g_fun,h_fun + + + call sedov_funcs(v,l_fun,dlamdv,f_fun,g_fun,h_fun) + sed_v_find = r2*l_fun - rwant + + return + end function sed_v_find + + + + + + real*16 function sed_r_find(r) + implicit none + save + +! given the similarity variable v, find the corresponding physical distance +! kamm equation 38 as a root find + +! declare the pass + real*16 :: r + + +! common block communication + include 'sedov3_qp.dek' + + +! local variables + real*16 :: l_fun,dlamdv,f_fun,g_fun,h_fun + + call sedov_funcs(vwant,l_fun,dlamdv,f_fun,g_fun,h_fun) + sed_r_find = r2*l_fun - r + + return + end function sed_r_find + + + + + + + + + real*16 function sed_lam_find(v) + implicit none + save + +! given the similarity variable v, find the corresponding physical distance +! kamm equation 38 as a root find + +! declare the pass + real*16 :: v + + +! common block communication + include 'sedov3_qp.dek' + + +! local variables + real*16 :: l_fun,dlamdv,f_fun,g_fun,h_fun + + + call sedov_funcs(v,l_fun,dlamdv,f_fun,g_fun,h_fun) + sed_lam_find = l_fun - xlam_want + + return + end function sed_lam_find + + + + + + + + subroutine sedov_funcs(v,l_fun,dlamdv,f_fun,g_fun,h_fun) + implicit none + save + +! given the similarity variable v, returns functions +! lambda, f, g, and h and the derivative of lambda with v dlamdv + +! although the ordinary differential equations are analytic, +! the sedov expressions appear to become singular for various +! combinations of parameters and at the lower limits of the integration +! range. all these singularies are removable and done so by this routine. + + +! declare the pass + real*16 :: v,l_fun,dlamdv,f_fun,g_fun,h_fun + + +! common block communication + include 'sedov3_qp.dek' + + +! local variables + real*16 :: x1,x2,x3,x4,dx1dv,dx2dv,dx3dv,dx4dv, & + cbag,ebag,beta0,pp1,pp2,pp3,pp4,c2,c6,y,z, & + dpp2dv + real*16, parameter :: eps = 1.0q-30 + + +! frequent combinations and their derivative with v +! kamm equation 29-32, x4 a bit different to save a divide +! x1 is book's F + + x1 = a_val * v + dx1dv = a_val + + cbag = max(eps, c_val * v - 1.0q0) + x2 = b_val * cbag + dx2dv = b_val * c_val + + ebag = 1.0q0 - e_val * v + x3 = d_val * ebag + dx3dv = -d_val * e_val + + x4 = b_val * (1.0q0 - 0.5q0 * xg2 *v) + dx4dv = -b_val * 0.5q0 * xg2 + + +! transition region between standard and vacuum cases +! kamm page 15 or equations 88-92 +! lambda = l_fun is book's zeta +! f_fun is books V, g_fun is book's D, h_fun is book's P + + if (ltransitional) then + l_fun = rwant/r2 + dlamdv = 0.0q0 + f_fun = l_fun + g_fun = l_fun**(xgeom - 2.0q0) + h_fun = l_fun**xgeom + + + +! for the vacuum case in the hole + else if (lvacuum .and. rwant .lt. rvv) then + + l_fun = 0.0q0 + dlamdv = 0.0q0 + f_fun = 0.0q0 + g_fun = 0.0q0 + h_fun = 0.0q0 + + + +! omega = omega2 = (2*(gamma -1) + xgeom)/gamma case, denom2 = 0 +! book expressions 20-22 + + else if (lomega2) then + + beta0 = 1.0q0/(2.0q0 * e_val) + pp1 = gamm1 * beta0 + c6 = 0.5q0 * gamp1 + c2 = c6/gamma + y = 1.0q0/(x1 - c2) + z = (1.0q0 - x1)*y + pp2 = gamp1 * beta0 * z + dpp2dv = -gamp1 * beta0 * dx1dv * y * (1.0q0 + z) + pp3 = (4.0q0 - xgeom - 2.0q0*gamma) * beta0 + pp4 = -xgeom * gamma * beta0 + + l_fun = x1**(-a0) * x2**(pp1) * exp(pp2) + dlamdv = (-a0*dx1dv/x1 + pp1*dx2dv/x2 + dpp2dv) * l_fun + f_fun = x1 * l_fun + g_fun = x1**(a0*omega) * x2**pp3 * x4**a5 * exp(-2.0q0*pp2) + h_fun = x1**(a0*xgeom) * x2**pp4 * x4**(1.0q0 + a5) + + + +! omega = omega3 = xgeom*(2 - gamma) case, denom3 = 0 +! book expressions 23-25 + + else if (lomega3) then + + beta0 = 1.0q0/(2.0q0 * e_val) + pp1 = a3 + omega * a2 + pp2 = 1.0q0 - 4.0q0 * beta0 + c6 = 0.5q0 * gamp1 + pp3 = -xgeom * gamma * gamp1 * beta0 * (1.0q0 - x1)/(c6 - x1) + pp4 = 2.0q0 * (xgeom * gamm1 - gamma) * beta0 + + l_fun = x1**(-a0) * x2**(-a2) * x4**(-a1) + dlamdv = -(a0*dx1dv/x1 + a2*dx2dv/x2 + a1*dx4dv/x4) * l_fun + f_fun = x1 * l_fun + g_fun = x1**(a0*omega) * x2**pp1 * x4**pp2 * exp(pp3) + h_fun = x1**(a0*xgeom) * x4**pp4 * exp(pp3) + + +! for the standard or vacuum case not in the hole +! kamm equations 38-41 + + else + l_fun = x1**(-a0) * x2**(-a2) * x3**(-a1) + dlamdv = -(a0*dx1dv/x1 + a2*dx2dv/x2 + a1*dx3dv/x3) * l_fun + f_fun = x1 * l_fun + if (x4 .eq. 0.0 .and. a5 .lt. 0.0) then + g_fun = 1.0q30 + else + g_fun = x1**(a0*omega)*x2**(a3+a2*omega)*x3**(a4+a1*omega)*x4**a5 + end if + h_fun = x1**(a0*xgeom)*x3**(a4+a1*(omega-2.0q0))*x4**(1.0q0 + a5) + end if + + return + end subroutine sedov_funcs + + + + + + + + + subroutine midpnt(func,a,b,s,n) + implicit none + save + +! this routine computes the n'th stage of refinement of an extended midpoint +! rule. func is input as the name of the function to be integrated between +! limits a and b. when called with n=1, the routine returns as s the crudest +! estimate of the integralof func from a to b. subsequent calls with n=2,3... +! improve the accuracy of s by adding 2/3*3**(n-1) addtional interior points. + +! declare the pass + external :: func + integer :: n + real*16 :: func,a,b,s + +! local variables + + integer :: it,j + real*16 :: tnm,del,ddel,x,asum + + if (n.eq.1) then + s = (b-a) * func(0.5q0*(a+b)) + else + it = 3**(n-2) + tnm = it + del = (b-a)/(3.0q0*tnm) + ddel = del + del + x = a + (0.5q0 * del) + asum = 0.0q0 + do j=1,it + asum = asum + func(x) + x = x + ddel + asum = asum + func(x) + x = x + del + enddo + s = (s + ((b-a) * asum/tnm)) / 3.0q0 + end if + return + end subroutine midpnt + + + + + + + + subroutine midpowl(funk,aa,bb,s,n) + implicit none + save + +! this routine is an exact replacement for midpnt, except that it allows for +! an integrable power-law singularity of the form (x - a)**(-gam_int) +! at the lower limit aa for 0 < gam_int < 1. + +! declare the pass + external :: funk + integer :: n + real*16 :: funk,aa,bb,s + +! local variables + integer :: it,j + real*16 :: func,a,b,tnm,del,ddel,x,asum + + +! common block communication + real*16 gam_int + common /cmidp/ gam_int + + +! statement function, recipe equation 4.4.3 + func(x) = 1.0q0/(1.0q0 - gam_int) * x**(gam_int/(1.0q0 - gam_int)) & + * funk(x**(1.0q0/(1.0q0 - gam_int)) + aa) + + b = (bb - aa)**(1.0q0 - gam_int) + a = 0.0q0 + +! now exactly as midpnt + if (n .eq. 1) then + s = (b-a) * func(0.5q0*(a+b)) + else + it = 3**(n-2) + tnm = it + del = (b-a)/(3.0q0*tnm) + ddel = del + del + x = a + (0.5q0 * del) + asum = 0.0q0 + do j=1,it + asum = asum + func(x) + x = x + ddel + asum = asum + func(x) + x = x + del + enddo + s = (s + ((b-a) * asum/tnm)) / 3.0q0 + end if + return + end subroutine midpowl + + + + + subroutine midpowl2(funk,aa,bb,s,n) + implicit none + save + +! this routine is an exact replacement for midpnt, except that it allows for +! an integrable power-law singularity of the form (a - x)**(-gam_int) +! at the lower limit aa for 0 < gam_int < 1. + +! declare the pass + external :: funk + integer :: n + real*16 :: funk,aa,bb,s + +! local variables + + integer :: it,j + real*16 :: func,a,b,tnm,del,ddel,x,asum + +! common block communication + real*16 gam_int + common /cmidp/ gam_int + + +! statement function and recipe equation 4.4.3 + func(x) = 1.0q0/(gam_int - 1.0q0) * x**(gam_int/(1.0q0 - gam_int)) & + * funk(aa - x**(1.0q0/(1.0q0 - gam_int))) + + b = (aa - bb)**(1.0q0 - gam_int) + a = 0.0q0 + +! now exactly as midpnt + if (n .eq. 1) then + s = (b-a) * func(0.5q0*(a+b)) + else + it = 3**(n-2) + tnm = it + del = (b-a)/(3.0q0*tnm) + ddel = del + del + x = a + (0.5q0 * del) + + asum = 0.0q0 + do j=1,it + asum = asum + func(x) + x = x + ddel + asum = asum + func(x) + x = x + del + enddo + s = (s + ((b-a) * asum/tnm)) / 3.0q0 + end if + return + end subroutine midpowl2 + + + + + + + subroutine qromo(func,a,b,eps,ss,choose) + implicit none + save + +! this routine returns as s the integral of the function func from a to b +! with fractional accuracy eps. +! jmax limits the number of steps; nsteps = 3**(jmax-1) +! integration is done via romberg algorithm. + +! it is assumed the call to choose triples the number of steps on each call +! and that its error series contains only even powers of the number of steps. +! the external choose may be any of the above drivers, i.e midpnt,midinf... + +! declare the pass + external :: choose,func + real*16 :: a,b,ss,eps,func + +! local variables + integer :: i,j + integer, parameter :: jmax=15, jmaxp=jmax+1, k=5, km=k-1 + real*16 :: s(jmaxp),h(jmaxp),dss + + + h(1) = 1.0q0 + do j=1,jmax + call choose(func,a,b,s(j),j) + if (j .ge. k) then + call polint(h(j-km),s(j-km),k,0.0q0,ss,dss) + if (abs(dss) .le. eps*abs(ss)) return + end if + s(j+1) = s(j) + h(j+1) = h(j)/9.0q0 + enddo + write(6,*) 'too many steps in qromo' + return + end subroutine qromo + + + + + + + + + subroutine polint(xa,ya,n,x,y,dy) + implicit none + save + +! given arrays xa and ya of length n and a value x, this routine returns a +! value y and an error estimate dy. if p(x) is the polynomial of degree n-1 +! such that ya = p(xa) ya then the returned value is y = p(x) + +! declare the pass + integer :: n + real*16 :: xa(n),ya(n),x,y,dy + +! local variables + integer :: ns,i,m + integer, parameter :: nmax=20 + real*16 :: c(nmax),d(nmax),dif,dift,ho,hp,w,den + +! find the index ns of the closest table entry; initialize the c and d tables + ns = 1 + dif = abs(x - xa(1)) + do i=1,n + dift = abs(x - xa(i)) + if (dift .lt. dif) then + ns = i + dif = dift + end if + c(i) = ya(i) + d(i) = ya(i) + enddo + +! first guess for y + y = ya(ns) + +! for each column of the table, loop over the c's and d's and update them + ns = ns - 1 + do m=1,n-1 + do i=1,n-m + ho = xa(i) - x + hp = xa(i+m) - x + w = c(i+1) - d(i) + den = ho - hp + if (den .eq. 0.0) stop ' 2 xa entries are the same in polint' + den = w/den + d(i) = hp * den + c(i) = ho * den + enddo + +! after each column is completed, decide which correction c or d, to add +! to the accumulating value of y, that is, which path to take in the table +! by forking up or down. ns is updated as we go to keep track of where we +! are. the last dy added is the error indicator. + if (2*ns .lt. n-m) then + dy = c(ns+1) + else + dy = d(ns) + ns = ns - 1 + end if + y = y + dy + enddo + return + end subroutine polint + + + + + + + + real*16 function zeroin( ax, bx, f, tol) + implicit real*16 (a-h,o-z) + +!----------------------------------------------------------------------- +! +! This subroutine solves for a zero of the function f(x) in the +! interval ax,bx. +! +! input.. +! +! ax left endpoint of initial interval +! bx right endpoint of initial interval +! f function subprogram which evaluates f(x) for any x in +! the interval ax,bx +! tol desired length of the interval of uncertainty of the +! final result ( .ge. 0.0q0) +! +! +! output.. +! +! zeroin abcissa approximating a zero of f in the interval ax,bx +! +! +! it is assumed that f(ax) and f(bx) have opposite signs +! without a check. zeroin returns a zero x in the given interval +! ax,bx to within a tolerance 4*macheps*abs(x) + tol, where macheps +! is the relative machine precision. +! this function subprogram is a slightly modified translation of +! the algol 60 procedure zero given in richard brent, algorithms for +! minimization without derivatives, prentice - hall, inc. (1973). +! +!----------------------------------------------------------------------- + +! .. call list variables + + real*16 :: ax, bx, f, tol, a, b, c, d, e, eps, fa, fb, fc, tol1, xm, p, q, r, s + external :: f + +!---------------------------------------------------------------------- + +! +! compute eps, the relative machine precision +! + eps = 1.0q0 + 10 eps = eps/2.0q0 + tol1 = 1.0q0 + eps + if (tol1 .gt. 1.0q0) go to 10 +! +! initialization +! + a = ax + b = bx + fa = f(a) + fb = f(b) + +! write(6,'(a,1p2e40.32)') ' fa fb ',fa, fb + +! +! begin step +! + 20 c = a + fc = fa + d = b - a + e = d + 30 if ( abs(fc) .ge. abs(fb)) go to 40 + a = b + b = c + c = a + fa = fb + fb = fc + fc = fa +! +! convergence test +! + 40 tol1 = 2.0q0*eps*abs(b) + 0.5q0*tol + xm = 0.5q0*(c - b) + +! write(6,'(a,1p2e40.32)') ' xm tol1 ',abs(xm), tol1 +! write(6,*) + + if (abs(xm) .le. tol1) go to 90 + if (fb .eq. 0.0q0) go to 90 +! +! is bisection necessary? +! + if (abs(e) .lt. tol1) go to 70 + if (abs(fa) .le. abs(fb)) go to 70 +! +! is quadratic interpolation possible? +! + if (a .ne. c) go to 50 +! +! linear interpolation +! + s = fb/fa + p = 2.0q0*xm*s + q = 1.0q0 - s + go to 60 +! +! inverse quadratic interpolation +! + 50 q = fa/fc + r = fb/fc + s = fb/fa + p = s*(2.0q0*xm*q*(q - r) - (b - a)*(r - 1.0q0)) + q = (q - 1.0q0)*(r - 1.0q0)*(s - 1.0q0) +! +! adjust signs +! + 60 if (p .gt. 0.0q0) q = -q + p = abs(p) +! +! is interpolation acceptable? +! + if ((2.0q0*p) .ge. (3.0q0*xm*q - abs(tol1*q))) go to 70 + if (p .ge. abs(0.5q0*e*q)) go to 70 + e = d + d = p/q + go to 80 +! +! bisection +! + 70 d = xm + e = d +! +! complete step +! + 80 a = b + fa = fb + if (abs(d) .gt. tol1) b = b + d + if (abs(d) .le. tol1) b = b + Sign(tol1, xm) + fb = f(b) + if ((fb*(fc/abs(fc))) .gt. 0.0q0) go to 20 + go to 30 +! +! done +! + 90 zeroin = b + + return + end function zeroin + + + + + + + + + + real*16 function value(string) + implicit none + save + + +! this routine takes a character string and converts it to a real number. +! on error during the conversion, a fortran stop is issued + +! declare + logical :: pflag + character*(*) :: string + integer :: noblnk,long,ipoint,power,psign,j,z,i + integer, parameter :: iten = 10 + real*16 :: x,sign,factor,temp + real*16, parameter :: rten = 10.0q0 + character*1, parameter :: plus = '+' , minus = '-' , decmal = '.' , & + blank = ' ' , & + se = 'e' , sd = 'd' , sq = 'q' , & + se1 = 'E' , sd1 = 'D' , sq1 ='Q' +! initialize + x = 0.0q0 + sign = 1.0q0 + factor = rten + pflag = .false. + noblnk = 0 + power = 0 + psign = 1 + long = len(string) + + +! remove any leading blanks and get the sign of the number + do z = 1,7 + noblnk = noblnk + 1 + if ( string(noblnk:noblnk) .eq. blank) then + if (noblnk .gt. 6 ) goto 30 + else + if (string(noblnk:noblnk) .eq. plus) then + noblnk = noblnk + 1 + else if (string(noblnk:noblnk) .eq. minus) then + noblnk = noblnk + 1 + sign = -1.0q0 + end if + goto 10 + end if + enddo + + +! main number conversion loop + 10 continue + do i = noblnk,long + ipoint = i + 1 + + +! if a blank character then we are done + if ( string(i:i) .eq. blank ) then + x = x * sign + value = x + return + + +! if an exponent character, process the whole exponent, and return + else if (string(i:i) .eq. se .or. & + string(i:i) .eq. sd .or. & + string(i:i) .eq. sq .or. & + string(i:i) .eq. se1 .or. & + string(i:i) .eq. sd1 .or. & + string(i:i) .eq. sq1 ) then + if (x .eq. 0.0 .and. ipoint.eq.2) x = 1.0q0 + if (sign .eq. -1.0 .and. ipoint.eq.3) x = 1.0q0 + if (string(ipoint:ipoint) .eq. plus) ipoint = ipoint + 1 + if (string(ipoint:ipoint) .eq. minus) then + ipoint = ipoint + 1 + psign = -1 + end if + do z = ipoint,long + if (string(z:z) .eq. blank) then + x = sign * x * rten**(power*psign) + value = x + return + else + j = ichar(string(z:z)) - 48 + if ( (j.lt.0) .or. (j.gt.9) ) goto 30 + power= (power * iten) + j + end if + enddo + + +! if an ascii number character, process ie + else if (string(i:i) .ne. decmal) then + j = ichar(string(i:i)) - 48 + if ( (j.lt.0) .or. (j.gt.9) ) goto 30 + if (.not.(pflag) ) then + x = (x*rten) + j + else + temp = j + x = x + (temp/factor) + factor = factor * rten + goto 20 + end if + +! must be a decimal point if none of the above +! check that there are not two decimal points! + else + if (pflag) goto 30 + pflag = .true. + end if + 20 continue + end do + +! if we got through the do loop ok, then we must be done + x = x * sign + value = x + return + + +! error processing the number + 30 write(6,40) long,string(1:long) + 40 format(' error converting the ',i4,' characters ',/, & + ' >',a,'< ',/, & + ' into a real number in function value') + stop ' error in routine value' + end function value + diff --git a/src/sedov_blast/sedov_plot.py b/src/sedov_blast/sedov_plot.py new file mode 100755 index 0000000..533e0f4 --- /dev/null +++ b/src/sedov_blast/sedov_plot.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python3 + +# -*- coding: utf-8 -*- + +import sys, getopt +import argparse + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import rc +rc('text', usetex=True) + +import configparser +import subprocess + +def sedov_plot(ini_filename): + + config = configparser.ConfigParser() + config.read(ini_filename) + + nbins=config.get('blast','num_radial_bins', fallback=1000) + tEnd=config.get('run','tEnd', fallback=0.1) + + omega=0 + gamma=1.4 + output_filename='sedov_blast_analytical_2d.dat' + Eblast = 0.311357 + dim=2 + + print('sedov args : nbins={} Eblast={} dim={} omega={} gamma={} tEnd={} output_filename={}'.format(nbins,Eblast,dim,omega,gamma,tEnd,output_filename)) + + + # compute analytical solution + subprocess.run(["./sedov_fortran/sedov3_qp", str(nbins), str(Eblast), str(dim), str(omega), str(gamma), str(tEnd), str(output_filename)]) + + # load analytical solution + sedov_ana=np.loadtxt(output_filename, skiprows=2) + sedov_ana_r = sedov_ana[:,1] + sedov_ana_density = sedov_ana[:,2] + + # load numerical solution + sedov_num_r = np.load('sedov_blast_radial_distances.npy') + sedov_num_density = np.load('sedov_blast_density_profile.npy') + + fig = plt.figure() + plt.plot(sedov_ana_r, sedov_ana_density, 'go-', label='analytical') + plt.plot(sedov_num_r, sedov_num_density, 'xb-', label='numerical') + plt.legend() + plt.title('Sedov blast at tEnd={} with gamma={}'.format(tEnd,gamma)) + plt.show() + +############################################################################### +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description='Display sedov density plots.') + parser.add_argument('--ini', type=str, default='test_sedov_blast_2d.ini', help='kanop ini parameter file') + args = parser.parse_args() + + sedov_plot(args.ini) diff --git a/src/sedov_blast/test_sedov_blast_2d.ini b/src/sedov_blast/test_sedov_blast_2d.ini new file mode 100644 index 0000000..1ea2977 --- /dev/null +++ b/src/sedov_blast/test_sedov_blast_2d.ini @@ -0,0 +1,64 @@ +[run] +tEnd=0.5 +nStepmax=3000 + +# noutput equals -1 means we dump data at every time steps +nOutput=1000 + +[mesh] +nx=256 +ny=256 + +xmin=0.0 +xmax=1.0 + +ymin=0.0 +ymax=1.0 + +boundary_type_xmin=2 +boundary_type_xmax=2 + +boundary_type_ymin=2 +boundary_type_ymax=2 + +[hydro] +gamma0=1.4 +cfl=0.8 +niter_riemann=10 +iorder=2 +slope_type=2 +problem=blast +riemann=hllc + +[blast] +density_in=1.0 +density_out=1.0 +# +# pressure inside should be +# (gamma0-1)*epsilon_sedov/volume +# +# note that volume must be adjusted to the discretize ball (if we don't results will be biased) +# +# where epsilon_sedov is +# - 1d : 0.0673185, volume = 2r +# - 2d : 0.311357, volume = pi r^2 +# - 3d : 0.851072, volume = 4/3*pi*r^3 +# +# for reference see: +# - https://cococubed.com/papers/kamm_2000.pdf +# - https://cococubed.com/research_pages/sedov.shtml +# +# Please note that pressure inside will be adjusted so that total energy inside +# given by total_energy_inside +pressure_in=1 +pressure_out=1e-7 +radius=0.01 +total_energy_inside=0.311357 +compute_radial_profile=yes +num_radial_bins=200 + +[output] +outputPrefix=sedov_blast_2d + +[other] +implementationVersion=0