diff --git a/include/openmc/tallies/tally_filter_legendre.h b/include/openmc/tallies/tally_filter_legendre.h new file mode 100644 index 0000000000..f70cdf663d --- /dev/null +++ b/include/openmc/tallies/tally_filter_legendre.h @@ -0,0 +1,57 @@ +#ifndef OPENMC_TALLY_FILTER_LEGENDRE_H +#define OPENMC_TALLY_FILTER_LEGENDRE_H + +#include + +#include "openmc/cell.h" +#include "openmc/error.h" +#include "openmc/math_functions.h" +#include "openmc/tallies/tally_filter.h" + + +namespace openmc { + +class LegendreFilter : public TallyFilter +{ +public: + virtual std::string type() const override {return "legendre";} + + virtual ~LegendreFilter() override = default; + + virtual void + from_xml(pugi::xml_node node) override + { + order_ = std::stoi(get_node_value(node, "order")); + n_bins_ = order_ + 1; + } + + virtual void + get_all_bins(Particle* p, int estimator, TallyFilterMatch& match) + const override + { + double wgt[n_bins_]; + calc_pn_c(order_, p->mu, wgt); + for (int i = 0; i < n_bins_; i++) { + match.bins.push_back(i + 1); + match.weights.push_back(wgt[i]); + } + } + + virtual void + to_statepoint(hid_t filter_group) const override + { + TallyFilter::to_statepoint(filter_group); + write_dataset(filter_group, "order", order_); + } + + virtual std::string + text_label(int bin) const override + { + return "Legendre expansion, P" + std::to_string(bin - 1); + } + + int order_; +}; + +} // namespace openmc +#endif // OPENMC_TALLY_FILTER_LEGENDRE_H diff --git a/include/openmc/tallies/tally_filter_sph_harm.h b/include/openmc/tallies/tally_filter_sph_harm.h new file mode 100644 index 0000000000..7c7ef00d82 --- /dev/null +++ b/include/openmc/tallies/tally_filter_sph_harm.h @@ -0,0 +1,106 @@ +#ifndef OPENMC_TALLY_FILTER_SPHERICAL_HARMONICS_H +#define OPENMC_TALLY_FILTER_SPHERICAL_HARMONICS_H + +#include + +#include "openmc/error.h" +#include "openmc/tallies/tally_filter.h" + + +namespace openmc { + +//TODO: those integer values are not needed when Fortran interop is removed +enum class SphericalHarmonicsCosine { + scatter = 1, particle = 2 +}; + +class SphericalHarmonicsFilter : public TallyFilter +{ +public: + virtual std::string type() const override {return "sphericalharmonics";} + + virtual ~SphericalHarmonicsFilter() override = default; + + virtual void + from_xml(pugi::xml_node node) override + { + order_ = std::stoi(get_node_value(node, "order")); + n_bins_ = (order_ + 1) * (order_ + 1); + + if (check_for_node(node, "cosine")) { + auto cos = get_node_value(node, "cosine", true); + std::cout << cos << "\n"; + std::cout << (cos == "scatter") << " " << (cos == "particle") << "\n"; + if (cos == "scatter") { + cosine_ = SphericalHarmonicsCosine::scatter; + } else if (cos == "particle") { + cosine_ = SphericalHarmonicsCosine::particle; + } else { + std::stringstream err_msg; + err_msg << "Unrecognized cosine type, \"" << cos + << "\" in spherical harmonics filter"; + fatal_error(err_msg); + } + } + } + + virtual void + get_all_bins(Particle* p, int estimator, TallyFilterMatch& match) + const override + { + // Determine cosine term for scatter expansion if necessary + double wgt[order_ + 1]; + if (cosine_ == SphericalHarmonicsCosine::scatter) { + calc_pn_c(order_, p->mu, wgt); + } else { + for (int i = 0; i < order_ + 1; i++) wgt[i] = 1; + } + + // Find the Rn,m values + double rn[n_bins_]; + calc_rn_c(order_, p->last_uvw, rn); + + int j = 0; + for (int n = 0; n < order_ + 1; n++) { + // Calculate n-th order spherical harmonics for (u,v,w) + int num_nm = 2*n + 1; + + // Append the matching (bin,weight) for each moment + for (int i = 0; i < num_nm; i++) { + match.weights.push_back(wgt[n] * rn[j]); + match.bins.push_back(++j); + } + } + } + + virtual void + to_statepoint(hid_t filter_group) const override + { + TallyFilter::to_statepoint(filter_group); + write_dataset(filter_group, "order", order_); + if (cosine_ == SphericalHarmonicsCosine::scatter) { + write_dataset(filter_group, "cosine", "scatter"); + } else { + write_dataset(filter_group, "cosine", "particle"); + } + } + + virtual std::string + text_label(int bin) const override + { + std::stringstream out; + for (int n = 0; n < order_ + 1; n++) { + if (bin <= (n + 1) * (n + 1)) { + int m = (bin - n*n - 1) - n; + out << "Spherical harmonic expansion, Y" << n << "," << m; + return out.str(); + } + } + } + + int order_; + SphericalHarmonicsCosine cosine_ {SphericalHarmonicsCosine::particle}; +}; + +} // namespace openmc +#endif // OPENMC_TALLY_FILTER_SPHERICAL_HARMONICS_H diff --git a/src/tallies/tally_filter.cpp b/src/tallies/tally_filter.cpp index d66505d470..b3b2a72321 100644 --- a/src/tallies/tally_filter.cpp +++ b/src/tallies/tally_filter.cpp @@ -8,11 +8,13 @@ #include "openmc/tallies/tally_filter_cellborn.h" #include "openmc/tallies/tally_filter_cellfrom.h" #include "openmc/tallies/tally_filter_distribcell.h" +#include "openmc/tallies/tally_filter_legendre.h" #include "openmc/tallies/tally_filter_material.h" #include "openmc/tallies/tally_filter_mesh.h" #include "openmc/tallies/tally_filter_meshsurface.h" #include "openmc/tallies/tally_filter_mu.h" #include "openmc/tallies/tally_filter_polar.h" +#include "openmc/tallies/tally_filter_sph_harm.h" #include "openmc/tallies/tally_filter_surface.h" #include "openmc/tallies/tally_filter_universe.h" @@ -96,6 +98,8 @@ extern "C" { tally_filters.push_back(new CellFromFilter()); } else if (type_ == "distribcell") { tally_filters.push_back(new DistribcellFilter()); + } else if (type_ == "legendre") { + tally_filters.push_back(new LegendreFilter()); } else if (type_ == "material") { tally_filters.push_back(new MaterialFilter()); } else if (type_ == "mesh") { @@ -108,6 +112,8 @@ extern "C" { tally_filters.push_back(new PolarFilter()); } else if (type_ == "surface") { tally_filters.push_back(new SurfaceFilter()); + } else if (type_ == "sphericalharmonics") { + tally_filters.push_back(new SphericalHarmonicsFilter()); } else if (type_ == "universe") { tally_filters.push_back(new UniverseFilter()); } else { @@ -149,6 +155,12 @@ extern "C" { *n = filt->cells_.size(); } + int legendre_filter_get_order(LegendreFilter* filt) + {return filt->order_;} + + void legendre_filter_set_order(LegendreFilter* filt, int order) + {filt->order_ = order;} + void material_filter_get_bins(MaterialFilter* filt, int32_t** bins, int32_t* n) { @@ -164,7 +176,7 @@ extern "C" { for (int i = 0; i < n; i++) filt->materials_[i] = bins[i]; filt->n_bins_ = filt->materials_.size(); filt->map_.clear(); - for (int i = 0; i < n; i++) filt->map_[bins[i]] = i; + for (int i = 0; i < n; i++) filt->map_[filt->materials_[i]] = i; } int mesh_filter_get_mesh(MeshFilter* filt) {return filt->mesh_;} @@ -186,6 +198,22 @@ extern "C" { filt->n_bins_ = 4 * meshes[mesh]->n_dimension_; for (auto dim : meshes[mesh]->shape_) filt->n_bins_ *= dim; } + + int sphharm_filter_get_order(SphericalHarmonicsFilter* filt) + {return filt->order_;} + + int sphharm_filter_get_cosine(SphericalHarmonicsFilter* filt) + {return static_cast(filt->cosine_);} + + void + sphharm_filter_set_order(SphericalHarmonicsFilter* filt, int order) + { + filt->order_ = order; + filt->n_bins_ = (order + 1) * (order + 1); + } + + void sphharm_filter_set_cosine(SphericalHarmonicsFilter* filt, int cosine) + {filt->cosine_ = static_cast(cosine);} } } // namespace openmc diff --git a/src/tallies/tally_filter_legendre.F90 b/src/tallies/tally_filter_legendre.F90 index 32d1b45066..60bb59c765 100644 --- a/src/tallies/tally_filter_legendre.F90 +++ b/src/tallies/tally_filter_legendre.F90 @@ -2,14 +2,8 @@ module tally_filter_legendre use, intrinsic :: ISO_C_BINDING - use constants use error - use hdf5_interface - use math, only: calc_pn - use particle_header, only: Particle - use string, only: to_str use tally_filter_header - use xml_interface implicit none private @@ -20,63 +14,11 @@ module tally_filter_legendre ! LEGENDREFILTER gives Legendre moments of the change in scattering angle !=============================================================================== - type, public, extends(TallyFilter) :: LegendreFilter - integer(C_INT) :: order - contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label + type, public, extends(CppTallyFilter) :: LegendreFilter end type LegendreFilter contains -!=============================================================================== -! LegendreFilter methods -!=============================================================================== - - subroutine from_xml(this, node) - class(LegendreFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - ! Get specified order - call get_node_value(node, "order", this % order) - this % n_bins = this % order + 1 - end subroutine from_xml - - subroutine get_all_bins(this, p, estimator, match) - class(LegendreFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - real(C_DOUBLE) :: wgt(this % n_bins) - - call calc_pn(this % order, p % mu, wgt) - do i = 1, this % n_bins - call match % bins_push_back(i) - call match % weights_push_back(wgt(i)) - end do - end subroutine get_all_bins - - subroutine to_statepoint(this, filter_group) - class(LegendreFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "legendre") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "order", this % order) - end subroutine to_statepoint - - function text_label(this, bin) result(label) - class(LegendreFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Legendre expansion, P" // trim(to_str(bin - 1)) - end function text_label - !=============================================================================== ! C API FUNCTIONS !=============================================================================== @@ -87,11 +29,19 @@ contains integer(C_INT), intent(out) :: order integer(C_INT) :: err + interface + function legendre_filter_get_order(filt) result(order) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT) :: order + end function + end interface + err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (LegendreFilter) - order = f % order + order = legendre_filter_get_order(f % ptr) class default err = E_INVALID_TYPE call set_errmsg("Tried to get order on a non-expansion filter.") @@ -106,12 +56,20 @@ contains integer(C_INT), value :: order integer(C_INT) :: err + interface + subroutine legendre_filter_set_order(filt, order) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT), value :: order + end subroutine legendre_filter_set_order + end interface + err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (LegendreFilter) - f % order = order - f % n_bins = order + 1 + call legendre_filter_set_order(f % ptr, order) + f % n_bins = f % n_bins_cpp() class default err = E_INVALID_TYPE call set_errmsg("Tried to set order on a non-expansion filter.") diff --git a/src/tallies/tally_filter_sph_harm.F90 b/src/tallies/tally_filter_sph_harm.F90 index 8fba715648..59120940c7 100644 --- a/src/tallies/tally_filter_sph_harm.F90 +++ b/src/tallies/tally_filter_sph_harm.F90 @@ -2,14 +2,9 @@ module tally_filter_sph_harm use, intrinsic :: ISO_C_BINDING - use constants use error - use hdf5_interface - use math, only: calc_pn, calc_rn - use particle_header, only: Particle use string, only: to_str, to_lower, to_f_string use tally_filter_header - use xml_interface implicit none private @@ -26,112 +21,25 @@ module tally_filter_sph_harm ! tally score !=============================================================================== - type, public, extends(TallyFilter) :: SphericalHarmonicsFilter - integer :: order - integer :: cosine + type, public, extends(CppTallyFilter) :: SphericalHarmonicsFilter contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label + procedure :: cosine end type SphericalHarmonicsFilter contains -!=============================================================================== -! SphericalHarmonicsFilter methods -!=============================================================================== - - subroutine from_xml(this, node) - class(SphericalHarmonicsFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - character(MAX_WORD_LEN) :: temp_str - - ! Get specified order - call get_node_value(node, "order", this % order) - this % n_bins = (this % order + 1)**2 - - ! Determine how cosine term is to be treated - if (check_for_node(node, "cosine")) then - call get_node_value(node, "cosine", temp_str) - select case (to_lower(temp_str)) - case ('scatter') - this % cosine = COSINE_SCATTER - case ('particle') - this % cosine = COSINE_PARTICLE - end select - else - this % cosine = COSINE_PARTICLE - end if - end subroutine from_xml - - subroutine get_all_bins(this, p, estimator, match) - class(SphericalHarmonicsFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i, j, n - integer :: num_nm - real(C_DOUBLE) :: wgt(this % order + 1) - real(C_DOUBLE) :: rn(this % n_bins) - - ! Determine cosine term for scatter expansion if necessary - if (this % cosine == COSINE_SCATTER) then - call calc_pn(this % order, p % mu, wgt) - else - wgt = ONE - end if - - ! Find the Rn,m values - call calc_rn(this % order, p % last_uvw, rn) - - j = 0 - do n = 0, this % order - ! Calculate n-th order spherical harmonics for (u,v,w) - num_nm = 2*n + 1 - - ! Append matching (bin,weight) for each moment - do i = 1, num_nm - j = j + 1 - call match % bins_push_back(j) - call match % weights_push_back(wgt(n + 1) * rn(j)) - end do - end do - - end subroutine get_all_bins - - subroutine to_statepoint(this, filter_group) - class(SphericalHarmonicsFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "sphericalharmonics") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "order", this % order) - if (this % cosine == COSINE_SCATTER) then - call write_dataset(filter_group, "cosine", "scatter") - else - call write_dataset(filter_group, "cosine", "particle") - end if - end subroutine to_statepoint - - function text_label(this, bin) result(label) - class(SphericalHarmonicsFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - integer :: n, m - - do n = 0, this % order - if (bin <= (n + 1)**2) then - m = (bin - n**2 - 1) - n - label = "Spherical harmonic expansion, Y" // trim(to_str(n)) // & - "," // trim(to_str(m)) - exit - end if - end do - end function text_label + function cosine(this) result(val) + class(SphericalHarmonicsFilter) :: this + integer :: val + interface + function sphharm_filter_get_cosine(filt) result(val) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT) :: val + end function + end interface + val = sphharm_filter_get_cosine(this % ptr) + end function cosine !=============================================================================== ! C API FUNCTIONS @@ -143,11 +51,19 @@ contains integer(C_INT), intent(out) :: order integer(C_INT) :: err + interface + function sphharm_filter_get_order(filt) result(order) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT) :: order + end function sphharm_filter_get_order + end interface + err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (SphericalHarmonicsFilter) - order = f % order + order = sphharm_filter_get_order(f % ptr) class default err = E_INVALID_TYPE call set_errmsg("Not a spherical harmonics filter.") @@ -169,7 +85,7 @@ contains if (err == 0) then select type (f => filters(index) % obj) type is (SphericalHarmonicsFilter) - select case (f % cosine) + select case (f % cosine()) case (COSINE_SCATTER) cosine_ = 'scatter' case (COSINE_PARTICLE) @@ -196,12 +112,19 @@ contains integer(C_INT), value :: order integer(C_INT) :: err + interface + subroutine sphharm_filter_set_order(filt, order) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT), value :: order + end subroutine + end interface + err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (SphericalHarmonicsFilter) - f % order = order - f % n_bins = (order + 1)**2 + call sphharm_filter_set_order(f % ptr, order) class default err = E_INVALID_TYPE call set_errmsg("Not a spherical harmonics filter.") @@ -218,6 +141,14 @@ contains character(:), allocatable :: cosine_ + interface + subroutine sphharm_filter_set_cosine(filt, cosine) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT), value :: cosine + end subroutine + end interface + ! Convert C string to Fortran string cosine_ = to_f_string(cosine) @@ -227,9 +158,9 @@ contains type is (SphericalHarmonicsFilter) select case (cosine_) case ('scatter') - f % cosine = COSINE_SCATTER + call sphharm_filter_set_cosine(f % ptr, COSINE_SCATTER) case ('particle') - f % cosine = COSINE_PARTICLE + call sphharm_filter_set_cosine(f % ptr, COSINE_PARTICLE) end select class default diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 75779d1ad3..4bb504b8c8 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -344,7 +344,7 @@ contains this % estimator = ESTIMATOR_ANALOG type is (SphericalHarmonicsFilter) j = FILTER_SPH_HARMONICS - if (filt % cosine == COSINE_SCATTER) then + if (filt % cosine() == COSINE_SCATTER) then this % estimator = ESTIMATOR_ANALOG end if type is (SpatialLegendreFilter)