crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Move mesh and meshsurface filter internals to C++

Browse source
This commit is contained in:
Sterling Harper 2018-10-07 15:12:04 -04:00
parent 7f2c2904e7
commit 5a73710110
7 changed files with 107 additions and 195 deletions
Download patch file Download diff file Expand all files Collapse all files

26
include/openmc/tallies/tally_filter_mesh.h
View file

@ -53,7 +53,31 @@ public:
}
}
virtual std::string text_label(int bin) const {};
virtual void
to_statepoint(hid_t filter_group) const override
{
write_dataset(filter_group, "type", "mesh");
write_dataset(filter_group, "n_bins", n_bins_);
write_dataset(filter_group, "bins", meshes[mesh_]->id_);
}
virtual std::string
text_label(int bin) const override
{
auto& mesh = *meshes[mesh_];
int n_dim = mesh.n_dimension_;
int ijk[n_dim];
mesh.get_indices_from_bin(bin, ijk);
std::stringstream out;
out << "Mesh Index (" << ijk[0];
if (n_dim > 1) out << ", " << ijk[1];
if (n_dim > 2) out << ", " << ijk[2];
out << ")";
return out.str();
}
int32_t mesh_;
};

58
include/openmc/tallies/tally_filter_meshsurface.h
View file

@ -4,6 +4,7 @@
#include <cstdint>
#include <sstream>
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/mesh.h"
#include "openmc/tallies/tally_filter.h"
@ -45,7 +46,62 @@ public:
for (auto b : match.bins) match.weights.push_back(1.0);
}
virtual std::string text_label(int bin) const {};
virtual void
to_statepoint(hid_t filter_group) const override
{
write_dataset(filter_group, "type", "meshsurface");
write_dataset(filter_group, "n_bins", n_bins_);
write_dataset(filter_group, "bins", meshes[mesh_]->id_);
}
virtual std::string
text_label(int bin) const override
{
auto& mesh = *meshes[mesh_];
int n_dim = mesh.n_dimension_;
// Get flattend mesh index and surface index.
int i_mesh = (bin - 1) / (4 * n_dim) + 1;
int i_surf = ((bin - 1) % (4 * n_dim)) + 1;
// Get mesh index part of label.
int ijk[n_dim];
mesh.get_indices_from_bin(i_mesh, ijk);
std::stringstream out;
out << "Mesh Index (" << ijk[0];
if (n_dim > 1) out << ", " << ijk[1];
if (n_dim > 2) out << ", " << ijk[2];
out << ")";
// Get surface part of label.
if (i_surf == OUT_LEFT) {
out << " Outgoing, x-min";
} else if (i_surf == IN_LEFT) {
out << " Incoming, x-min";
} else if (i_surf == OUT_RIGHT) {
out << " Outgoing, x-max";
} else if (i_surf == IN_RIGHT) {
out << " Incoming, x-max";
} else if (i_surf == OUT_BACK) {
out << " Outgoing, y-min";
} else if (i_surf == IN_BACK) {
out << " Incoming, y-min";
} else if (i_surf == OUT_FRONT) {
out << " Outgoing, y-max";
} else if (i_surf == IN_FRONT) {
out << " Incoming, y-max";
} else if (i_surf == OUT_BOTTOM) {
out << " Outgoing, z-min";
} else if (i_surf == IN_BOTTOM) {
out << " Incoming, z-min";
} else if (i_surf == OUT_TOP) {
out << " Outgoing, z-max";
} else if (i_surf == IN_TOP) {
out << " Incoming, z-max";
}
return out.str();
}
int32_t mesh_;
};

10
src/cmfd_data.F90
View file

@ -82,6 +82,14 @@ contains
real(8) :: flux ! temp variable for flux
type(RegularMesh) :: m ! pointer for mesh object
interface
function mesh_filter_get_mesh(filt) result(index_mesh) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: filt
integer(C_INT) :: index_mesh
end function mesh_filter_get_mesh
end interface
! Extract spatial and energy indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
@ -99,7 +107,7 @@ contains
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m = meshes(filt % mesh)
m = meshes(filt % mesh())
end select
! Set mesh widths

1
src/tallies/tally_filter_header.F90
View file

@ -379,6 +379,7 @@ contains
class(CppTallyFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
call this % from_xml_cpp_inner(node)
this % n_bins = this % n_bins_cpp()
end subroutine from_xml_cpp_default
subroutine get_all_bins_cpp_default(this, p, estimator, match)

94
src/tallies/tally_filter_mesh.F90
View file

@ -3,14 +3,9 @@ module tally_filter_mesh
use, intrinsic :: ISO_C_BINDING
use constants
use dict_header, only: EMPTY
use error
use mesh_header
use hdf5_interface
use particle_header, only: Particle
use string, only: to_str
use tally_filter_header
use xml_interface
implicit none
private
@ -24,77 +19,24 @@ module tally_filter_mesh
!===============================================================================
type, public, extends(CppTallyFilter) :: MeshFilter
integer :: mesh
contains
procedure :: from_xml
procedure :: to_statepoint => to_statepoint_mesh
procedure :: text_label => text_label_mesh
procedure :: mesh => get_mesh
end type MeshFilter
contains
subroutine from_xml(this, node)
class(MeshFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
integer :: id
integer :: n
integer(C_INT) :: err
call this % from_xml_cpp_inner(node)
n = node_word_count(node, "bins")
if (n /= 1) call fatal_error("Only one mesh can be &
&specified per mesh filter.")
! Determine id of mesh
call get_node_value(node, "bins", id)
! Get pointer to mesh
err = openmc_get_mesh_index(id, this % mesh)
if (err /= 0) then
call fatal_error("Could not find mesh " // trim(to_str(id)) &
// " specified on filter.")
end if
! Determine number of bins
this % n_bins = this % n_bins_cpp()
end subroutine from_xml
subroutine to_statepoint_mesh(this, filter_group)
class(MeshFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
type(RegularMesh) :: m
m = meshes(this % mesh)
call write_dataset(filter_group, "type", "mesh")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", m % id())
end subroutine to_statepoint_mesh
function text_label_mesh(this, bin) result(label)
class(MeshFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
integer, allocatable :: ijk(:)
type(RegularMesh) :: m
m = meshes(this % mesh)
allocate(ijk(m % n_dimension()))
call m % get_indices_from_bin(bin, ijk)
if (m % n_dimension() == 1) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
elseif (m % n_dimension() == 2) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension() == 3) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
end function text_label_mesh
function get_mesh(this) result(mesh)
class(MeshFilter) :: this
integer :: mesh
interface
function mesh_filter_get_mesh(filt) result(index_mesh) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: filt
integer(C_INT) :: index_mesh
end function mesh_filter_get_mesh
end interface
mesh = mesh_filter_get_mesh(this % ptr)
end function get_mesh
!===============================================================================
! C API FUNCTIONS
@ -106,19 +48,11 @@ contains
integer(C_INT32_T), intent(out) :: index_mesh
integer(C_INT) :: err
interface
function mesh_filter_get_mesh(filt) result(index_mesh) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: filt
integer(C_INT) :: index_mesh
end function mesh_filter_get_mesh
end interface
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (MeshFilter)
index_mesh = mesh_filter_get_mesh(f % ptr)
index_mesh = f % mesh()
class default
err = E_INVALID_TYPE
call set_errmsg("Tried to set mesh on a non-mesh filter.")

111
src/tallies/tally_filter_meshsurface.F90
View file

@ -5,11 +5,7 @@ module tally_filter_meshsurface
use constants
use error
use mesh_header
use hdf5_interface
use particle_header, only: Particle
use string, only: to_str
use tally_filter_header
use xml_interface
implicit none
private
@ -23,116 +19,10 @@ module tally_filter_meshsurface
!===============================================================================
type, public, extends(CppTallyFilter) :: MeshSurfaceFilter
integer :: mesh
contains
procedure :: from_xml
procedure :: to_statepoint
procedure :: text_label
end type MeshSurfaceFilter
contains
subroutine from_xml(this, node)
class(MeshSurfaceFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
integer :: id
integer :: n
integer(C_INT) :: err
call this % from_xml_cpp_inner(node)
n = node_word_count(node, "bins")
if (n /= 1) call fatal_error("Only one mesh can be &
&specified per meshsurface filter.")
! Determine id of mesh
call get_node_value(node, "bins", id)
! Get pointer to mesh
err = openmc_get_mesh_index(id, this % mesh)
if (err /= 0) then
call fatal_error("Could not find mesh " // trim(to_str(id)) &
// " specified on filter.")
end if
! Determine number of bins
this % n_bins = this % n_bins_cpp()
end subroutine from_xml
subroutine to_statepoint(this, filter_group)
class(MeshSurfaceFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
type(RegularMesh) :: m
m = meshes(this % mesh)
call write_dataset(filter_group, "type", "meshsurface")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", m % id())
end subroutine to_statepoint
function text_label(this, bin) result(label)
class(MeshSurfaceFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
integer :: i_mesh
integer :: i_surf
integer :: n_dim
integer, allocatable :: ijk(:)
type(RegularMesh) :: m
m = meshes(this % mesh)
n_dim = m % n_dimension()
allocate(ijk(n_dim))
! Get flattend mesh index and surface index
i_mesh = (bin - 1) / (4*n_dim) + 1
i_surf = mod(bin - 1, 4*n_dim) + 1
! Get mesh index part of label
call m % get_indices_from_bin(i_mesh, ijk)
if (m % n_dimension() == 1) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
elseif (m % n_dimension() == 2) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension() == 3) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
! Get surface part of label
select case (i_surf)
case (OUT_LEFT)
label = trim(label) // " Outgoing, x-min"
case (IN_LEFT)
label = trim(label) // " Incoming, x-min"
case (OUT_RIGHT)
label = trim(label) // " Outgoing, x-max"
case (IN_RIGHT)
label = trim(label) // " Incoming, x-max"
case (OUT_BACK)
label = trim(label) // " Outgoing, y-min"
case (IN_BACK)
label = trim(label) // " Incoming, y-min"
case (OUT_FRONT)
label = trim(label) // " Outgoing, y-max"
case (IN_FRONT)
label = trim(label) // " Incoming, y-max"
case (OUT_BOTTOM)
label = trim(label) // " Outgoing, z-min"
case (IN_BOTTOM)
label = trim(label) // " Incoming, z-min"
case (OUT_TOP)
label = trim(label) // " Outgoing, z-max"
case (IN_TOP)
label = trim(label) // " Incoming, z-max"
end select
end function text_label
!===============================================================================
! C API FUNCTIONS
!===============================================================================
@ -155,7 +45,6 @@ contains
if (err == 0) then
select type (f => filters(index) % obj)
type is (MeshSurfaceFilter)
index_mesh = f % mesh
index_mesh = meshsurface_filter_get_mesh(f % ptr)
class default
err = E_INVALID_TYPE

2
src/tallies/trigger.F90
View file

@ -264,7 +264,7 @@ contains
i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m = meshes(filt % mesh)
m = meshes(filt % mesh())
end select
! initialize bins array