From 967161cee87bde5a9b73e643c888276d5abd8e53 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Fri, 15 Mar 2019 12:16:49 +0300 Subject: [PATCH 01/31] test changes --- .cproject | 108 +++++++++++++++++++++ .project | 27 ++++++ .settings/org.eclipse.cdt.codan.core.prefs | 75 ++++++++++++++ include/openmc/constants.h | 2 +- include/openmc/lattice.h | 1 + src/initialize.cpp | 4 + 6 files changed, 216 insertions(+), 1 deletion(-) create mode 100644 .cproject create mode 100644 .project create mode 100644 .settings/org.eclipse.cdt.codan.core.prefs diff --git a/.cproject b/.cproject new file mode 100644 index 000000000..bdc88d84b --- /dev/null +++ b/.cproject @@ -0,0 +1,108 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/.project b/.project new file mode 100644 index 000000000..5eb05adae --- /dev/null +++ b/.project @@ -0,0 +1,27 @@ + + + openmccc + + + + + + org.eclipse.cdt.managedbuilder.core.genmakebuilder + clean,full,incremental, + + + + + org.eclipse.cdt.managedbuilder.core.ScannerConfigBuilder + full,incremental, + + + + + + org.eclipse.cdt.core.cnature + org.eclipse.cdt.core.ccnature + org.eclipse.cdt.managedbuilder.core.managedBuildNature + org.eclipse.cdt.managedbuilder.core.ScannerConfigNature + + diff --git a/.settings/org.eclipse.cdt.codan.core.prefs b/.settings/org.eclipse.cdt.codan.core.prefs new file mode 100644 index 000000000..77496f558 --- /dev/null +++ b/.settings/org.eclipse.cdt.codan.core.prefs @@ -0,0 +1,75 @@ +eclipse.preferences.version=1 +org.eclipse.cdt.codan.checkers.errnoreturn=Warning +org.eclipse.cdt.codan.checkers.errnoreturn.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"No return\\")",implicit\=>false} +org.eclipse.cdt.codan.checkers.errreturnvalue=Error +org.eclipse.cdt.codan.checkers.errreturnvalue.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Unused return value\\")"} +org.eclipse.cdt.codan.checkers.nocommentinside=-Error +org.eclipse.cdt.codan.checkers.nocommentinside.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Nesting comments\\")"} +org.eclipse.cdt.codan.checkers.nolinecomment=-Error +org.eclipse.cdt.codan.checkers.nolinecomment.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Line comments\\")"} +org.eclipse.cdt.codan.checkers.noreturn=Error +org.eclipse.cdt.codan.checkers.noreturn.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"No return value\\")",implicit\=>false} +org.eclipse.cdt.codan.internal.checkers.AbstractClassCreation=Error +org.eclipse.cdt.codan.internal.checkers.AbstractClassCreation.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Abstract class cannot be instantiated\\")"} +org.eclipse.cdt.codan.internal.checkers.AmbiguousProblem=Error +org.eclipse.cdt.codan.internal.checkers.AmbiguousProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Ambiguous problem\\")"} +org.eclipse.cdt.codan.internal.checkers.AssignmentInConditionProblem=Warning +org.eclipse.cdt.codan.internal.checkers.AssignmentInConditionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Assignment in condition\\")"} +org.eclipse.cdt.codan.internal.checkers.AssignmentToItselfProblem=Error +org.eclipse.cdt.codan.internal.checkers.AssignmentToItselfProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Assignment to itself\\")"} +org.eclipse.cdt.codan.internal.checkers.CaseBreakProblem=Warning +org.eclipse.cdt.codan.internal.checkers.CaseBreakProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"No break at end of case\\")",no_break_comment\=>"no break",last_case_param\=>false,empty_case_param\=>false,enable_fallthrough_quickfix_param\=>false} +org.eclipse.cdt.codan.internal.checkers.CatchByReference=Warning +org.eclipse.cdt.codan.internal.checkers.CatchByReference.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Catching by reference is recommended\\")",unknown\=>false,exceptions\=>()} +org.eclipse.cdt.codan.internal.checkers.CircularReferenceProblem=Error +org.eclipse.cdt.codan.internal.checkers.CircularReferenceProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Circular inheritance\\")"} +org.eclipse.cdt.codan.internal.checkers.ClassMembersInitialization=Warning +org.eclipse.cdt.codan.internal.checkers.ClassMembersInitialization.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Class members should be properly initialized\\")",skip\=>true} +org.eclipse.cdt.codan.internal.checkers.DecltypeAutoProblem=Error +org.eclipse.cdt.codan.internal.checkers.DecltypeAutoProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid 'decltype(auto)' specifier\\")"} +org.eclipse.cdt.codan.internal.checkers.FieldResolutionProblem=Error +org.eclipse.cdt.codan.internal.checkers.FieldResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Field cannot be resolved\\")"} +org.eclipse.cdt.codan.internal.checkers.FunctionResolutionProblem=Error +org.eclipse.cdt.codan.internal.checkers.FunctionResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Function cannot be resolved\\")"} +org.eclipse.cdt.codan.internal.checkers.InvalidArguments=Error +org.eclipse.cdt.codan.internal.checkers.InvalidArguments.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid arguments\\")"} +org.eclipse.cdt.codan.internal.checkers.InvalidTemplateArgumentsProblem=Error +org.eclipse.cdt.codan.internal.checkers.InvalidTemplateArgumentsProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid template argument\\")"} +org.eclipse.cdt.codan.internal.checkers.LabelStatementNotFoundProblem=Error +org.eclipse.cdt.codan.internal.checkers.LabelStatementNotFoundProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Label statement not found\\")"} +org.eclipse.cdt.codan.internal.checkers.MemberDeclarationNotFoundProblem=Error +org.eclipse.cdt.codan.internal.checkers.MemberDeclarationNotFoundProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Member declaration not found\\")"} +org.eclipse.cdt.codan.internal.checkers.MethodResolutionProblem=Error +org.eclipse.cdt.codan.internal.checkers.MethodResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Method cannot be resolved\\")"} +org.eclipse.cdt.codan.internal.checkers.NamingConventionFunctionChecker=-Info +org.eclipse.cdt.codan.internal.checkers.NamingConventionFunctionChecker.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Name convention for function\\")",pattern\=>"^[a-z]",macro\=>true,exceptions\=>()} +org.eclipse.cdt.codan.internal.checkers.NonVirtualDestructorProblem=Warning +org.eclipse.cdt.codan.internal.checkers.NonVirtualDestructorProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Class has a virtual method and non-virtual destructor\\")"} +org.eclipse.cdt.codan.internal.checkers.OverloadProblem=Error +org.eclipse.cdt.codan.internal.checkers.OverloadProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid overload\\")"} +org.eclipse.cdt.codan.internal.checkers.RedeclarationProblem=Error +org.eclipse.cdt.codan.internal.checkers.RedeclarationProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid redeclaration\\")"} +org.eclipse.cdt.codan.internal.checkers.RedefinitionProblem=Error +org.eclipse.cdt.codan.internal.checkers.RedefinitionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid redefinition\\")"} +org.eclipse.cdt.codan.internal.checkers.ReturnStyleProblem=-Warning +org.eclipse.cdt.codan.internal.checkers.ReturnStyleProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Return with parenthesis\\")"} +org.eclipse.cdt.codan.internal.checkers.ScanfFormatStringSecurityProblem=-Warning +org.eclipse.cdt.codan.internal.checkers.ScanfFormatStringSecurityProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Format String Vulnerability\\")"} +org.eclipse.cdt.codan.internal.checkers.StatementHasNoEffectProblem=Warning +org.eclipse.cdt.codan.internal.checkers.StatementHasNoEffectProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Statement has no effect\\")",macro\=>true,exceptions\=>()} +org.eclipse.cdt.codan.internal.checkers.SuggestedParenthesisProblem=Warning +org.eclipse.cdt.codan.internal.checkers.SuggestedParenthesisProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Suggested parenthesis around expression\\")",paramNot\=>false} +org.eclipse.cdt.codan.internal.checkers.SuspiciousSemicolonProblem=Warning +org.eclipse.cdt.codan.internal.checkers.SuspiciousSemicolonProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Suspicious semicolon\\")",else\=>false,afterelse\=>false} +org.eclipse.cdt.codan.internal.checkers.TypeResolutionProblem=Error +org.eclipse.cdt.codan.internal.checkers.TypeResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Type cannot be resolved\\")"} +org.eclipse.cdt.codan.internal.checkers.UnusedFunctionDeclarationProblem=Warning +org.eclipse.cdt.codan.internal.checkers.UnusedFunctionDeclarationProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Unused function declaration\\")",macro\=>true} +org.eclipse.cdt.codan.internal.checkers.UnusedStaticFunctionProblem=Warning +org.eclipse.cdt.codan.internal.checkers.UnusedStaticFunctionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Unused static function\\")",macro\=>true} +org.eclipse.cdt.codan.internal.checkers.UnusedVariableDeclarationProblem=Warning +org.eclipse.cdt.codan.internal.checkers.UnusedVariableDeclarationProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Unused variable declaration in file scope\\")",macro\=>true,exceptions\=>("@(\#)","$Id")} +org.eclipse.cdt.codan.internal.checkers.VariableResolutionProblem=Error +org.eclipse.cdt.codan.internal.checkers.VariableResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Symbol is not resolved\\")"} +org.eclipse.cdt.qt.core.qtproblem=Warning +org.eclipse.cdt.qt.core.qtproblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>false,RUN_ON_INC_BUILD\=>false,RUN_ON_FILE_OPEN\=>true,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>null} diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 87388c91a..c11f6535e 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -25,7 +25,7 @@ constexpr int VERSION_RELEASE {0}; constexpr std::array VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE}; // HDF5 data format -constexpr int HDF5_VERSION[] {2, 0}; +constexpr int HDF5_VERSION[] {2, 0};//DR change version HDF5 to work // Version numbers for binary files constexpr std::array VERSION_STATEPOINT {17, 0}; diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index 9c8b8af68..eb8e1c14f 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -266,6 +266,7 @@ public: void to_hdf5_inner(hid_t group_id) const; private: + int n_hextype_; //! Date: Fri, 15 Mar 2019 12:19:33 +0300 Subject: [PATCH 02/31] new ignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 383232296..f73309646 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,7 @@ *.pyc # Python distribution +.settings/ dist/ openmc.egg-info/ From 4b2bb55054deb234ad9eb1b6238138628f5803d9 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Mon, 8 Apr 2019 16:01:36 +0300 Subject: [PATCH 03/31] new hex first commit --- include/openmc/lattice.h | 14 +- openmc/lattice.py | 489 ++++++++++++++++++++++++++++++++++----- src/lattice.cpp | 423 +++++++++++++++++++++++++-------- src/relaxng/geometry.rnc | 1 + src/relaxng/geometry.rng | 10 + 5 files changed, 778 insertions(+), 159 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index eb8e1c14f..752e928c4 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -249,15 +249,25 @@ public: bool are_valid_indices(const int i_xyz[3]) const; + std::pair> distance(Position r, Direction u, const std::array& i_xyz) const; + //DR OX distance calculation std::pair> - distance(Position r, Direction u, const std::array& i_xyz) const; + distance_ox(Position r, Direction u, const std::array& i_xyz) const; + //DR OY distance calculation + std::pair> + distance_oy(Position r, Direction u, const std::array& i_xyz) const; std::array get_indices(Position r, Direction u) const; + Position get_local_position(Position r, const std::array i_xyz) const; bool is_valid_index(int indx) const; + //DR fill lattice with OY by default orientation + void fill_lattice_oy(std::vector univ_words); + //DR fill lattice with OX orientation + void fill_lattice_ox(std::vector univ_words); int32_t& offset(int map, const int i_xyz[3]); @@ -266,9 +276,9 @@ public: void to_hdf5_inner(hid_t group_id) const; private: - int n_hextype_; //! pitch_; //!< Lattice tile width and height }; diff --git a/openmc/lattice.py b/openmc/lattice.py index 8bf4e4757..699a7540d 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -141,6 +141,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): center = group['center'][...] pitch = group['pitch'][...] outer = group['outer'].value + orientation = group['orientation'].value universe_ids = group['universes'][...] @@ -152,63 +153,121 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): # If the Universe specified outer the Lattice is not void if outer >= 0: lattice.outer = universes[outer] + if (orientation.lower() is "oy"): + # Build array of Universe pointers for the Lattice. Note that + # we need to convert between the HDF5's square array of + # (x, alpha, z) to the Python API's format of a ragged nested + # list of (z, ring, theta). + uarray = [] + for z in range(n_axial): + # Add a list for this axial level. + uarray.append([]) + x = n_rings - 1 + a = 2*n_rings - 2 + for r in range(n_rings - 1, 0, -1): + # Add a list for this ring. + uarray[-1].append([]) - # Build array of Universe pointers for the Lattice. Note that - # we need to convert between the HDF5's square array of - # (x, alpha, z) to the Python API's format of a ragged nested - # list of (z, ring, theta). - uarray = [] - for z in range(n_axial): - # Add a list for this axial level. - uarray.append([]) - x = n_rings - 1 - a = 2*n_rings - 2 - for r in range(n_rings - 1, 0, -1): - # Add a list for this ring. - uarray[-1].append([]) + # Climb down the top-right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + x += 1 + a -= 1 - # Climb down the top-right. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - x += 1 - a -= 1 + # Climb down the right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + a -= 1 - # Climb down the right. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - a -= 1 + # Climb down the bottom-right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + x -= 1 - # Climb down the bottom-right. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - x -= 1 + # Climb up the bottom-left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + x -= 1 + a += 1 - # Climb up the bottom-left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - x -= 1 - a += 1 + # Climb up the left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + a += 1 - # Climb up the left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - a += 1 + # Climb up the top-left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + x += 1 - # Climb up the top-left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - x += 1 + # Move down to the next ring. + a -= 1 - # Move down to the next ring. - a -= 1 + # Convert the ids into Universe objects. + uarray[-1][-1] = [universes[u_id] + for u_id in uarray[-1][-1]] - # Convert the ids into Universe objects. - uarray[-1][-1] = [universes[u_id] - for u_id in uarray[-1][-1]] + # Handle the degenerate center ring separately. + u_id = universe_ids[z, a, x] + uarray[-1].append([universes[u_id]]) + else: + # Build array of Universe pointers for the Lattice. Note that + # we need to convert between the HDF5's square array of + # (alpha, y, z) to the Python API's format of a ragged nested + # list of (z, ring, theta). + uarray = [] + for z in range(n_axial): + # Add a list for this axial level. + uarray.append([]) + a = 2*n_rings - 2 + y = n_rings - 1 + for r in range(n_rings - 1, 0, -1): + # Add a list for this ring. + uarray[-1].append([]) - # Handle the degenerate center ring separately. - u_id = universe_ids[z, a, x] - uarray[-1].append([universes[u_id]]) + # Climb up the top-left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a -= 1 + y += 1 + + + # Climb the left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a -= 1 + + # Climb down the bottom-left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + y -= 1 + + # Climb down the bottom-right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a += 1 + y -= 1 + + # Climb up the right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a += 1 + + # Climb up the top-right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + y += 1 + + # Move down to the next ring. + a -= 1 + + # Convert the ids into Universe objects. + uarray[-1][-1] = [universes[u_id] + for u_id in uarray[-1][-1]] + + # Handle the degenerate center ring separately. + u_id = universe_ids[z, y, a] + uarray[-1].append([universes[u_id]]) # Add the universes to the lattice. if len(pitch) == 2: @@ -815,10 +874,10 @@ class HexLattice(Lattice): Most methods for this class use a natural indexing scheme wherein elements are assigned an index corresponding to their position relative to skewed - :math:`(x,\alpha,z)` axes as described fully in - :ref:`hexagonal_indexing`. However, note that when universes are assigned to - lattice elements using the :attr:`HexLattice.universes` property, the array - indices do not correspond to natural indices. + :math:`(x,\alpha,z)` - OY orientation or :math:`(\alpha,y,z)` - OX orientation //DR + axes as described fully in :ref:`hexagonal_indexing`. However, note that when + universes are assigned to lattice elements using the :attr:`HexLattice.universes` property, + the array indices do not correspond to natural indices. Parameters ---------- @@ -855,6 +914,8 @@ class HexLattice(Lattice): possible, where z is the axial index, r is in the ring index (starting from the outermost ring), and i is the index with a ring starting from the top and proceeding clockwise. + orientation : str by default 'OY' orientation of main lattice diagonal another option + - 'OX'//DR num_rings : int Number of radial ring positions in the xy-plane num_axial : int @@ -869,11 +930,13 @@ class HexLattice(Lattice): self._num_rings = None self._num_axial = None self._center = None + self._hextype = 0 # by Default orientaion is OY def __repr__(self): string = 'HexLattice\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) + string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', "OX" if self._hextype else "OY" )#//DR string += '{0: <16}{1}{2}\n'.format('\t# Rings', '=\t', self._num_rings) string += '{0: <16}{1}{2}\n'.format('\t# Axial', '=\t', self._num_axial) string += '{0: <16}{1}{2}\n'.format('\tCenter', '=\t', @@ -901,7 +964,12 @@ class HexLattice(Lattice): @property def num_rings(self): return self._num_rings - + @property + def orientation(self):#//DR + if self._hextype: + return "OX" + else: + return "OY" @property def num_axial(self): return self._num_axial @@ -955,6 +1023,12 @@ class HexLattice(Lattice): cv.check_length('lattice center', center, 2, 3) self._center = center + @orientation.setter + def orientation(self, orientation): + cv.check_type('orientation', orientation, str) + if orientation.lower() == "ox": + self._hextype = 1 + @Lattice.pitch.setter def pitch(self, pitch): cv.check_type('lattice pitch', pitch, Iterable, Real) @@ -1045,7 +1119,7 @@ class HexLattice(Lattice): ------- 3-tuple of int Indices of corresponding lattice element in :math:`(x,\alpha,z)` - bases + or :math:`(\alpha,y,z)`basesis numpy.ndarray Carestian coordinates of the point in the corresponding lattice element coordinate system @@ -1059,15 +1133,25 @@ class HexLattice(Lattice): else: z = point[2] - self.center[2] iz = floor(z/self.pitch[1] + 0.5*self.num_axial) - alpha = y - x/sqrt(3.) - ix = floor(x/(sqrt(0.75) * self.pitch[0])) - ia = floor(alpha/self.pitch[0]) + if self._hextype is 0:#//DR + alpha = y - x/sqrt(3.) + ix = floor(x/(sqrt(0.75) * self.pitch[0])) + ia = floor(alpha/self.pitch[0]) + else: + alpha = y - x*sqrt(3.) + ia = floor(-alpha/(sqrt(3.0) * self.pitch[0])) + iy = floor(y/(sqrt(0.75) * self.pitch[0])) + # Check four lattice elements to see which one is closest based on local # coordinates + #//DR + numbersOY=[(ix, ia, iz), (ix + 1, ia, iz), (ix, ia + 1, iz), \ + (ix + 1, ia + 1, iz)]# in case OY default orientation + numbersOX=[(ia, iy, iz), (ia + 1, iy, iz), (ia, iy + 1, iz), \ + (ia + 1, iy + 1, iz)]# in case OX orientation d_min = np.inf - for idx in [(ix, ia, iz), (ix + 1, ia, iz), (ix, ia + 1, iz), - (ix + 1, ia + 1, iz)]: + for idx in (numbersOX if self._hextype else numbersOY): #//DR p = self.get_local_coordinates(point, idx) d = p[0]**2 + p[1]**2 if d < d_min: @@ -1094,16 +1178,40 @@ class HexLattice(Lattice): system """ - x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0]) - y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0]) + if self._hextype:#//DR + x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0]) + y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0]) + else: + x = point[0] - (self.center[0] + self.pitch[0]*idx[0] + 0.5*self.pitch[0]*idx[1]) + y = point[1] - (self.center[1] + (sqrt(0.75)*self.pitch[0]*idx[1])) + if self._num_axial is None: z = point[2] else: z = point[2] - (self.center[2] + (idx[2] + 0.5 - 0.5*self.num_axial)* self.pitch[1]) return (x, y, z) + def get_universe_index(self, idx): #//DR + r"""Return index in the universes array corresponding to a lattice element index - def get_universe_index(self, idx): + Parameters + ---------- + idx : Iterable of int + Lattice element indices in the :math:`(x,\alpha,z)` coordinate + system in 'OY' orientation case, or indices in the :math:`(\alpha,y,z)` coordinate + system in 'OX' one + + Returns + ------- + 2- or 3-tuple of int + Indices used when setting the :attr:`HexLattice.universes` property + + """ + if self._hextype: + return self.get_universe_index_ox(idx) + else: + return self.get_universe_index_oy(idx) + def get_universe_index_oy(self, idx):#//DR r"""Return index in the universes array corresponding to a lattice element index Parameters @@ -1143,13 +1251,55 @@ class HexLattice(Lattice): else: return (idx[2], i_ring, i_within) + + def get_universe_index_ox(self, idx):#//DR + r"""Return index in the universes array corresponding to a lattice element index + numeration opposite clockwise + + Parameters + ---------- + idx : Iterable of int + Lattice element indices in the :math:`(\alpha,y,z)` coordinate + system + + Returns + ------- + 2- or 3-tuple of int + Indices used when setting the :attr:`HexLattice.universes` property + + """ + + # First we determine which ring the index corresponds to. + a = idx[0] + y = idx[1] + z = -a - y + g = max(abs(y), abs(a), abs(z)) + + # Next we use a clever method to figure out where along the ring we are. + i_ring = self._num_rings - 1 - g + if y >= 0: + if a >= 0: + i_within = y + else: + i_within = 2*g + z + else: + if a <= 0: + i_within = 3*g - y + else: + i_within = 5*g - z + + if self.num_axial is None: + return (i_ring, i_within) + else: + return (idx[2], i_ring, i_within) + def is_valid_index(self, idx): r"""Determine whether lattice element index is within defined range Parameters ---------- idx : Iterable of int - Lattice element indices in the :math:`(x,\alpha,z)` coordinate + Lattice element indices in the both :math:`(x,\alpha,z)` and :math:`(\alpha,y,z)` coordinate system Returns @@ -1193,6 +1343,9 @@ class HexLattice(Lattice): self._outer.create_xml_subelement(xml_element) lattice_subelement.set("n_rings", str(self._num_rings)) + # If orientation is "OX" export it to XML + if self._hextype: + lattice_subelement.set("orient", "OX") if self._num_axial is not None: lattice_subelement.set("n_axial", str(self._num_axial)) @@ -1315,10 +1468,123 @@ class HexLattice(Lattice): j += 1 lat.universes = univs return lat - def _repr_axial_slice(self, universes): """Return string representation for the given 2D group of universes. + The 'universes' argument should be a list of lists of universes where + each sub-list represents a single ring. The first list should be the + outer ring. + """ + if self._hextype: + return self._repr_axial_slice_ox(universes) + else: + return self._repr_axial_slice_oy(universes) + + def _repr_axial_slice_ox(self, universes): + """Return string representation for the given 2D group of universes in 'OX' orientation case. + + The 'universes' argument should be a list of lists of universes where + each sub-list represents a single ring. The first list should be the + outer ring. + """ + + # Find the largest universe ID and count the number of digits so we can + # properly pad the output string later. + largest_id = max([max([univ._id for univ in ring]) + for ring in universes]) + n_digits = len(str(largest_id)) + pad = ' '*n_digits + id_form = '{: ^' + str(n_digits) + 'd}' + + # Initialize the list for each row. + rows = [[] for i in range(2*self._num_rings - 1)] + middle = self._num_rings - 1 + + # Start with the degenerate first ring. + universe = universes[-1][0] + rows[middle] = [id_form.format(universe._id)] + + # Add universes one ring at a time. + for r in range(1, self._num_rings): + # r_prime increments down while r increments up. + r_prime = self._num_rings - 1 - r + theta = 0 + y = middle + + # Climb up the top-right. + for i in range(r): + # Add the universe. + universe = universes[r_prime][theta] + rows[y].append(id_form.format(universe._id)) + + # Translate the indices. + y += 1 + theta += 1 + + # Climb left the top-left. + for i in range(r): + # Add the universe. + universe = universes[r_prime][theta] + rows[y].insert(0,id_form.format(universe._id)) + + # Translate the indices. + theta += 1 + + # Climb down the middle left. + for i in range(r): + # Add the universe. + universe = universes[r_prime][theta] + rows[y].insert(0,id_form.format(universe._id)) + + # Translate the indices. + y -= 1 + theta += 1 + + # Climb down the down left. + for i in range(r): + # Add the universe. + universe = universes[r_prime][theta] + rows[y].insert(0, id_form.format(universe._id)) + + # Translate the indices. + y -= 1 + theta += 1 + + # Climb right the down right. + for i in range(r): + # Add the universe. + universe = universes[r_prime][theta] + rows[y].append(id_form.format(universe._id)) + + # Translate the indices. + theta += 1 + + # Climb up the middle-right. + for i in range(r): + # Add the universe. + universe = universes[r_prime][theta] + rows[y].append(id_form.format(universe._id)) + + # Translate the indices. + y += 1 + theta += 1 + + # Flip the rows and join each row into a single string. + rows = [pad.join(x) for x in rows[::-1]] + + # Pad the beginning of the rows so they line up properly. + for y in range(self._num_rings - 1): + rows[y] = (self._num_rings - 1 - y)*pad + rows[y] + rows[-1 - y] = (self._num_rings - 1 - y)*pad + rows[-1 - y] + + # Join the rows together and return the string. + universe_ids = '\n'.join(rows) + return universe_ids + + + def _repr_axial_slice_oy(self, universes): + """Return string representation for the given 2D group of universes in 'OY' orientation case.. + The 'universes' argument should be a list of lists of universes where each sub-list represents a single ring. The first list should be the outer ring. @@ -1527,3 +1793,100 @@ class HexLattice(Lattice): # Join the rows together and return the string. return '\n'.join(rows) + @staticmethod + def show_indices_ox(num_rings):#//DR + """Return a diagram of the hexagonal lattice with OX orientation layout with indices. + + This method can be used to show the proper indices to be used when + setting the :attr:`HexLattice.universes` property. For example, running + this method with num_rings=3 will return the similar diagram:: + + (0, 4) (0, 3) (0, 2) + + (0, 5) (1, 2) (1, 1) (0, 1) + + (0, 6) (1, 3) (2, 0) (1, 0) (0, 0) + + (0, 7) (1, 4) (1, 5) (0,11) + + (0, 8) (0, 9) (0,10) + + Parameters + ---------- + num_rings : int + Number of rings in the hexagonal lattice + + Returns + ------- + str + Diagram of the hexagonal lattice showing indices in OX orientation + + """ + + # Find the largest string and count the number of digits so we can + # properly pad the output string later + largest_index = 6*(num_rings - 1) + n_digits_index = len(str(largest_index)) + n_digits_ring = len(str(num_rings - 1)) + str_form = '({{:{}}},{{:{}}})'.format(n_digits_ring, n_digits_index) + pad = ' '*(n_digits_index + n_digits_ring + 3) + + # Initialize the list for each row. + rows = [[] for i in range(2*num_rings - 1)] + middle = num_rings - 1 + + # Start with the degenerate first ring. + rows[middle] = [str_form.format(num_rings - 1, 0)] + + # Add universes one ring at a time. + for r in range(1, num_rings): + # r_prime increments down while r increments up. + r_prime = num_rings - 1 - r + theta = 0 + y = middle + + for i in range(r): + # Climb up the top-right. + rows[y].append(str_form.format(r_prime, theta)) + y += 1 + theta += 1 + + for i in range(r): + # Climb left the top-left. + rows[y].insert(0,str_form.format(r_prime, theta)) + theta += 1 + + for i in range(r): + # Climb down the middle left. + rows[y].insert(0,str_form.format(r_prime, theta)) + y -= 1 + theta += 1 + + for i in range(r): + # Climb down the down left. + rows[y].insert(0, str_form.format(r_prime, theta)) + y -= 1 + theta += 1 + + for i in range(r): + # Climb right the down right. + rows[y].append(str_form.format(r_prime, theta)) + theta += 1 + + for i in range(r): + # Climb up the middle-right. + rows[y].append(str_form.format(r_prime, theta)) + y += 1 + theta += 1 + + # Flip the rows and join each row into a single string. + rows = [pad.join(x) for x in rows[::-1]] + + # Pad the beginning of the rows so they line up properly. + for y in range(num_rings - 1): + rows[y] = (num_rings - 1 - y)*pad + rows[y] + rows[-1 - y] = (num_rings - 1 - y)*pad + rows[-1 - y] + + # Join the rows together and return the string. + return '\n\n'.join(rows) + diff --git a/src/lattice.cpp b/src/lattice.cpp index 6e0f2331f..49161195d 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -429,6 +429,8 @@ HexLattice::HexLattice(pugi::xml_node lat_node) { type_ = LatticeType::hex; + hextype = 0;//DR a default OY orientation of lattice + // Read the number of lattice cells in each dimension. n_rings_ = std::stoi(get_node_value(lat_node, "n_rings")); if (check_for_node(lat_node, "n_axial")) { @@ -439,6 +441,12 @@ HexLattice::HexLattice(pugi::xml_node lat_node) is_3d_ = false; } + //DR Check if "OX" orientation desc in node + if (check_for_node(lat_node, "orient")) { + std :: string orientation = get_node_value(lat_node, "orient"); + if (orientation == "OX") {hextype=1;} + } + // Read the lattice center. std::string center_str {get_node_value(lat_node, "center")}; std::vector center_words {split(center_str)}; @@ -481,96 +489,173 @@ HexLattice::HexLattice(pugi::xml_node lat_node) // Parse the universes. // Universes in hexagonal lattices are stored in a manner that represents - // a skewed coordinate system: (x, alpha) rather than (x, y). There is + // a skewed coordinate system: (x, alpha) in case of OY orientation + // and (alpha,y) in OX one rather than (x, y). There is // no obvious, direct relationship between the order of universes in the // input and the order that they will be stored in the skewed array so // the following code walks a set of index values across the skewed array // in a manner that matches the input order. Note that i_x = 0, i_a = 0 - // corresponds to the center of the hexagonal lattice. + // or i_a=0 , i_y = 0 corresponds to the center of the hexagonal lattice. universes_.resize((2*n_rings_-1) * (2*n_rings_-1) * n_axial_, C_NONE); - int input_index = 0; - for (int m = 0; m < n_axial_; m++) { - // Initialize lattice indecies. - int i_x = 1; - int i_a = n_rings_ - 1; - // Map upper triangular region of hexagonal lattice which is found in the - // first n_rings-1 rows of the input. - for (int k = 0; k < n_rings_-1; k++) { - // Walk the index to lower-left neighbor of last row start. - i_x -= 1; - - // Iterate over the input columns. - for (int j = 0; j < k+1; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_a+n_rings_-1) - + (i_x+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - // Walk the index to the right neighbor (which is not adjacent). - i_x += 2; - i_a -= 1; - } - - // Return the lattice index to the start of the current row. - i_x -= 2 * (k+1); - i_a += (k+1); - } - - // Map the middle square region of the hexagonal lattice which is found in - // the next 2*n_rings-1 rows of the input. - for (int k = 0; k < 2*n_rings_-1; k++) { - if ((k % 2) == 0) { - // Walk the index to the lower-left neighbor of the last row start. - i_x -= 1; - } else { - // Walk the index to the lower-right neighbor of the last row start. - i_x += 1; - i_a -= 1; - } - - // Iterate over the input columns. - for (int j = 0; j < n_rings_ - (k % 2); j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_a+n_rings_-1) - + (i_x+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - // Walk the index to the right neighbor (which is not adjacent). - i_x += 2; - i_a -= 1; - } - - // Return the lattice index to the start of the current row. - i_x -= 2*(n_rings_ - (k % 2)); - i_a += n_rings_ - (k % 2); - } - - // Map the lower triangular region of the hexagonal lattice. - for (int k = 0; k < n_rings_-1; k++) { - // Walk the index to the lower-right neighbor of the last row start. - i_x += 1; - i_a -= 1; - - // Iterate over the input columns. - for (int j = 0; j < n_rings_-k-1; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_a+n_rings_-1) - + (i_x+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - // Walk the index to the right neighbor (which is not adjacent). - i_x += 2; - i_a -= 1; - } - - // Return lattice index to start of current row. - i_x -= 2*(n_rings_ - k - 1); - i_a += n_rings_ - k - 1; - } + if (hextype==0) { + fill_lattice_oy(univ_words); + } else + { + fill_lattice_ox(univ_words); } + + } +// DR fill universes in OX orientation +//============================================================================== +void +HexLattice::fill_lattice_ox(std::vector univ_words){ + + int input_index = 0; + for (int m = 0; m < n_axial_; m++) { + // Initialize lattice indecies. + int i_a = -(n_rings_ - 1); + int i_y = n_rings_ - 1; + + + //DR Map upper region of hexagonal lattice which is found in the + // first n_rings-1 rows of the input. + for (int k = 0; k < n_rings_-1; k++) { + + // Iterate over the input columns. + for (int j = 0; j < k+n_rings_; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_y+n_rings_-1) + + (i_a+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + //DR Move to the next right neighbour cell + i_a += 1; + + } + + // Return the lattice index to the start of the current row. + i_a = -(n_rings_ - 1); + i_y -= 1; + } + + // Map the lower region from the centerline of cart to down side + for (int k = 0; k < n_rings_; k++) { + // Walk the index to the lower-right neighbor of the last row start. + i_a = -(n_rings_ - 1) + k; + + // Iterate over the input columns. + for (int j = 0; j < 2*n_rings_-k-1; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_y+n_rings_-1) + + (i_a+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + //DR Move to the next right neighbour cell + i_a += 1; + + } + + // Return lattice index to start of current row. + i_y -= 1; + } + } +} + +//============================================================================== + + +// DR fill universes in OY orientation +//============================================================================== +void +HexLattice::fill_lattice_oy(std::vector univ_words){ + + int input_index = 0; + for (int m = 0; m < n_axial_; m++) { + // Initialize lattice indecies. + int i_x = 1; + int i_a = n_rings_ - 1; + + // Map upper triangular region of hexagonal lattice which is found in the + // first n_rings-1 rows of the input. + for (int k = 0; k < n_rings_-1; k++) { + // Walk the index to lower-left neighbor of last row start. + i_x -= 1; + + // Iterate over the input columns. + for (int j = 0; j < k+1; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_a+n_rings_-1) + + (i_x+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return the lattice index to the start of the current row. + i_x -= 2 * (k+1); + i_a += (k+1); + } + + // Map the middle square region of the hexagonal lattice which is found in + // the next 2*n_rings-1 rows of the input. + for (int k = 0; k < 2*n_rings_-1; k++) { + if ((k % 2) == 0) { + // Walk the index to the lower-left neighbor of the last row start. + i_x -= 1; + } else { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; + } + + // Iterate over the input columns. + for (int j = 0; j < n_rings_ - (k % 2); j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_a+n_rings_-1) + + (i_x+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return the lattice index to the start of the current row. + i_x -= 2*(n_rings_ - (k % 2)); + i_a += n_rings_ - (k % 2); + } + + // Map the lower triangular region of the hexagonal lattice. + for (int k = 0; k < n_rings_-1; k++) { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; + + // Iterate over the input columns. + for (int j = 0; j < n_rings_-k-1; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_a+n_rings_-1) + + (i_x+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return lattice index to start of current row. + i_x -= 2*(n_rings_ - k - 1); + i_a += n_rings_ - k - 1; + } + } +} + +//============================================================================== //============================================================================== @@ -603,11 +688,127 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const && (i_xyz[2] < n_axial_)); } -//============================================================================== - +//DR wrapper HexLattice function both for OX and OY orientation std::pair> HexLattice::distance(Position r, Direction u, const std::array& i_xyz) const +{ + if (hextype==0){ + return distance_oy(r,u,i_xyz); + } else { + return distance_ox(r,u,i_xyz); + } +} + +//DR OX distance calculation +std::pair> +HexLattice::distance_ox(Position r, Direction u, const std::array& i_xyz) +const +{ + // Compute the direction on the hexagonal basis. + // DR for OX angle=60 g + double beta_dir = u.x*0.5 + u.y * std::sqrt(3.0) / 2.0; + double gamma_dir = u.x*0.5 - u.y * std::sqrt(3.0) / 2.0; + + // Note that hexagonal lattice distance calculations are performed + // using the particle's coordinates relative to the neighbor lattice + // cells, not relative to the particle's current cell. This is done + // because there is significant disagreement between neighboring cells + // on where the lattice boundary is due to finite precision issues. + + // Upper-right and lower-left sides. + double d {INFTY}; + std::array lattice_trans; + double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge + Position r_t; + if (beta_dir > 0) { + const std::array i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double beta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0; + if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { + d = (edge - beta) / beta_dir; + if (beta_dir > 0) { + lattice_trans = {0, 1, 0}; + } else { + lattice_trans = {0, -1, 0}; + } + } + + // Lower-right and upper-left sides. + edge = -copysign(0.5*pitch_[0], gamma_dir); + if (gamma_dir > 0) { + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0; + if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { + double this_d = (edge - gamma) / gamma_dir; + if (this_d < d) { + if (gamma_dir > 0) { + lattice_trans = {1, -1, 0}; + } else { + lattice_trans = {-1, 1, 0}; + } + d = this_d; + } + } + + // right and left sides. + edge = -copysign(0.5*pitch_[0], u.x); + if (u.x > 0) { + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + if ((std::abs(r_t.x - edge) > FP_PRECISION) && u.x != 0) { + double this_d = (edge - r_t.x) / u.x; + if (this_d < d) { + if (u.x > 0) { + lattice_trans = {1, 0, 0}; + } else { + lattice_trans = {-1, 0, 0}; + } + d = this_d; + } + } + + // Top and bottom sides + if (is_3d_) { + double z = r.z; + double z0 {copysign(0.5 * pitch_[1], u.z)}; + if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { + double this_d = (z0 - z) / u.z; + if (this_d < d) { + d = this_d; + if (u.z > 0) { + lattice_trans = {0, 0, 1}; + } else { + lattice_trans = {0, 0, -1}; + } + d = this_d; + } + } + } + + return {d, lattice_trans}; +} + +//============================================================================== + +//============================================================================== +//DR OY distance calculation +std::pair> +HexLattice::distance_oy(Position r, Direction u, const std::array& i_xyz) +const { // Compute the direction on the hexagonal basis. double beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; @@ -705,6 +906,7 @@ const return {d, lattice_trans}; } +//============================================================================== //============================================================================== std::array @@ -714,7 +916,8 @@ HexLattice::get_indices(Position r, Direction u) const // information. As a result, we move the position slightly forward to // determine what lattice index the particle is most likely to be in. r += TINY_BIT * u; - + float dbg1; + float dbg2; // Offset the xyz by the lattice center. Position r_o {r.x - center_.x, r.y - center_.y, r.z}; if (is_3d_) {r_o.z -= center_.z;} @@ -727,11 +930,21 @@ HexLattice::get_indices(Position r, Direction u) const out[2] = 0; } + if (hextype==0){//DR type implement OY default behaviour // Convert coordinates into skewed bases. The (x, alpha) basis is used to // find the index of the global coordinates to within 4 cells. double alpha = r_o.y - r_o.x / std::sqrt(3.0); out[0] = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0])); - out[1] = std::floor(alpha / pitch_[0]); + out[1] = std::floor(alpha / pitch_[0]);} + else{//DR type implement OX (alpha,y) basis + // Convert coordinates into skewed bases. The (alpha,y) basis is used to + // find the index of the global coordinates to within 4 cells. + double alpha = r_o.y - r_o.x * std::sqrt(3.0); + dbg1=-alpha / (std::sqrt(3.0) * pitch_[0]); + dbg2=r_o.y / 0.5*std::sqrt(3.0) * pitch_[0]; + out[0] = std::floor(-alpha / (std::sqrt(3.0) * pitch_[0])); + out[1] = std::floor(r_o.y / (0.5*std::sqrt(3.0) * pitch_[0])); + } // Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but // the array is offset so that the indices never go below 0). @@ -781,16 +994,32 @@ Position HexLattice::get_local_position(Position r, const std::array i_xyz) const { - // x_l = x_g - (center + pitch_x*cos(30)*index_x) - r.x -= center_.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0]; - // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) - r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0] - + (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0); - if (is_3d_) { - r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; - } + // DR for OY implementation + if (hextype == 0){ + // x_l = x_g - (center + pitch_x*cos(30)*index_x) + r.x -= center_.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0]; + // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) + r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0] + + (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0); + if (is_3d_) { + r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; + } - return r; + } + else {// DR for OX implementation + + // x_l = x_g - (center + pitch_x*index_a + pitch_y*sin(30)*index_y) + r.x -= (center_.x + (i_xyz[0] - n_rings_ + 1) * pitch_[0] + + (i_xyz[1] - n_rings_ + 1) * pitch_[0] / 2.0); + // y_l = y_g - (center + pitch_y*cos(30)*index_y) + r.y -= center_.y + std::sqrt(3.0)/2.0 * (i_xyz[1] - n_rings_ + 1) * pitch_[0]; + + if (is_3d_) { + r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; + } + + } + return r; } //============================================================================== @@ -848,6 +1077,12 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const write_string(lat_group, "type", "hexagonal", false); write_dataset(lat_group, "n_rings", n_rings_); write_dataset(lat_group, "n_axial", n_axial_); + if (hextype == 0){ + write_string(lat_group, "orientation", "oy",false); + } + else { + write_string(lat_group, "orientation", "ox",false); + } if (is_3d_) { write_dataset(lat_group, "pitch", pitch_); write_dataset(lat_group, "center", center_); diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index c9c53f6d6..30fdf2c8a 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -36,6 +36,7 @@ element geometry { attribute dimension { list { xsd:positiveInteger+ } }) & (element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) & (element pitch { list { xsd:double+ } } | attribute pitch { list { xsd:double+ } }) & + (element orient { ( "OX" ) } | attribute orient { ( "OX" ) })? & (element universes { list { xsd:int+ } } | attribute universes { list { xsd:int+ } }) & (element outer { xsd:int } | attribute outer { xsd:int })? }* diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng index d2c1f635e..e6c5f6ca1 100644 --- a/src/relaxng/geometry.rng +++ b/src/relaxng/geometry.rng @@ -294,6 +294,16 @@ + + + + OX + + + OX + + + From 4ef5ad64e31631869158b369b8c4ae727e6a7139 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Wed, 10 Apr 2019 10:58:35 +0300 Subject: [PATCH 04/31] garbage remove from commit --- .cproject | 108 --------------------- .project | 27 ------ .settings/org.eclipse.cdt.codan.core.prefs | 75 -------------- 3 files changed, 210 deletions(-) delete mode 100644 .cproject delete mode 100644 .project delete mode 100644 .settings/org.eclipse.cdt.codan.core.prefs diff --git a/.cproject b/.cproject deleted file mode 100644 index bdc88d84b..000000000 --- a/.cproject +++ /dev/null @@ -1,108 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/.project b/.project deleted file mode 100644 index 5eb05adae..000000000 --- a/.project +++ /dev/null @@ -1,27 +0,0 @@ - - - openmccc - - - - - - org.eclipse.cdt.managedbuilder.core.genmakebuilder - clean,full,incremental, - - - - - org.eclipse.cdt.managedbuilder.core.ScannerConfigBuilder - full,incremental, - - - - - - org.eclipse.cdt.core.cnature - org.eclipse.cdt.core.ccnature - org.eclipse.cdt.managedbuilder.core.managedBuildNature - org.eclipse.cdt.managedbuilder.core.ScannerConfigNature - - diff --git a/.settings/org.eclipse.cdt.codan.core.prefs b/.settings/org.eclipse.cdt.codan.core.prefs deleted file mode 100644 index 77496f558..000000000 --- a/.settings/org.eclipse.cdt.codan.core.prefs +++ /dev/null @@ -1,75 +0,0 @@ -eclipse.preferences.version=1 -org.eclipse.cdt.codan.checkers.errnoreturn=Warning -org.eclipse.cdt.codan.checkers.errnoreturn.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"No return\\")",implicit\=>false} -org.eclipse.cdt.codan.checkers.errreturnvalue=Error -org.eclipse.cdt.codan.checkers.errreturnvalue.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Unused return value\\")"} -org.eclipse.cdt.codan.checkers.nocommentinside=-Error -org.eclipse.cdt.codan.checkers.nocommentinside.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Nesting comments\\")"} -org.eclipse.cdt.codan.checkers.nolinecomment=-Error -org.eclipse.cdt.codan.checkers.nolinecomment.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Line comments\\")"} -org.eclipse.cdt.codan.checkers.noreturn=Error -org.eclipse.cdt.codan.checkers.noreturn.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"No return value\\")",implicit\=>false} -org.eclipse.cdt.codan.internal.checkers.AbstractClassCreation=Error -org.eclipse.cdt.codan.internal.checkers.AbstractClassCreation.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Abstract class cannot be instantiated\\")"} -org.eclipse.cdt.codan.internal.checkers.AmbiguousProblem=Error -org.eclipse.cdt.codan.internal.checkers.AmbiguousProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Ambiguous problem\\")"} -org.eclipse.cdt.codan.internal.checkers.AssignmentInConditionProblem=Warning -org.eclipse.cdt.codan.internal.checkers.AssignmentInConditionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Assignment in condition\\")"} -org.eclipse.cdt.codan.internal.checkers.AssignmentToItselfProblem=Error -org.eclipse.cdt.codan.internal.checkers.AssignmentToItselfProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Assignment to itself\\")"} -org.eclipse.cdt.codan.internal.checkers.CaseBreakProblem=Warning -org.eclipse.cdt.codan.internal.checkers.CaseBreakProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"No break at end of case\\")",no_break_comment\=>"no break",last_case_param\=>false,empty_case_param\=>false,enable_fallthrough_quickfix_param\=>false} -org.eclipse.cdt.codan.internal.checkers.CatchByReference=Warning -org.eclipse.cdt.codan.internal.checkers.CatchByReference.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Catching by reference is recommended\\")",unknown\=>false,exceptions\=>()} -org.eclipse.cdt.codan.internal.checkers.CircularReferenceProblem=Error -org.eclipse.cdt.codan.internal.checkers.CircularReferenceProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Circular inheritance\\")"} -org.eclipse.cdt.codan.internal.checkers.ClassMembersInitialization=Warning -org.eclipse.cdt.codan.internal.checkers.ClassMembersInitialization.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Class members should be properly initialized\\")",skip\=>true} -org.eclipse.cdt.codan.internal.checkers.DecltypeAutoProblem=Error -org.eclipse.cdt.codan.internal.checkers.DecltypeAutoProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid 'decltype(auto)' specifier\\")"} -org.eclipse.cdt.codan.internal.checkers.FieldResolutionProblem=Error -org.eclipse.cdt.codan.internal.checkers.FieldResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Field cannot be resolved\\")"} -org.eclipse.cdt.codan.internal.checkers.FunctionResolutionProblem=Error -org.eclipse.cdt.codan.internal.checkers.FunctionResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Function cannot be resolved\\")"} -org.eclipse.cdt.codan.internal.checkers.InvalidArguments=Error -org.eclipse.cdt.codan.internal.checkers.InvalidArguments.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid arguments\\")"} -org.eclipse.cdt.codan.internal.checkers.InvalidTemplateArgumentsProblem=Error -org.eclipse.cdt.codan.internal.checkers.InvalidTemplateArgumentsProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid template argument\\")"} -org.eclipse.cdt.codan.internal.checkers.LabelStatementNotFoundProblem=Error -org.eclipse.cdt.codan.internal.checkers.LabelStatementNotFoundProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Label statement not found\\")"} -org.eclipse.cdt.codan.internal.checkers.MemberDeclarationNotFoundProblem=Error -org.eclipse.cdt.codan.internal.checkers.MemberDeclarationNotFoundProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Member declaration not found\\")"} -org.eclipse.cdt.codan.internal.checkers.MethodResolutionProblem=Error -org.eclipse.cdt.codan.internal.checkers.MethodResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Method cannot be resolved\\")"} -org.eclipse.cdt.codan.internal.checkers.NamingConventionFunctionChecker=-Info -org.eclipse.cdt.codan.internal.checkers.NamingConventionFunctionChecker.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Name convention for function\\")",pattern\=>"^[a-z]",macro\=>true,exceptions\=>()} -org.eclipse.cdt.codan.internal.checkers.NonVirtualDestructorProblem=Warning -org.eclipse.cdt.codan.internal.checkers.NonVirtualDestructorProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Class has a virtual method and non-virtual destructor\\")"} -org.eclipse.cdt.codan.internal.checkers.OverloadProblem=Error -org.eclipse.cdt.codan.internal.checkers.OverloadProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid overload\\")"} -org.eclipse.cdt.codan.internal.checkers.RedeclarationProblem=Error -org.eclipse.cdt.codan.internal.checkers.RedeclarationProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid redeclaration\\")"} -org.eclipse.cdt.codan.internal.checkers.RedefinitionProblem=Error -org.eclipse.cdt.codan.internal.checkers.RedefinitionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Invalid redefinition\\")"} -org.eclipse.cdt.codan.internal.checkers.ReturnStyleProblem=-Warning -org.eclipse.cdt.codan.internal.checkers.ReturnStyleProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Return with parenthesis\\")"} -org.eclipse.cdt.codan.internal.checkers.ScanfFormatStringSecurityProblem=-Warning -org.eclipse.cdt.codan.internal.checkers.ScanfFormatStringSecurityProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Format String Vulnerability\\")"} -org.eclipse.cdt.codan.internal.checkers.StatementHasNoEffectProblem=Warning -org.eclipse.cdt.codan.internal.checkers.StatementHasNoEffectProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Statement has no effect\\")",macro\=>true,exceptions\=>()} -org.eclipse.cdt.codan.internal.checkers.SuggestedParenthesisProblem=Warning -org.eclipse.cdt.codan.internal.checkers.SuggestedParenthesisProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Suggested parenthesis around expression\\")",paramNot\=>false} -org.eclipse.cdt.codan.internal.checkers.SuspiciousSemicolonProblem=Warning -org.eclipse.cdt.codan.internal.checkers.SuspiciousSemicolonProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Suspicious semicolon\\")",else\=>false,afterelse\=>false} -org.eclipse.cdt.codan.internal.checkers.TypeResolutionProblem=Error -org.eclipse.cdt.codan.internal.checkers.TypeResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Type cannot be resolved\\")"} -org.eclipse.cdt.codan.internal.checkers.UnusedFunctionDeclarationProblem=Warning -org.eclipse.cdt.codan.internal.checkers.UnusedFunctionDeclarationProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Unused function declaration\\")",macro\=>true} -org.eclipse.cdt.codan.internal.checkers.UnusedStaticFunctionProblem=Warning -org.eclipse.cdt.codan.internal.checkers.UnusedStaticFunctionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Unused static function\\")",macro\=>true} -org.eclipse.cdt.codan.internal.checkers.UnusedVariableDeclarationProblem=Warning -org.eclipse.cdt.codan.internal.checkers.UnusedVariableDeclarationProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Unused variable declaration in file scope\\")",macro\=>true,exceptions\=>("@(\#)","$Id")} -org.eclipse.cdt.codan.internal.checkers.VariableResolutionProblem=Error -org.eclipse.cdt.codan.internal.checkers.VariableResolutionProblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>true,RUN_ON_INC_BUILD\=>true,RUN_ON_FILE_OPEN\=>false,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>"@suppress(\\"Symbol is not resolved\\")"} -org.eclipse.cdt.qt.core.qtproblem=Warning -org.eclipse.cdt.qt.core.qtproblem.params={launchModes\=>{RUN_ON_FULL_BUILD\=>false,RUN_ON_INC_BUILD\=>false,RUN_ON_FILE_OPEN\=>true,RUN_ON_FILE_SAVE\=>false,RUN_AS_YOU_TYPE\=>true,RUN_ON_DEMAND\=>true},suppression_comment\=>null} From 7813b284ad5a58add2f2bc9a9f27cb4385f9ee1d Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Thu, 11 Apr 2019 12:23:44 +0300 Subject: [PATCH 05/31] code was corrected,tabs were replaced by spaces, results for tallies test get from current version --- include/openmc/lattice.h | 2 +- openmc/lattice.py | 27 +- src/lattice.cpp | 290 +++++++++--------- .../regression_tests/tallies/results_true.dat | 2 +- 4 files changed, 164 insertions(+), 157 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index 752e928c4..7c06117b5 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -278,7 +278,7 @@ public: private: int n_rings_; //!< Number of radial tile positions int n_axial_; //!< Number of axial tile positions - int hextype; //!DR type of hexagonal lattice orientation //0 - OY by default;1-OX + int hextype; //!DR type of hexagonal lattice orientation //0 - OY by default;1-OX Position center_; //!< Global center of lattice std::array pitch_; //!< Lattice tile width and height }; diff --git a/openmc/lattice.py b/openmc/lattice.py index 699a7540d..1f0ec5466 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -141,7 +141,11 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): center = group['center'][...] pitch = group['pitch'][...] outer = group['outer'].value - orientation = group['orientation'].value + #DR Temp added for compatibility with previous version + if ('orientation' in group.keys()): + orientation = group['orientation'].value + else: + orientation = "oy" universe_ids = group['universes'][...] @@ -1133,14 +1137,16 @@ class HexLattice(Lattice): else: z = point[2] - self.center[2] iz = floor(z/self.pitch[1] + 0.5*self.num_axial) - if self._hextype is 0:#//DR + if self._hextype:#//DR + alpha = y - x*sqrt(3.) + ix = 0 + ia = floor(-alpha/(sqrt(3.0) * self.pitch[0])) + iy = floor(y/(sqrt(0.75) * self.pitch[0])) + else: alpha = y - x/sqrt(3.) ix = floor(x/(sqrt(0.75) * self.pitch[0])) ia = floor(alpha/self.pitch[0]) - else: - alpha = y - x*sqrt(3.) - ia = floor(-alpha/(sqrt(3.0) * self.pitch[0])) - iy = floor(y/(sqrt(0.75) * self.pitch[0])) + iy = 0 # Check four lattice elements to see which one is closest based on local @@ -1169,7 +1175,8 @@ class HexLattice(Lattice): point : Iterable of float Cartesian coordinates of point idx : Iterable of int - Indices of lattice element in :math:`(x,\alpha,z)` bases + Indices of lattice element in :math:`(x,\alpha,z)` + or :math:`(\alpha,y,z)` basesis Returns ------- @@ -1179,11 +1186,11 @@ class HexLattice(Lattice): """ if self._hextype:#//DR - x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0]) - y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0]) - else: x = point[0] - (self.center[0] + self.pitch[0]*idx[0] + 0.5*self.pitch[0]*idx[1]) y = point[1] - (self.center[1] + (sqrt(0.75)*self.pitch[0]*idx[1])) + else: + x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0]) + y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0]) if self._num_axial is None: z = point[2] diff --git a/src/lattice.cpp b/src/lattice.cpp index 49161195d..028932010 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -443,8 +443,8 @@ HexLattice::HexLattice(pugi::xml_node lat_node) //DR Check if "OX" orientation desc in node if (check_for_node(lat_node, "orient")) { - std :: string orientation = get_node_value(lat_node, "orient"); - if (orientation == "OX") {hextype=1;} + std :: string orientation = get_node_value(lat_node, "orient"); + if (orientation == "OX") {hextype=1;} } // Read the lattice center. @@ -500,10 +500,10 @@ HexLattice::HexLattice(pugi::xml_node lat_node) universes_.resize((2*n_rings_-1) * (2*n_rings_-1) * n_axial_, C_NONE); if (hextype==0) { - fill_lattice_oy(univ_words); + fill_lattice_oy(univ_words); } else { - fill_lattice_ox(univ_words); + fill_lattice_ox(univ_words); } @@ -513,55 +513,55 @@ HexLattice::HexLattice(pugi::xml_node lat_node) void HexLattice::fill_lattice_ox(std::vector univ_words){ - int input_index = 0; - for (int m = 0; m < n_axial_; m++) { - // Initialize lattice indecies. - int i_a = -(n_rings_ - 1); - int i_y = n_rings_ - 1; + int input_index = 0; + for (int m = 0; m < n_axial_; m++) { + // Initialize lattice indecies. + int i_a = -(n_rings_ - 1); + int i_y = n_rings_ - 1; - //DR Map upper region of hexagonal lattice which is found in the - // first n_rings-1 rows of the input. - for (int k = 0; k < n_rings_-1; k++) { + //DR Map upper region of hexagonal lattice which is found in the + // first n_rings-1 rows of the input. + for (int k = 0; k < n_rings_-1; k++) { - // Iterate over the input columns. - for (int j = 0; j < k+n_rings_; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_y+n_rings_-1) - + (i_a+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - //DR Move to the next right neighbour cell - i_a += 1; + // Iterate over the input columns. + for (int j = 0; j < k+n_rings_; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_y+n_rings_-1) + + (i_a+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + //DR Move to the next right neighbour cell + i_a += 1; - } + } - // Return the lattice index to the start of the current row. - i_a = -(n_rings_ - 1); - i_y -= 1; - } + // Return the lattice index to the start of the current row. + i_a = -(n_rings_ - 1); + i_y -= 1; + } - // Map the lower region from the centerline of cart to down side - for (int k = 0; k < n_rings_; k++) { - // Walk the index to the lower-right neighbor of the last row start. - i_a = -(n_rings_ - 1) + k; + // Map the lower region from the centerline of cart to down side + for (int k = 0; k < n_rings_; k++) { + // Walk the index to the lower-right neighbor of the last row start. + i_a = -(n_rings_ - 1) + k; - // Iterate over the input columns. - for (int j = 0; j < 2*n_rings_-k-1; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_y+n_rings_-1) - + (i_a+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - //DR Move to the next right neighbour cell - i_a += 1; + // Iterate over the input columns. + for (int j = 0; j < 2*n_rings_-k-1; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_y+n_rings_-1) + + (i_a+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + //DR Move to the next right neighbour cell + i_a += 1; - } + } - // Return lattice index to start of current row. - i_y -= 1; - } - } + // Return lattice index to start of current row. + i_y -= 1; + } + } } //============================================================================== @@ -572,87 +572,87 @@ HexLattice::fill_lattice_ox(std::vector univ_words){ void HexLattice::fill_lattice_oy(std::vector univ_words){ - int input_index = 0; - for (int m = 0; m < n_axial_; m++) { - // Initialize lattice indecies. - int i_x = 1; - int i_a = n_rings_ - 1; + int input_index = 0; + for (int m = 0; m < n_axial_; m++) { + // Initialize lattice indecies. + int i_x = 1; + int i_a = n_rings_ - 1; - // Map upper triangular region of hexagonal lattice which is found in the - // first n_rings-1 rows of the input. - for (int k = 0; k < n_rings_-1; k++) { - // Walk the index to lower-left neighbor of last row start. - i_x -= 1; + // Map upper triangular region of hexagonal lattice which is found in the + // first n_rings-1 rows of the input. + for (int k = 0; k < n_rings_-1; k++) { + // Walk the index to lower-left neighbor of last row start. + i_x -= 1; - // Iterate over the input columns. - for (int j = 0; j < k+1; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_a+n_rings_-1) - + (i_x+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - // Walk the index to the right neighbor (which is not adjacent). - i_x += 2; - i_a -= 1; - } + // Iterate over the input columns. + for (int j = 0; j < k+1; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_a+n_rings_-1) + + (i_x+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } - // Return the lattice index to the start of the current row. - i_x -= 2 * (k+1); - i_a += (k+1); - } + // Return the lattice index to the start of the current row. + i_x -= 2 * (k+1); + i_a += (k+1); + } - // Map the middle square region of the hexagonal lattice which is found in - // the next 2*n_rings-1 rows of the input. - for (int k = 0; k < 2*n_rings_-1; k++) { - if ((k % 2) == 0) { - // Walk the index to the lower-left neighbor of the last row start. - i_x -= 1; - } else { - // Walk the index to the lower-right neighbor of the last row start. - i_x += 1; - i_a -= 1; - } + // Map the middle square region of the hexagonal lattice which is found in + // the next 2*n_rings-1 rows of the input. + for (int k = 0; k < 2*n_rings_-1; k++) { + if ((k % 2) == 0) { + // Walk the index to the lower-left neighbor of the last row start. + i_x -= 1; + } else { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; + } - // Iterate over the input columns. - for (int j = 0; j < n_rings_ - (k % 2); j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_a+n_rings_-1) - + (i_x+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - // Walk the index to the right neighbor (which is not adjacent). - i_x += 2; - i_a -= 1; - } + // Iterate over the input columns. + for (int j = 0; j < n_rings_ - (k % 2); j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_a+n_rings_-1) + + (i_x+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } - // Return the lattice index to the start of the current row. - i_x -= 2*(n_rings_ - (k % 2)); - i_a += n_rings_ - (k % 2); - } + // Return the lattice index to the start of the current row. + i_x -= 2*(n_rings_ - (k % 2)); + i_a += n_rings_ - (k % 2); + } - // Map the lower triangular region of the hexagonal lattice. - for (int k = 0; k < n_rings_-1; k++) { - // Walk the index to the lower-right neighbor of the last row start. - i_x += 1; - i_a -= 1; + // Map the lower triangular region of the hexagonal lattice. + for (int k = 0; k < n_rings_-1; k++) { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; - // Iterate over the input columns. - for (int j = 0; j < n_rings_-k-1; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_a+n_rings_-1) - + (i_x+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - // Walk the index to the right neighbor (which is not adjacent). - i_x += 2; - i_a -= 1; - } + // Iterate over the input columns. + for (int j = 0; j < n_rings_-k-1; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_a+n_rings_-1) + + (i_x+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } - // Return lattice index to start of current row. - i_x -= 2*(n_rings_ - k - 1); - i_a += n_rings_ - k - 1; - } - } + // Return lattice index to start of current row. + i_x -= 2*(n_rings_ - k - 1); + i_a += n_rings_ - k - 1; + } + } } //============================================================================== @@ -693,11 +693,11 @@ std::pair> HexLattice::distance(Position r, Direction u, const std::array& i_xyz) const { - if (hextype==0){ - return distance_oy(r,u,i_xyz); - } else { - return distance_ox(r,u,i_xyz); - } + if (hextype==0){ + return distance_oy(r,u,i_xyz); + } else { + return distance_ox(r,u,i_xyz); + } } //DR OX distance calculation @@ -994,32 +994,32 @@ Position HexLattice::get_local_position(Position r, const std::array i_xyz) const { - // DR for OY implementation - if (hextype == 0){ - // x_l = x_g - (center + pitch_x*cos(30)*index_x) - r.x -= center_.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0]; - // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) - r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0] - + (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0); - if (is_3d_) { - r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; - } + // DR for OY implementation + if (hextype == 0){ + // x_l = x_g - (center + pitch_x*cos(30)*index_x) + r.x -= center_.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0]; + // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) + r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0] + + (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0); + if (is_3d_) { + r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; + } - } - else {// DR for OX implementation + } + else {// DR for OX implementation - // x_l = x_g - (center + pitch_x*index_a + pitch_y*sin(30)*index_y) - r.x -= (center_.x + (i_xyz[0] - n_rings_ + 1) * pitch_[0] - + (i_xyz[1] - n_rings_ + 1) * pitch_[0] / 2.0); - // y_l = y_g - (center + pitch_y*cos(30)*index_y) - r.y -= center_.y + std::sqrt(3.0)/2.0 * (i_xyz[1] - n_rings_ + 1) * pitch_[0]; + // x_l = x_g - (center + pitch_x*index_a + pitch_y*sin(30)*index_y) + r.x -= (center_.x + (i_xyz[0] - n_rings_ + 1) * pitch_[0] + + (i_xyz[1] - n_rings_ + 1) * pitch_[0] / 2.0); + // y_l = y_g - (center + pitch_y*cos(30)*index_y) + r.y -= center_.y + std::sqrt(3.0)/2.0 * (i_xyz[1] - n_rings_ + 1) * pitch_[0]; - if (is_3d_) { - r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; - } + if (is_3d_) { + r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; + } - } - return r; + } + return r; } //============================================================================== @@ -1078,10 +1078,10 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const write_dataset(lat_group, "n_rings", n_rings_); write_dataset(lat_group, "n_axial", n_axial_); if (hextype == 0){ - write_string(lat_group, "orientation", "oy",false); + write_string(lat_group, "orientation", "oy",false); } else { - write_string(lat_group, "orientation", "ox",false); + write_string(lat_group, "orientation", "ox",false); } if (is_3d_) { write_dataset(lat_group, "pitch", pitch_); diff --git a/tests/regression_tests/tallies/results_true.dat b/tests/regression_tests/tallies/results_true.dat index 9830ff872..de1a18eff 100644 --- a/tests/regression_tests/tallies/results_true.dat +++ b/tests/regression_tests/tallies/results_true.dat @@ -1 +1 @@ -de773a799f84348241ebd65bf7beae8114861d8674b652bfd443d578c146a7d37ca38f2efc5d05124fdbb1cf12829907768351e6b2d6b5f4bd2c121417757507 \ No newline at end of file +c37e0468f1684f7810429332721762884f3cdf0429d38caa99500bbde04ad84cf1e1639adc46b6db44b28b6f31028ed029f85920c8164d77463b6d89336043a8 \ No newline at end of file From 7059b56e79095d43d707cc7ce2c7fe04d3175608 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Thu, 23 May 2019 12:29:19 +0300 Subject: [PATCH 06/31] one distance function --- include/openmc/lattice.h | 8 +- src/lattice.cpp | 319 +++++++++++++++------------------------ 2 files changed, 119 insertions(+), 208 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index e2547afd0..afec39b22 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -249,12 +249,6 @@ public: bool are_valid_indices(const int i_xyz[3]) const; std::pair> distance(Position r, Direction u, const std::array& i_xyz) const; - //DR OX distance calculation - std::pair> - distance_ox(Position r, Direction u, const std::array& i_xyz) const; - //DR OY distance calculation - std::pair> - distance_oy(Position r, Direction u, const std::array& i_xyz) const; std::array get_indices(Position r, Direction u) const; @@ -271,7 +265,7 @@ public: int32_t& offset(int map, const int i_xyz[3]); std::string index_to_string(int indx) const; - + //TODO: find a desicion for compatible old hdf5 file without orientation void to_hdf5_inner(hid_t group_id) const; private: diff --git a/src/lattice.cpp b/src/lattice.cpp index 820d4594e..af025033a 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -694,217 +694,134 @@ std::pair> HexLattice::distance(Position r, Direction u, const std::array& i_xyz) const { + double cosa,sina,u_xy; + // Index + 1 of translation direction + int i_trans {1}; + xt :: xarray translation_matrix; if (hextype==0){ - return distance_oy(r,u,i_xyz); + cosa = std::sqrt(3.0) / 2.0; + sina = 0.5; + u_xy = u.y; + //Translation matrix 4x3 {{beta},{gamma},{oy},{oz}} for OY case + translation_matrix = {{1, 0, 0}, + {1, -1, 0}, + {0, 1, 0}, + {0, 0, 1}}; } else { - return distance_ox(r,u,i_xyz); + cosa = 0.5; + sina = std::sqrt(3.0) / 2.0; + u_xy = u.x; + //Translation matrix 4x3 {{beta},{gamma},{ox},{oz}} for OX case + translation_matrix = {{0, 1, 0}, + {1, -1, 0}, + {1, 0, 0}, + {0, 0, 1}}; } -} + double beta_dir = u.x * cosa + u.y * sina; + double gamma_dir = u.x * cosa - u.y * sina; -//DR OX distance calculation -std::pair> -HexLattice::distance_ox(Position r, Direction u, const std::array& i_xyz) -const -{ - // Compute the direction on the hexagonal basis. - // DR for OX angle=60 g - double beta_dir = u.x*0.5 + u.y * std::sqrt(3.0) / 2.0; - double gamma_dir = u.x*0.5 - u.y * std::sqrt(3.0) / 2.0; + // Note that hexagonal lattice distance calculations are performed + // using the particle's coordinates relative to the neighbor lattice + // cells, not relative to the particle's current cell. This is done + // because there is significant disagreement between neighboring cells + // on where the lattice boundary is due to finite precision issues. - // Note that hexagonal lattice distance calculations are performed - // using the particle's coordinates relative to the neighbor lattice - // cells, not relative to the particle's current cell. This is done - // because there is significant disagreement between neighboring cells - // on where the lattice boundary is due to finite precision issues. - - // Upper-right and lower-left sides. - double d {INFTY}; - std::array lattice_trans; - double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge - Position r_t; - if (beta_dir > 0) { - const std::array i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double beta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0; - if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { - d = (edge - beta) / beta_dir; - if (beta_dir > 0) { - lattice_trans = {0, 1, 0}; - } else { - lattice_trans = {0, -1, 0}; - } - } - - // Lower-right and upper-left sides. - edge = -copysign(0.5*pitch_[0], gamma_dir); - if (gamma_dir > 0) { - const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0; - if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { - double this_d = (edge - gamma) / gamma_dir; - if (this_d < d) { - if (gamma_dir > 0) { - lattice_trans = {1, -1, 0}; - } else { - lattice_trans = {-1, 1, 0}; - } - d = this_d; - } - } - - // right and left sides. - edge = -copysign(0.5*pitch_[0], u.x); - if (u.x > 0) { - const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - if ((std::abs(r_t.x - edge) > FP_PRECISION) && u.x != 0) { - double this_d = (edge - r_t.x) / u.x; - if (this_d < d) { - if (u.x > 0) { - lattice_trans = {1, 0, 0}; - } else { - lattice_trans = {-1, 0, 0}; - } - d = this_d; - } - } - - // Top and bottom sides - if (is_3d_) { - double z = r.z; - double z0 {copysign(0.5 * pitch_[1], u.z)}; - if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { - double this_d = (z0 - z) / u.z; - if (this_d < d) { - d = this_d; - if (u.z > 0) { - lattice_trans = {0, 0, 1}; + //beta direction + double d {INFTY}; + std::array lattice_trans; + double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge + Position r_t; + if (beta_dir > 0) { + const std::array i_xyz_t {i_xyz[0]+translation_matrix(0,0), i_xyz[1]+translation_matrix(0,1), i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); } else { - lattice_trans = {0, 0, -1}; + const std::array i_xyz_t {i_xyz[0]-translation_matrix(0,0), i_xyz[1]-translation_matrix(0,1), i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); } - d = this_d; - } - } - } - - return {d, lattice_trans}; -} - -//============================================================================== - -//============================================================================== -//DR OY distance calculation -std::pair> -HexLattice::distance_oy(Position r, Direction u, const std::array& i_xyz) -const -{ - // Compute the direction on the hexagonal basis. - double beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; - double gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0; - - // Note that hexagonal lattice distance calculations are performed - // using the particle's coordinates relative to the neighbor lattice - // cells, not relative to the particle's current cell. This is done - // because there is significant disagreement between neighboring cells - // on where the lattice boundary is due to finite precision issues. - - // Upper-right and lower-left sides. - double d {INFTY}; - std::array lattice_trans; - double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge - Position r_t; - if (beta_dir > 0) { - const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0; - if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { - d = (edge - beta) / beta_dir; - if (beta_dir > 0) { - lattice_trans = {1, 0, 0}; - } else { - lattice_trans = {-1, 0, 0}; - } - } - - // Lower-right and upper-left sides. - edge = -copysign(0.5*pitch_[0], gamma_dir); - if (gamma_dir > 0) { - const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0; - if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { - double this_d = (edge - gamma) / gamma_dir; - if (this_d < d) { - if (gamma_dir > 0) { - lattice_trans = {1, -1, 0}; - } else { - lattice_trans = {-1, 1, 0}; - } - d = this_d; - } - } - - // Upper and lower sides. - edge = -copysign(0.5*pitch_[0], u.y); - if (u.y > 0) { - const std::array i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - if ((std::abs(r_t.y - edge) > FP_PRECISION) && u.y != 0) { - double this_d = (edge - r_t.y) / u.y; - if (this_d < d) { - if (u.y > 0) { - lattice_trans = {0, 1, 0}; - } else { - lattice_trans = {0, -1, 0}; - } - d = this_d; - } - } - - // Top and bottom sides - if (is_3d_) { - double z = r.z; - double z0 {copysign(0.5 * pitch_[1], u.z)}; - if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { - double this_d = (z0 - z) / u.z; - if (this_d < d) { - d = this_d; - if (u.z > 0) { - lattice_trans = {0, 0, 1}; + double beta = r_t.x * cosa + r_t.y * sina; + if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { + d = (edge - beta) / beta_dir; + if (beta_dir > 0) { + i_trans = 1; + } else { + i_trans = -1; + } + } + //gamma direction + edge = -copysign(0.5*pitch_[0], gamma_dir); + if (gamma_dir > 0) { + const std::array i_xyz_t {i_xyz[0]+translation_matrix(1,0), i_xyz[1]+translation_matrix(1,1), i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); } else { - lattice_trans = {0, 0, -1}; + const std::array i_xyz_t {i_xyz[0]-translation_matrix(1,0), i_xyz[1]-translation_matrix(1,1), i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); } - d = this_d; - } - } - } + double gamma = r_t.x * cosa - r_t.y * sina; + if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { + double this_d = (edge - gamma) / gamma_dir; + if (this_d < d) { + if (gamma_dir > 0) { + i_trans = 2; + } else { + i_trans = -2; + } + d = this_d; + } + } + //y or x direction + edge = -copysign(0.5*pitch_[0], u_xy); + if (u_xy > 0) { + const std::array i_xyz_t {i_xyz[0]+translation_matrix(2,0), i_xyz[1]+translation_matrix(2,1), i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0]-translation_matrix(2,0), i_xyz[1]-translation_matrix(2,1), i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double r_txy; + if (hextype==0){ + r_txy = r_t.y; + } + else { + r_txy = r_t.x; + } + if ((std::abs(r_txy - edge) > FP_PRECISION) && u_xy != 0) { + double this_d = (edge - r_txy) / u_xy; + if (this_d < d) { + if (u_xy > 0) { + i_trans = 3; + } else { + i_trans = -3; + } + d = this_d; + } + } + // Top and bottom sides + if (is_3d_) { + double z = r.z; + double z0 {copysign(0.5 * pitch_[1], u.z)}; + if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { + double this_d = (z0 - z) / u.z; + if (this_d < d) { + d = this_d; + if (u.z > 0) { + i_trans = 4; + } else { + i_trans = -4; + } + d = this_d; + } + } + } - return {d, lattice_trans}; + for (int i=0;i<3;i++){ + if (i_trans < 0){ + lattice_trans[i] = -translation_matrix(-(i_trans+1),i); + }else{ + lattice_trans[i] = translation_matrix((i_trans-1),i); + } + } + + return {d, lattice_trans}; } //============================================================================== From 4ef3e307559d8404396b837ccf3ec61caa1f6931 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Mon, 27 May 2019 10:54:44 +0300 Subject: [PATCH 07/31] lattice translation with S.Harper basis --- src/lattice.cpp | 251 ++++++++++++++++++++++++------------------------ 1 file changed, 127 insertions(+), 124 deletions(-) diff --git a/src/lattice.cpp b/src/lattice.cpp index af025033a..5a3fea58a 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -694,134 +694,137 @@ std::pair> HexLattice::distance(Position r, Direction u, const std::array& i_xyz) const { - double cosa,sina,u_xy; - // Index + 1 of translation direction - int i_trans {1}; - xt :: xarray translation_matrix; - if (hextype==0){ - cosa = std::sqrt(3.0) / 2.0; - sina = 0.5; - u_xy = u.y; - //Translation matrix 4x3 {{beta},{gamma},{oy},{oz}} for OY case - translation_matrix = {{1, 0, 0}, - {1, -1, 0}, - {0, 1, 0}, - {0, 0, 1}}; - } else { - cosa = 0.5; - sina = std::sqrt(3.0) / 2.0; - u_xy = u.x; - //Translation matrix 4x3 {{beta},{gamma},{ox},{oz}} for OX case - translation_matrix = {{0, 1, 0}, - {1, -1, 0}, - {1, 0, 0}, - {0, 0, 1}}; - } - double beta_dir = u.x * cosa + u.y * sina; - double gamma_dir = u.x * cosa - u.y * sina; + //Short description + //OY - orientation: + // basis0 = (1, 0) + // basis1 = (-1/sqrt(3), 1) = +120 degrees from basis0 + // beta = (sqrt(3)/2, 1/2) = +30 degrees from basis0 + // gamma = (sqrt(3)/2, -1/2) = -60 degrees from beta + // delta = (0, 1) = +60 degrees from beta + //OX - orientation: + // basis0 = (1/sqrt(3), -1) + // basis1 = (0, 1) = +120 degrees from basis0 + // beta = (1, 0) = +30 degrees from basis0 + // gamma = (1/2, -sqrt(3)/2) = -60 degrees from beta + // delta = (1/2, sqrt(3)/2) = +60 degrees from beta + //OZ be considered separetly + double beta_dir; + double gamma_dir; + double delta_dir; + if (hextype==0){ + beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; + gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0; + delta_dir = u.y; + } + else{ + beta_dir = u.x; + gamma_dir = u.x / 2.0 - u.y * std::sqrt(3.0) / 2.0; + delta_dir = u.x / 2.0 + u.y * std::sqrt(3.0) / 2.0; + } - // Note that hexagonal lattice distance calculations are performed - // using the particle's coordinates relative to the neighbor lattice - // cells, not relative to the particle's current cell. This is done - // because there is significant disagreement between neighboring cells - // on where the lattice boundary is due to finite precision issues. + // Note that hexagonal lattice distance calculations are performed + // using the particle's coordinates relative to the neighbor lattice + // cells, not relative to the particle's current cell. This is done + // because there is significant disagreement between neighboring cells + // on where the lattice boundary is due to finite precision issues. - //beta direction - double d {INFTY}; - std::array lattice_trans; - double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge - Position r_t; - if (beta_dir > 0) { - const std::array i_xyz_t {i_xyz[0]+translation_matrix(0,0), i_xyz[1]+translation_matrix(0,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-translation_matrix(0,0), i_xyz[1]-translation_matrix(0,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double beta = r_t.x * cosa + r_t.y * sina; - if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { - d = (edge - beta) / beta_dir; - if (beta_dir > 0) { - i_trans = 1; - } else { - i_trans = -1; - } - } - //gamma direction - edge = -copysign(0.5*pitch_[0], gamma_dir); - if (gamma_dir > 0) { - const std::array i_xyz_t {i_xyz[0]+translation_matrix(1,0), i_xyz[1]+translation_matrix(1,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-translation_matrix(1,0), i_xyz[1]-translation_matrix(1,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double gamma = r_t.x * cosa - r_t.y * sina; - if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { - double this_d = (edge - gamma) / gamma_dir; - if (this_d < d) { - if (gamma_dir > 0) { - i_trans = 2; - } else { - i_trans = -2; - } - d = this_d; - } - } - //y or x direction - edge = -copysign(0.5*pitch_[0], u_xy); - if (u_xy > 0) { - const std::array i_xyz_t {i_xyz[0]+translation_matrix(2,0), i_xyz[1]+translation_matrix(2,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-translation_matrix(2,0), i_xyz[1]-translation_matrix(2,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double r_txy; - if (hextype==0){ - r_txy = r_t.y; - } - else { - r_txy = r_t.x; - } - if ((std::abs(r_txy - edge) > FP_PRECISION) && u_xy != 0) { - double this_d = (edge - r_txy) / u_xy; - if (this_d < d) { - if (u_xy > 0) { - i_trans = 3; - } else { - i_trans = -3; - } - d = this_d; - } - } - // Top and bottom sides - if (is_3d_) { - double z = r.z; - double z0 {copysign(0.5 * pitch_[1], u.z)}; - if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { - double this_d = (z0 - z) / u.z; - if (this_d < d) { - d = this_d; - if (u.z > 0) { - i_trans = 4; - } else { - i_trans = -4; - } - d = this_d; - } - } - } + //beta direction + double d {INFTY}; + std::array lattice_trans; + double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge + double beta; + Position r_t; + if (beta_dir > 0) { + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + if (hextype==0){ + beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0;} + else{ + beta = r_t.x ;} + if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { + d = (edge - beta) / beta_dir; + if (beta_dir > 0) { + lattice_trans = {1, 0, 0}; + } else { + lattice_trans = {-1, 0, 0}; + } + } - for (int i=0;i<3;i++){ - if (i_trans < 0){ - lattice_trans[i] = -translation_matrix(-(i_trans+1),i); - }else{ - lattice_trans[i] = translation_matrix((i_trans-1),i); - } - } + // gamma direction. + edge = -copysign(0.5*pitch_[0], gamma_dir); + if (gamma_dir > 0) { + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double gamma; + if (hextype==0){ + gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0;} + else{ + gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0;} + if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { + double this_d = (edge - gamma) / gamma_dir; + if (this_d < d) { + if (gamma_dir > 0) { + lattice_trans = {1, -1, 0}; + } else { + lattice_trans = {-1, 1, 0}; + } + d = this_d; + } + } - return {d, lattice_trans}; + // delta directions. + edge = -copysign(0.5*pitch_[0], delta_dir); + if (delta_dir > 0) { + const std::array i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double delta; + if (hextype==0){ + delta = r_t.y ;} + else{ + delta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0;} + if ((std::abs(delta - edge) > FP_PRECISION) && delta_dir != 0) { + double this_d = (edge - delta) / delta_dir; + if (this_d < d) { + if (delta_dir > 0) { + lattice_trans = {0, 1, 0}; + } else { + lattice_trans = {0, -1, 0}; + } + d = this_d; + } + } + + // Top and bottom sides + if (is_3d_) { + double z = r.z; + double z0 {copysign(0.5 * pitch_[1], u.z)}; + if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { + double this_d = (z0 - z) / u.z; + if (this_d < d) { + d = this_d; + if (u.z > 0) { + lattice_trans = {0, 0, 1}; + } else { + lattice_trans = {0, 0, -1}; + } + d = this_d; + } + } + } + + return {d, lattice_trans}; } //============================================================================== From f8690b87a0f901b212689b2c817f2d3f9ad61837 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 27 May 2019 14:59:05 -0400 Subject: [PATCH 08/31] Modify code to comply with the style guide --- include/openmc/lattice.h | 18 +- src/lattice.cpp | 627 +++++++++++++++++++-------------------- 2 files changed, 322 insertions(+), 323 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index afec39b22..449b92164 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -248,30 +248,32 @@ public: bool are_valid_indices(const int i_xyz[3]) const; - std::pair> distance(Position r, Direction u, const std::array& i_xyz) const; + std::pair> + distance(Position r, Direction u, const std::array& i_xyz) const; std::array get_indices(Position r, Direction u) const; - Position get_local_position(Position r, const std::array i_xyz) const; bool is_valid_index(int indx) const; - //DR fill lattice with OY by default orientation - void fill_lattice_oy(std::vector univ_words); - //DR fill lattice with OX orientation - void fill_lattice_ox(std::vector univ_words); int32_t& offset(int map, const int i_xyz[3]); std::string index_to_string(int indx) const; - //TODO: find a desicion for compatible old hdf5 file without orientation + void to_hdf5_inner(hid_t group_id) const; private: + //! Fill universes_ vector for OY orientation + void fill_lattice_oy(std::vector univ_words); + + //! Fill universes_ vector for OX orientation + void fill_lattice_ox(std::vector univ_words); + int n_rings_; //!< Number of radial tile positions int n_axial_; //!< Number of axial tile positions - int hextype; //!DR type of hexagonal lattice orientation //0 - OY by default;1-OX + int hextype; //!< Type of hexagonal lattice orientation //0 - OY by default;1-OX Position center_; //!< Global center of lattice std::array pitch_; //!< Lattice tile width and height }; diff --git a/src/lattice.cpp b/src/lattice.cpp index 5a3fea58a..314670939 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -430,8 +430,6 @@ HexLattice::HexLattice(pugi::xml_node lat_node) { type_ = LatticeType::hex; - hextype = 0;//DR a default OY orientation of lattice - // Read the number of lattice cells in each dimension. n_rings_ = std::stoi(get_node_value(lat_node, "n_rings")); if (check_for_node(lat_node, "n_axial")) { @@ -442,10 +440,19 @@ HexLattice::HexLattice(pugi::xml_node lat_node) is_3d_ = false; } - //DR Check if "OX" orientation desc in node + // Read the lattice orientation. Default to OY. if (check_for_node(lat_node, "orient")) { - std :: string orientation = get_node_value(lat_node, "orient"); - if (orientation == "OX") {hextype=1;} + std::string orientation = get_node_value(lat_node, "orient"); + if (orientation == "OY") { + hextype = 0; + } else if (orientation == "OX") { + hextype = 1; + } else { + fatal_error("Unrecognized orientation '" + orientation + + "' for lattice " + std::to_string(id_)); + } + } else { + hextype = 0; } // Read the lattice center. @@ -496,167 +503,159 @@ HexLattice::HexLattice(pugi::xml_node lat_node) // input and the order that they will be stored in the skewed array so // the following code walks a set of index values across the skewed array // in a manner that matches the input order. Note that i_x = 0, i_a = 0 - // or i_a=0 , i_y = 0 corresponds to the center of the hexagonal lattice. + // or i_a = 0, i_y = 0 corresponds to the center of the hexagonal lattice. universes_.resize((2*n_rings_-1) * (2*n_rings_-1) * n_axial_, C_NONE); if (hextype==0) { fill_lattice_oy(univ_words); - } else - { + } else { fill_lattice_ox(univ_words); } } -// DR fill universes in OX orientation + //============================================================================== + void -HexLattice::fill_lattice_ox(std::vector univ_words){ +HexLattice::fill_lattice_ox(std::vector univ_words) +{ + int input_index = 0; + for (int m = 0; m < n_axial_; m++) { + // Initialize lattice indecies. + int i_a = -(n_rings_ - 1); + int i_y = n_rings_ - 1; - int input_index = 0; - for (int m = 0; m < n_axial_; m++) { - // Initialize lattice indecies. - int i_a = -(n_rings_ - 1); - int i_y = n_rings_ - 1; + // Map upper region of hexagonal lattice which is found in the + // first n_rings-1 rows of the input. + for (int k = 0; k < n_rings_-1; k++) { - - //DR Map upper region of hexagonal lattice which is found in the - // first n_rings-1 rows of the input. - for (int k = 0; k < n_rings_-1; k++) { - - // Iterate over the input columns. - for (int j = 0; j < k+n_rings_; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_y+n_rings_-1) - + (i_a+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - //DR Move to the next right neighbour cell - i_a += 1; - - } - - // Return the lattice index to the start of the current row. - i_a = -(n_rings_ - 1); - i_y -= 1; - } - - // Map the lower region from the centerline of cart to down side - for (int k = 0; k < n_rings_; k++) { - // Walk the index to the lower-right neighbor of the last row start. - i_a = -(n_rings_ - 1) + k; - - // Iterate over the input columns. - for (int j = 0; j < 2*n_rings_-k-1; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_y+n_rings_-1) - + (i_a+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - //DR Move to the next right neighbour cell - i_a += 1; - - } - - // Return lattice index to start of current row. - i_y -= 1; - } + // Iterate over the input columns. + for (int j = 0; j < k+n_rings_; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_y+n_rings_-1) + + (i_a+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Move to the next right neighbour cell + i_a += 1; } + + // Return the lattice index to the start of the current row. + i_a = -(n_rings_ - 1); + i_y -= 1; + } + + // Map the lower region from the centerline of cart to down side + for (int k = 0; k < n_rings_; k++) { + // Walk the index to the lower-right neighbor of the last row start. + i_a = -(n_rings_ - 1) + k; + + // Iterate over the input columns. + for (int j = 0; j < 2*n_rings_-k-1; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_y+n_rings_-1) + + (i_a+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Move to the next right neighbour cell + i_a += 1; + } + + // Return lattice index to start of current row. + ki_y -= 1; + } + } } //============================================================================== - -// DR fill universes in OY orientation -//============================================================================== void -HexLattice::fill_lattice_oy(std::vector univ_words){ +HexLattice::fill_lattice_oy(std::vector univ_words) +{ + int input_index = 0; + for (int m = 0; m < n_axial_; m++) { + // Initialize lattice indecies. + int i_x = 1; + int i_a = n_rings_ - 1; - int input_index = 0; - for (int m = 0; m < n_axial_; m++) { - // Initialize lattice indecies. - int i_x = 1; - int i_a = n_rings_ - 1; + // Map upper triangular region of hexagonal lattice which is found in the + // first n_rings-1 rows of the input. + for (int k = 0; k < n_rings_-1; k++) { + // Walk the index to lower-left neighbor of last row start. + i_x -= 1; - // Map upper triangular region of hexagonal lattice which is found in the - // first n_rings-1 rows of the input. - for (int k = 0; k < n_rings_-1; k++) { - // Walk the index to lower-left neighbor of last row start. - i_x -= 1; - - // Iterate over the input columns. - for (int j = 0; j < k+1; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_a+n_rings_-1) - + (i_x+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - // Walk the index to the right neighbor (which is not adjacent). - i_x += 2; - i_a -= 1; - } - - // Return the lattice index to the start of the current row. - i_x -= 2 * (k+1); - i_a += (k+1); - } - - // Map the middle square region of the hexagonal lattice which is found in - // the next 2*n_rings-1 rows of the input. - for (int k = 0; k < 2*n_rings_-1; k++) { - if ((k % 2) == 0) { - // Walk the index to the lower-left neighbor of the last row start. - i_x -= 1; - } else { - // Walk the index to the lower-right neighbor of the last row start. - i_x += 1; - i_a -= 1; - } - - // Iterate over the input columns. - for (int j = 0; j < n_rings_ - (k % 2); j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_a+n_rings_-1) - + (i_x+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - // Walk the index to the right neighbor (which is not adjacent). - i_x += 2; - i_a -= 1; - } - - // Return the lattice index to the start of the current row. - i_x -= 2*(n_rings_ - (k % 2)); - i_a += n_rings_ - (k % 2); - } - - // Map the lower triangular region of the hexagonal lattice. - for (int k = 0; k < n_rings_-1; k++) { - // Walk the index to the lower-right neighbor of the last row start. - i_x += 1; - i_a -= 1; - - // Iterate over the input columns. - for (int j = 0; j < n_rings_-k-1; j++) { - int indx = (2*n_rings_-1)*(2*n_rings_-1) * m - + (2*n_rings_-1) * (i_a+n_rings_-1) - + (i_x+n_rings_-1); - universes_[indx] = std::stoi(univ_words[input_index]); - input_index++; - // Walk the index to the right neighbor (which is not adjacent). - i_x += 2; - i_a -= 1; - } - - // Return lattice index to start of current row. - i_x -= 2*(n_rings_ - k - 1); - i_a += n_rings_ - k - 1; - } + // Iterate over the input columns. + for (int j = 0; j < k+1; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_a+n_rings_-1) + + (i_x+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; } -} -//============================================================================== + // Return the lattice index to the start of the current row. + i_x -= 2 * (k+1); + i_a += (k+1); + } + + // Map the middle square region of the hexagonal lattice which is found in + // the next 2*n_rings-1 rows of the input. + for (int k = 0; k < 2*n_rings_-1; k++) { + if ((k % 2) == 0) { + // Walk the index to the lower-left neighbor of the last row start. + i_x -= 1; + } else { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; + } + + // Iterate over the input columns. + for (int j = 0; j < n_rings_ - (k % 2); j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_a+n_rings_-1) + + (i_x+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return the lattice index to the start of the current row. + i_x -= 2*(n_rings_ - (k % 2)); + i_a += n_rings_ - (k % 2); + } + + // Map the lower triangular region of the hexagonal lattice. + for (int k = 0; k < n_rings_-1; k++) { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; + + // Iterate over the input columns. + for (int j = 0; j < n_rings_-k-1; j++) { + int indx = (2*n_rings_-1)*(2*n_rings_-1) * m + + (2*n_rings_-1) * (i_a+n_rings_-1) + + (i_x+n_rings_-1); + universes_[indx] = std::stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return lattice index to start of current row. + i_x -= 2*(n_rings_ - k - 1); + i_a += n_rings_ - k - 1; + } + } +} //============================================================================== @@ -689,145 +688,148 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const && (i_xyz[2] < n_axial_)); } -//DR wrapper HexLattice function both for OX and OY orientation +//============================================================================== + std::pair> HexLattice::distance(Position r, Direction u, const std::array& i_xyz) const { - //Short description - //OY - orientation: - // basis0 = (1, 0) - // basis1 = (-1/sqrt(3), 1) = +120 degrees from basis0 - // beta = (sqrt(3)/2, 1/2) = +30 degrees from basis0 - // gamma = (sqrt(3)/2, -1/2) = -60 degrees from beta - // delta = (0, 1) = +60 degrees from beta - //OX - orientation: - // basis0 = (1/sqrt(3), -1) - // basis1 = (0, 1) = +120 degrees from basis0 - // beta = (1, 0) = +30 degrees from basis0 - // gamma = (1/2, -sqrt(3)/2) = -60 degrees from beta - // delta = (1/2, sqrt(3)/2) = +60 degrees from beta - //OZ be considered separetly - double beta_dir; - double gamma_dir; - double delta_dir; - if (hextype==0){ - beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; - gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0; - delta_dir = u.y; - } - else{ - beta_dir = u.x; - gamma_dir = u.x / 2.0 - u.y * std::sqrt(3.0) / 2.0; - delta_dir = u.x / 2.0 + u.y * std::sqrt(3.0) / 2.0; - } + // Short description of the direction vectors used here. The beta, gamma, and + // delta vectors point towards the flat sides of each hexagonal tile. + // OY - orientation: + // basis0 = (1, 0) + // basis1 = (-1/sqrt(3), 1) = +120 degrees from basis0 + // beta = (sqrt(3)/2, 1/2) = +30 degrees from basis0 + // gamma = (sqrt(3)/2, -1/2) = -60 degrees from beta + // delta = (0, 1) = +60 degrees from beta + // OX - orientation: + // basis0 = (1/sqrt(3), -1) + // basis1 = (0, 1) = +120 degrees from basis0 + // beta = (1, 0) = +30 degrees from basis0 + // gamma = (1/2, -sqrt(3)/2) = -60 degrees from beta + // delta = (1/2, sqrt(3)/2) = +60 degrees from beta + // OZ is considered separately + double beta_dir; + double gamma_dir; + double delta_dir; + if (hextype == 0){ + beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; + gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0; + delta_dir = u.y; + } else { + beta_dir = u.x; + gamma_dir = u.x / 2.0 - u.y * std::sqrt(3.0) / 2.0; + delta_dir = u.x / 2.0 + u.y * std::sqrt(3.0) / 2.0; + } - // Note that hexagonal lattice distance calculations are performed - // using the particle's coordinates relative to the neighbor lattice - // cells, not relative to the particle's current cell. This is done - // because there is significant disagreement between neighboring cells - // on where the lattice boundary is due to finite precision issues. + // Note that hexagonal lattice distance calculations are performed + // using the particle's coordinates relative to the neighbor lattice + // cells, not relative to the particle's current cell. This is done + // because there is significant disagreement between neighboring cells + // on where the lattice boundary is due to finite precision issues. - //beta direction - double d {INFTY}; - std::array lattice_trans; - double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge - double beta; - Position r_t; - if (beta_dir > 0) { - const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - if (hextype==0){ - beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0;} - else{ - beta = r_t.x ;} - if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { - d = (edge - beta) / beta_dir; - if (beta_dir > 0) { - lattice_trans = {1, 0, 0}; - } else { - lattice_trans = {-1, 0, 0}; - } - } + // beta direction + double d {INFTY}; + std::array lattice_trans; + double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge + Position r_t; + if (beta_dir > 0) { + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double beta; + if (hextype == 0) { + beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0; + } else { + beta = r_t.x; + } + if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { + d = (edge - beta) / beta_dir; + if (beta_dir > 0) { + lattice_trans = {1, 0, 0}; + } else { + lattice_trans = {-1, 0, 0}; + } + } - // gamma direction. - edge = -copysign(0.5*pitch_[0], gamma_dir); - if (gamma_dir > 0) { - const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double gamma; - if (hextype==0){ - gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0;} - else{ - gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0;} - if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { - double this_d = (edge - gamma) / gamma_dir; - if (this_d < d) { - if (gamma_dir > 0) { - lattice_trans = {1, -1, 0}; - } else { - lattice_trans = {-1, 1, 0}; - } - d = this_d; - } - } + // gamma direction + edge = -copysign(0.5*pitch_[0], gamma_dir); + if (gamma_dir > 0) { + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double gamma; + if (hextype == 0) { + gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0; + } else { + gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0; + } + if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { + double this_d = (edge - gamma) / gamma_dir; + if (this_d < d) { + if (gamma_dir > 0) { + lattice_trans = {1, -1, 0}; + } else { + lattice_trans = {-1, 1, 0}; + } + d = this_d; + } + } - // delta directions. - edge = -copysign(0.5*pitch_[0], delta_dir); - if (delta_dir > 0) { - const std::array i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double delta; - if (hextype==0){ - delta = r_t.y ;} - else{ - delta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0;} - if ((std::abs(delta - edge) > FP_PRECISION) && delta_dir != 0) { - double this_d = (edge - delta) / delta_dir; - if (this_d < d) { - if (delta_dir > 0) { - lattice_trans = {0, 1, 0}; - } else { - lattice_trans = {0, -1, 0}; - } - d = this_d; - } - } + // delta direction + edge = -copysign(0.5*pitch_[0], delta_dir); + if (delta_dir > 0) { + const std::array i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double delta; + if (hextype == 0){ + delta = r_t.y; + } else { + delta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0; + } + if ((std::abs(delta - edge) > FP_PRECISION) && delta_dir != 0) { + double this_d = (edge - delta) / delta_dir; + if (this_d < d) { + if (delta_dir > 0) { + lattice_trans = {0, 1, 0}; + } else { + lattice_trans = {0, -1, 0}; + } + d = this_d; + } + } - // Top and bottom sides - if (is_3d_) { - double z = r.z; - double z0 {copysign(0.5 * pitch_[1], u.z)}; - if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { - double this_d = (z0 - z) / u.z; - if (this_d < d) { - d = this_d; - if (u.z > 0) { - lattice_trans = {0, 0, 1}; - } else { - lattice_trans = {0, 0, -1}; - } - d = this_d; - } - } - } + // Top and bottom sides + if (is_3d_) { + double z = r.z; + double z0 {copysign(0.5 * pitch_[1], u.z)}; + if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { + double this_d = (z0 - z) / u.z; + if (this_d < d) { + d = this_d; + if (u.z > 0) { + lattice_trans = {0, 0, 1}; + } else { + lattice_trans = {0, 0, -1}; + } + d = this_d; + } + } + } - return {d, lattice_trans}; + return {d, lattice_trans}; } -//============================================================================== //============================================================================== std::array @@ -855,25 +857,23 @@ HexLattice::get_indices(Position r, Direction u) const } } - int i1{};// in case of OY orientation such index means x - coordinate in OY - alpha - int i2{};// in case of OX orientation such index means alpha - coordinate in OY - y - if (hextype==0){//DR type implement OY default behaviour - // Convert coordinates into skewed bases. The (x, alpha) basis is used to - // find the index of the global coordinates to within 4 cells. - double alpha = r_o.y - r_o.x / std::sqrt(3.0); - i1 = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0])); - i2 = std::floor(alpha / pitch_[0]);} - else{//DR type implement OX (alpha,y) basis - // Convert coordinates into skewed bases. The (alpha,y) basis is used to - // find the index of the global coordinates to within 4 cells. - double alpha = r_o.y - r_o.x * std::sqrt(3.0); - - i1 = std::floor(-alpha / (std::sqrt(3.0) * pitch_[0])); - i2 = std::floor(r_o.y / (0.5*std::sqrt(3.0) * pitch_[0])); + int i1, i2; + if (hextype == 0) { + // Convert coordinates into skewed bases. The (x, alpha) basis is used to + // find the index of the global coordinates to within 4 cells. + double alpha = r_o.y - r_o.x / std::sqrt(3.0); + i1 = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0])); + i2 = std::floor(alpha / pitch_[0]); + } else { + // Convert coordinates into skewed bases. The (alpha, y) basis is used to + // find the index of the global coordinates to within 4 cells. + double alpha = r_o.y - r_o.x * std::sqrt(3.0); + i1 = std::floor(-alpha / (std::sqrt(3.0) * pitch_[0])); + i2 = std::floor(r_o.y / (0.5*std::sqrt(3.0) * pitch_[0])); } - // Add offset to indices (the center cell is (i_x/i_alpha, i_alpha/i_y) = (0, 0) for OY/OX - // orientation but the array is offset so that the indices never go below 0). + // Add offset to indices (the center cell is (i1, i2) = (0, 0) but + // the array is offset so that the indices never go below 0). i1 += n_rings_-1; i2 += n_rings_-1; @@ -936,32 +936,30 @@ Position HexLattice::get_local_position(Position r, const std::array i_xyz) const { - // DR for OY implementation - if (hextype == 0){ - // x_l = x_g - (center + pitch_x*cos(30)*index_x) - r.x -= center_.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0]; - // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) - r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0] - + (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0); - if (is_3d_) { - r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; - } + if (hextype == 0) { + // x_l = x_g - (center + pitch_x*cos(30)*index_x) + r.x -= center_.x + + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0]; + // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) + r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0] + + (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0); + if (is_3d_) { + r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; + } - } - else {// DR for OX implementation + } else { + // x_l = x_g - (center + pitch_x*index_a + pitch_y*sin(30)*index_y) + r.x -= (center_.x + (i_xyz[0] - n_rings_ + 1) * pitch_[0] + + (i_xyz[1] - n_rings_ + 1) * pitch_[0] / 2.0); + // y_l = y_g - (center + pitch_y*cos(30)*index_y) + r.y -= center_.y + + std::sqrt(3.0)/2.0 * (i_xyz[1] - n_rings_ + 1) * pitch_[0]; + if (is_3d_) { + r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; + } + } - // x_l = x_g - (center + pitch_x*index_a + pitch_y*sin(30)*index_y) - r.x -= (center_.x + (i_xyz[0] - n_rings_ + 1) * pitch_[0] - + (i_xyz[1] - n_rings_ + 1) * pitch_[0] / 2.0); - // y_l = y_g - (center + pitch_y*cos(30)*index_y) - r.y -= center_.y + std::sqrt(3.0)/2.0 * (i_xyz[1] - n_rings_ + 1) * pitch_[0]; - - if (is_3d_) { - r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; - } - - } - return r; + return r; } //============================================================================== @@ -1018,11 +1016,10 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const write_string(lat_group, "type", "hexagonal", false); write_dataset(lat_group, "n_rings", n_rings_); write_dataset(lat_group, "n_axial", n_axial_); - if (hextype == 0){ - write_string(lat_group, "orientation", "oy",false); - } - else { - write_string(lat_group, "orientation", "ox",false); + if (hextype == 0) { + write_string(lat_group, "orientation", "oy", false); + } else { + write_string(lat_group, "orientation", "ox", false); } if (is_3d_) { write_dataset(lat_group, "pitch", pitch_); From 2394920f7dc11a12294167f3638c3ed0749a88e3 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 27 May 2019 15:03:00 -0400 Subject: [PATCH 09/31] Remove TINY_BIT to reestablish #1211 changes --- src/lattice.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/lattice.cpp b/src/lattice.cpp index 314670939..40f30af2c 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -835,12 +835,6 @@ const std::array HexLattice::get_indices(Position r, Direction u) const { - - // The implementation for HexLattice currently doesn't use direction - // information. As a result, we move the position slightly forward to - // determine what lattice index the particle is most likely to be in. - r += TINY_BIT * u; - // Offset the xyz by the lattice center. Position r_o {r.x - center_.x, r.y - center_.y, r.z}; if (is_3d_) {r_o.z -= center_.z;} From e442f521029df4b63c494b2b453a669832426bab Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 27 May 2019 15:16:37 -0400 Subject: [PATCH 10/31] Use enum for hex lattice orientation --- include/openmc/lattice.h | 6 +++++- src/lattice.cpp | 26 +++++++++++++------------- 2 files changed, 18 insertions(+), 14 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index 449b92164..6b30486b8 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -27,6 +27,10 @@ enum class LatticeType { rect, hex }; +enum class HexOrientation { + oy, ox +}; + //============================================================================== // Global variables //============================================================================== @@ -273,7 +277,7 @@ private: int n_rings_; //!< Number of radial tile positions int n_axial_; //!< Number of axial tile positions - int hextype; //!< Type of hexagonal lattice orientation //0 - OY by default;1-OX + HexOrientation orientation_; //!< Flat side up (oy) vs. sharp side up (ox) Position center_; //!< Global center of lattice std::array pitch_; //!< Lattice tile width and height }; diff --git a/src/lattice.cpp b/src/lattice.cpp index 40f30af2c..b98ba6366 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -444,15 +444,15 @@ HexLattice::HexLattice(pugi::xml_node lat_node) if (check_for_node(lat_node, "orient")) { std::string orientation = get_node_value(lat_node, "orient"); if (orientation == "OY") { - hextype = 0; + orientation_ = HexOrientation::oy; } else if (orientation == "OX") { - hextype = 1; + orientation_ = HexOrientation::ox; } else { fatal_error("Unrecognized orientation '" + orientation + "' for lattice " + std::to_string(id_)); } } else { - hextype = 0; + orientation_ = HexOrientation::oy; } // Read the lattice center. @@ -507,7 +507,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node) universes_.resize((2*n_rings_-1) * (2*n_rings_-1) * n_axial_, C_NONE); - if (hextype==0) { + if (orientation_ == HexOrientation::oy) { fill_lattice_oy(univ_words); } else { fill_lattice_ox(univ_words); @@ -564,7 +564,7 @@ HexLattice::fill_lattice_ox(std::vector univ_words) } // Return lattice index to start of current row. - ki_y -= 1; + i_y -= 1; } } } @@ -708,11 +708,11 @@ const // beta = (1, 0) = +30 degrees from basis0 // gamma = (1/2, -sqrt(3)/2) = -60 degrees from beta // delta = (1/2, sqrt(3)/2) = +60 degrees from beta - // OZ is considered separately + // The z-axis is considered separately. double beta_dir; double gamma_dir; double delta_dir; - if (hextype == 0){ + if (orientation_ == HexOrientation::oy) { beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0; delta_dir = u.y; @@ -741,7 +741,7 @@ const r_t = get_local_position(r, i_xyz_t); } double beta; - if (hextype == 0) { + if (orientation_ == HexOrientation::oy) { beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0; } else { beta = r_t.x; @@ -765,7 +765,7 @@ const r_t = get_local_position(r, i_xyz_t); } double gamma; - if (hextype == 0) { + if (orientation_ == HexOrientation::oy) { gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0; } else { gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0; @@ -792,7 +792,7 @@ const r_t = get_local_position(r, i_xyz_t); } double delta; - if (hextype == 0){ + if (orientation_ == HexOrientation::oy) { delta = r_t.y; } else { delta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0; @@ -852,7 +852,7 @@ HexLattice::get_indices(Position r, Direction u) const } int i1, i2; - if (hextype == 0) { + if (orientation_ == HexOrientation::oy) { // Convert coordinates into skewed bases. The (x, alpha) basis is used to // find the index of the global coordinates to within 4 cells. double alpha = r_o.y - r_o.x / std::sqrt(3.0); @@ -930,7 +930,7 @@ Position HexLattice::get_local_position(Position r, const std::array i_xyz) const { - if (hextype == 0) { + if (orientation_ == HexOrientation::oy) { // x_l = x_g - (center + pitch_x*cos(30)*index_x) r.x -= center_.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0]; @@ -1010,7 +1010,7 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const write_string(lat_group, "type", "hexagonal", false); write_dataset(lat_group, "n_rings", n_rings_); write_dataset(lat_group, "n_axial", n_axial_); - if (hextype == 0) { + if (orientation_ == HexOrientation::oy) { write_string(lat_group, "orientation", "oy", false); } else { write_string(lat_group, "orientation", "ox", false); From 79f792c09a7e0ef72584f350fe79ad96bdeae847 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Thu, 30 May 2019 10:54:08 +0300 Subject: [PATCH 11/31] New regression test for lattice with OX orient --- .../lattice_hex_OX/__init__.py | 0 .../lattice_hex_OX/inputs_true.dat | 126 +++++++++++++ .../lattice_hex_OX/results_true.dat | 2 + tests/regression_tests/lattice_hex_OX/test.py | 171 ++++++++++++++++++ 4 files changed, 299 insertions(+) create mode 100644 tests/regression_tests/lattice_hex_OX/__init__.py create mode 100644 tests/regression_tests/lattice_hex_OX/inputs_true.dat create mode 100644 tests/regression_tests/lattice_hex_OX/results_true.dat create mode 100644 tests/regression_tests/lattice_hex_OX/test.py diff --git a/tests/regression_tests/lattice_hex_OX/__init__.py b/tests/regression_tests/lattice_hex_OX/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/lattice_hex_OX/inputs_true.dat b/tests/regression_tests/lattice_hex_OX/inputs_true.dat new file mode 100644 index 000000000..9130cf284 --- /dev/null +++ b/tests/regression_tests/lattice_hex_OX/inputs_true.dat @@ -0,0 +1,126 @@ + + + + + + + + + + + + + + + + + 1.235 5.0 + 4 +
0.0 0.0 5.0
+ + 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 + 1 1 1 1 1 3 1 1 1 1 3 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 3 1 1 1 1 3 1 1 1 3 1 1 1 1 + 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 +1 1 1 1 1 3 1 1 1 1 2 1 1 1 1 3 1 1 1 1 1 + 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 + 1 1 1 1 3 1 1 1 3 1 1 1 1 3 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 3 1 1 1 1 3 1 1 1 1 1 + 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 + 1 1 1 1 1 3 1 1 1 1 3 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 3 1 1 1 1 3 1 1 1 3 1 1 1 1 + 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 +1 1 1 1 1 3 1 1 1 1 2 1 1 1 1 3 1 1 1 1 1 + 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 + 1 1 1 1 3 1 1 1 3 1 1 1 1 3 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 3 1 1 1 1 3 1 1 1 1 1 + 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 5 + + + -13.62546635287517 -13.62546635287517 0.0 13.62546635287517 13.62546635287517 10.0 + + + + false + + 22 + diff --git a/tests/regression_tests/lattice_hex_OX/results_true.dat b/tests/regression_tests/lattice_hex_OX/results_true.dat new file mode 100644 index 000000000..a742aa0d5 --- /dev/null +++ b/tests/regression_tests/lattice_hex_OX/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.355663E+00 2.896562E-02 diff --git a/tests/regression_tests/lattice_hex_OX/test.py b/tests/regression_tests/lattice_hex_OX/test.py new file mode 100644 index 000000000..2eab458dc --- /dev/null +++ b/tests/regression_tests/lattice_hex_OX/test.py @@ -0,0 +1,171 @@ +from tests.testing_harness import PyAPITestHarness +import openmc +import numpy as np + +class HexLatticeOXTestHarness(PyAPITestHarness): + def _build_inputs(self): + materials = openmc.Materials() + + fuel_mat = openmc.Material(material_id=1,name="UO2") + fuel_mat.set_density('sum') + fuel_mat.add_nuclide('U235',0.87370e-03) + fuel_mat.add_nuclide('U238',1.87440e-02) + fuel_mat.add_nuclide('O16',3.92350e-02) + materials.append(fuel_mat) + + coolant = openmc.Material(material_id=2,name="borated H2O") + coolant.set_density('sum') + coolant.add_nuclide('H1',0.06694) + coolant.add_nuclide('O16',0.03347) + coolant.add_nuclide('B10',6.6262e-6) + coolant.add_nuclide('B11',2.6839e-5) + materials.append(coolant) + + absorber = openmc.Material(material_id=3,name="pellet B4C") + absorber.set_density('sum') + absorber.add_nuclide('C0',0.01966) + absorber.add_nuclide('B11',4.7344e-6) + absorber.add_nuclide('B10',1.9177e-5) + materials.append(absorber) + + zirc = openmc.Material(material_id=4,name="Zirc4") + zirc.set_density('sum') + zirc.add_element('Zr',4.23e-2) + materials.append(zirc) + + materials.export_to_xml() + + ### Geometry ### + pin_rad = 0.7 # cm + assembly_pitch = 1.235 # cm + hexagonal_pitch = 23.6 # cm + length = 10.0 # cm + #Fuel pin surfaces + cylfuelin = openmc.ZCylinder(surface_id=1,r=0.386) + cylfuelout = openmc.ZCylinder(surface_id=2,r=0.4582) + #Fuel cells + infcell = openmc.Cell(cell_id=1) + infcell.region = -cylfuelin + infcell.fill = fuel_mat + + clfcell = openmc.Cell(cell_id=2) + clfcell.region = -cylfuelout & +cylfuelin + clfcell.fill = zirc + + outfcell = openmc.Cell(cell_id=3) + outfcell.region = +cylfuelout + outfcell.fill = coolant + #Fuel universe + fuel_ch_univ = openmc.Universe(universe_id=1,name="Fuel channel",cells=[infcell, clfcell,outfcell]) + + #Central tube surfaces + cyltubein = openmc.ZCylinder(surface_id=3,r=0.45) + cyltubeout = openmc.ZCylinder(surface_id=4,r=0.5177) + #Central tube cells + inctcell = openmc.Cell(cell_id=4) + inctcell.region = -cyltubein + inctcell.fill = coolant + + clctcell = openmc.Cell(cell_id=5) + clctcell.region = -cyltubeout & +cyltubein + clctcell.fill = zirc + + outctcell = openmc.Cell(cell_id=6) + outctcell.region = +cyltubeout + outctcell.fill = coolant + #Central tubel universe + tube_ch_univ = openmc.Universe(universe_id=2,name="Central tube channel",cells=[inctcell, clctcell,outctcell]) + + #Absorber tube surfaces + cylabsin = openmc.ZCylinder(surface_id=5,r=0.35) + cylabsout = openmc.ZCylinder(surface_id=6,r=0.41) + cylabsclin = openmc.ZCylinder(surface_id=7,r=0.545) + cylabsclout = openmc.ZCylinder(surface_id=8,r=0.6323) + #Absorber tube cells + inabscell = openmc.Cell(cell_id=7) + inabscell.region = -cylabsin + inabscell.fill = absorber + + clabscell = openmc.Cell(cell_id=8) + clabscell.region = -cylabsout & +cylabsin + clabscell.fill = zirc + + interabscell = openmc.Cell(cell_id=9) + interabscell.region = -cylabsclin & +cylabsout + interabscell.fill = coolant + + clatcell = openmc.Cell(cell_id=10) + clatcell.region = -cylabsclout & +cylabsclin + clatcell.fill = zirc + + outabscell = openmc.Cell(cell_id=11) + outabscell.region = +cylabsclout + outabscell.fill = coolant + #Absorber tube universe + abs_ch_univ = openmc.Universe(universe_id=3,name="Central tube channel",cells=[inabscell,clabscell,interabscell,clatcell,outabscell]) + #Assembly surfaces + edge_length = (1./np.sqrt(3.0)) * hexagonal_pitch + fuel_bottom = openmc.ZPlane(surface_id=9,z0=0.0, boundary_type = 'reflective') + fuel_top = openmc.ZPlane(surface_id=10,z0=length, boundary_type = 'reflective') + # a hex surface for the core to go inside of + hexprism = openmc.model.get_hexagonal_prism(edge_length=edge_length, + origin=(0.0, 0.0), + boundary_type = 'reflective', + orientation='x') + region = hexprism & +fuel_bottom & -fuel_top + + inf_mat = openmc.Cell(cell_id=12) + inf_mat.fill = coolant + inf_mat_univ = openmc.Universe(universe_id=4,cells=[inf_mat,]) + # Fill lattice by channels + nring = 11 + universes = [] # array of rings + for ring in range(nring - 1,-1,-1): + arr = [] + arr.append(fuel_ch_univ) + for cell in range(ring * 6 - 1): + arr.append(fuel_ch_univ) + universes.append(arr) + universes[-1] = [tube_ch_univ] + channels = [(7,16),(7,13),(7,10),(7,7),(7,4),(7,1),(5,0),\ + (4,33),(5,25),(4,27),(5,20),(4,21),(5,15),\ + (4,15),(5,10),(4,9),(5,5),(4,3)] + for i,j in channels: + universes[i][j] = abs_ch_univ + lattice = openmc.HexLattice(name="regular fuel assembly") + lattice.orientation = "ox" + lattice.center = (0., 0., length/2.0) + lattice.pitch = (assembly_pitch,length/2.0) + lattice.universes = 2*[universes] + lattice.outer = inf_mat_univ + + assembly_cell = openmc.Cell(cell_id=13,name="container assembly cell") + assembly_cell.region = region + assembly_cell.fill = lattice + + root_univ = openmc.Universe(universe_id=5,name="root universe", cells=[assembly_cell,]) + + geom = openmc.Geometry(root_univ) + geom.export_to_xml() + + ### Settings ### + + settings = openmc.Settings() + settings.run_mode = 'eigenvalue' + + source = openmc.Source() + ll = [-edge_length, -edge_length, 0.0] + ur = [edge_length, edge_length, 10.0] + source.space = openmc.stats.Box(ll, ur) + source.strength = 1.0 + settings.source = source + settings.output = {'summary' : False} + settings.batches = 10 + settings.inactive = 5 + settings.particles = 1000 + settings.seed = 22 + settings.export_to_xml() + +def test_lattice_hex_ox_surf(): + harness = HexLatticeOXTestHarness('statepoint.10.h5') + harness.main() From 818a6c245b92a749b5b7d1ddd8279a4d34fe9bd9 Mon Sep 17 00:00:00 2001 From: dryuri92 <39188804+dryuri92@users.noreply.github.com> Date: Thu, 30 May 2019 11:02:42 +0300 Subject: [PATCH 12/31] Update constants.h --- include/openmc/constants.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 7abb3d7b9..d9307abe2 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -25,7 +25,7 @@ constexpr int VERSION_RELEASE {0}; constexpr std::array VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE}; // HDF5 data format -constexpr int HDF5_VERSION[] {2, 0};//DR change version HDF5 to work +constexpr int HDF5_VERSION[] {2, 0}; // Version numbers for binary files constexpr std::array VERSION_STATEPOINT {17, 0}; From 86dc85621bd809b58124bb11d079dfe29618e494 Mon Sep 17 00:00:00 2001 From: dryuri92 <39188804+dryuri92@users.noreply.github.com> Date: Thu, 30 May 2019 11:26:49 +0300 Subject: [PATCH 13/31] Update lattice.py with sharper recomendation --- openmc/lattice.py | 79 +++++++++++++++++++++++------------------------ 1 file changed, 39 insertions(+), 40 deletions(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index 45a3c1b81..9931f2bd2 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -136,15 +136,13 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): lattice.universes = uarray elif lattice_type == 'hexagonal': - n_rings = group['n_rings'][()] n_axial = group['n_axial'][()] center = group['center'][()] pitch = group['pitch'][()] outer = group['outer'][()] - #DR Temp added for compatibility with previous version if ('orientation' in group.keys()): - orientation = group['orientation'].value + orientation = group['orientation'][()] else: orientation = "oy" @@ -155,7 +153,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): lattice = openmc.HexLattice(lattice_id, name) lattice.center = center lattice.pitch = pitch - + lattice.orientation = orientation # If the Universe specified outer the Lattice is not void if outer >= 0: lattice.outer = universes[outer] @@ -880,11 +878,11 @@ class HexLattice(Lattice): Most methods for this class use a natural indexing scheme wherein elements are assigned an index corresponding to their position relative to skewed - :math:`(x,\alpha,z)` - OY orientation or :math:`(\alpha,y,z)` - OX orientation //DR - axes as described fully in :ref:`hexagonal_indexing`. However, note that when - universes are assigned to lattice elements using the :attr:`HexLattice.universes` property, - the array indices do not correspond to natural indices. - + :math:`(x,\alpha,z)` or :math:`(\alpha,y,z)` bases, depending on the lattice + orientation, as described fully in :ref:`hexagonal_indexing`. However, note + that when universes are assigned to lattice elements using the + :attr:`HexLattice.universes` property, the array indices do not correspond + to natural indices. Parameters ---------- lattice_id : int, optional @@ -921,7 +919,7 @@ class HexLattice(Lattice): from the outermost ring), and i is the index with a ring starting from the top and proceeding clockwise. orientation : str by default 'OY' orientation of main lattice diagonal another option - - 'OX'//DR + - 'OX' num_rings : int Number of radial ring positions in the xy-plane num_axial : int @@ -970,12 +968,14 @@ class HexLattice(Lattice): @property def num_rings(self): return self._num_rings + @property - def orientation(self):#//DR + def orientation(self): if self._hextype: return "OX" else: return "OY" + @property def num_axial(self): return self._num_axial @@ -1031,7 +1031,7 @@ class HexLattice(Lattice): @orientation.setter def orientation(self, orientation): - cv.check_type('orientation', orientation, str) + cv.check_value('orientation', orientation.lower(), ('ox', 'oy')) if orientation.lower() == "ox": self._hextype = 1 @@ -1125,7 +1125,7 @@ class HexLattice(Lattice): ------- 3-tuple of int Indices of corresponding lattice element in :math:`(x,\alpha,z)` - or :math:`(\alpha,y,z)`basesis + or :math:`(\alpha,y,z)`bases numpy.ndarray Carestian coordinates of the point in the corresponding lattice element coordinate system @@ -1139,27 +1139,20 @@ class HexLattice(Lattice): else: z = point[2] - self.center[2] iz = floor(z/self.pitch[1] + 0.5*self.num_axial) - if self._hextype:#//DR - alpha = y - x*sqrt(3.) - ix = 0 - ia = floor(-alpha/(sqrt(3.0) * self.pitch[0])) - iy = floor(y/(sqrt(0.75) * self.pitch[0])) + if self._hextype: + alpha = y - x*sqrt(3.) + i1 = floor(-alpha/(sqrt(3.0) * self.pitch[0])) + i2 = floor(y/(sqrt(0.75) * self.pitch[0])) else: - alpha = y - x/sqrt(3.) - ix = floor(x/(sqrt(0.75) * self.pitch[0])) - ia = floor(alpha/self.pitch[0]) - iy = 0 - - + alpha = y - x/sqrt(3.) + i1 = floor(x/(sqrt(0.75) * self.pitch[0])) + i2 = floor(alpha/self.pitch[0] # Check four lattice elements to see which one is closest based on local # coordinates - #//DR - numbersOY=[(ix, ia, iz), (ix + 1, ia, iz), (ix, ia + 1, iz), \ - (ix + 1, ia + 1, iz)]# in case OY default orientation - numbersOX=[(ia, iy, iz), (ia + 1, iy, iz), (ia, iy + 1, iz), \ - (ia + 1, iy + 1, iz)]# in case OX orientation + numbers=[(i1, i2, iz), (i1 + 1, i2, iz), (i1, i2 + 1, iz), \ + (i1 + 1, i2 + 1, iz)] d_min = np.inf - for idx in (numbersOX if self._hextype else numbersOY): #//DR + for idx in numbers: p = self.get_local_coordinates(point, idx) d = p[0]**2 + p[1]**2 if d < d_min: @@ -1187,12 +1180,16 @@ class HexLattice(Lattice): system """ - if self._hextype:#//DR - x = point[0] - (self.center[0] + self.pitch[0]*idx[0] + 0.5*self.pitch[0]*idx[1]) - y = point[1] - (self.center[1] + (sqrt(0.75)*self.pitch[0]*idx[1])) + if self._hextype: + x = point[0] - (self.center[0] + self.pitch[0]*idx[0] + + 0.5*self.pitch[0]*idx[1]) + y = point[1] - (self.center[1] + + (sqrt(0.75)*self.pitch[0]*idx[1])) else: - x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0]) - y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0]) + x = point[0] - (self.center[0] + + sqrt(0.75)*self.pitch[0]*idx[0]) + y = point[1] - (self.center[1] + + (0.5*idx[0] + idx[1])*self.pitch[0]) if self._num_axial is None: z = point[2] @@ -1200,7 +1197,8 @@ class HexLattice(Lattice): z = point[2] - (self.center[2] + (idx[2] + 0.5 - 0.5*self.num_axial)* self.pitch[1]) return (x, y, z) - def get_universe_index(self, idx): #//DR + + def get_universe_index(self, idx): r"""Return index in the universes array corresponding to a lattice element index Parameters @@ -1217,10 +1215,11 @@ class HexLattice(Lattice): """ if self._hextype: - return self.get_universe_index_ox(idx) + return self._get_universe_index_ox(idx) else: - return self.get_universe_index_oy(idx) - def get_universe_index_oy(self, idx):#//DR + return self._get_universe_index_oy(idx) + + def _get_universe_index_oy(self, idx): r"""Return index in the universes array corresponding to a lattice element index Parameters @@ -1261,7 +1260,7 @@ class HexLattice(Lattice): return (idx[2], i_ring, i_within) - def get_universe_index_ox(self, idx):#//DR + def _get_universe_index_ox(self, idx): r"""Return index in the universes array corresponding to a lattice element index numeration opposite clockwise From 8ae6c30d6b2880a12125d6e85bed00015ebfc878 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Thu, 30 May 2019 12:29:33 +0300 Subject: [PATCH 14/31] correct python code --- openmc/lattice.py | 8 ++++---- tests/regression_tests/lattice_hex_OX/test.py | 15 ++++++++++----- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index 9931f2bd2..281cc3add 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -1146,13 +1146,13 @@ class HexLattice(Lattice): else: alpha = y - x/sqrt(3.) i1 = floor(x/(sqrt(0.75) * self.pitch[0])) - i2 = floor(alpha/self.pitch[0] + i2 = floor(alpha/self.pitch[0]) # Check four lattice elements to see which one is closest based on local # coordinates - numbers=[(i1, i2, iz), (i1 + 1, i2, iz), (i1, i2 + 1, iz), \ - (i1 + 1, i2 + 1, iz)] + numbersaxy=[(i1, i2, iz), (i1 + 1, i2, iz), (i1, i2 + 1, iz), + (i1 + 1, i2 + 1, iz)] d_min = np.inf - for idx in numbers: + for idx in numbersaxy: p = self.get_local_coordinates(point, idx) d = p[0]**2 + p[1]**2 if d < d_min: diff --git a/tests/regression_tests/lattice_hex_OX/test.py b/tests/regression_tests/lattice_hex_OX/test.py index 2eab458dc..0af36e6f7 100644 --- a/tests/regression_tests/lattice_hex_OX/test.py +++ b/tests/regression_tests/lattice_hex_OX/test.py @@ -56,7 +56,8 @@ class HexLatticeOXTestHarness(PyAPITestHarness): outfcell.region = +cylfuelout outfcell.fill = coolant #Fuel universe - fuel_ch_univ = openmc.Universe(universe_id=1,name="Fuel channel",cells=[infcell, clfcell,outfcell]) + fuel_ch_univ = openmc.Universe(universe_id=1,name="Fuel channel", + cells=[infcell, clfcell,outfcell]) #Central tube surfaces cyltubein = openmc.ZCylinder(surface_id=3,r=0.45) @@ -74,7 +75,8 @@ class HexLatticeOXTestHarness(PyAPITestHarness): outctcell.region = +cyltubeout outctcell.fill = coolant #Central tubel universe - tube_ch_univ = openmc.Universe(universe_id=2,name="Central tube channel",cells=[inctcell, clctcell,outctcell]) + tube_ch_univ = openmc.Universe(universe_id=2,name="Central tube channel", + cells=[inctcell, clctcell,outctcell]) #Absorber tube surfaces cylabsin = openmc.ZCylinder(surface_id=5,r=0.35) @@ -102,11 +104,14 @@ class HexLatticeOXTestHarness(PyAPITestHarness): outabscell.region = +cylabsclout outabscell.fill = coolant #Absorber tube universe - abs_ch_univ = openmc.Universe(universe_id=3,name="Central tube channel",cells=[inabscell,clabscell,interabscell,clatcell,outabscell]) + abs_ch_univ = openmc.Universe(universe_id=3,name="Central tube channel",\ + cells=[inabscell,clabscell,interabscell,clatcell,outabscell]) #Assembly surfaces edge_length = (1./np.sqrt(3.0)) * hexagonal_pitch - fuel_bottom = openmc.ZPlane(surface_id=9,z0=0.0, boundary_type = 'reflective') - fuel_top = openmc.ZPlane(surface_id=10,z0=length, boundary_type = 'reflective') + fuel_bottom = openmc.ZPlane(surface_id=9,z0=0.0, + boundary_type = 'reflective') + fuel_top = openmc.ZPlane(surface_id=10,z0=length, + boundary_type = 'reflective') # a hex surface for the core to go inside of hexprism = openmc.model.get_hexagonal_prism(edge_length=edge_length, origin=(0.0, 0.0), From 00f346d40cb6f4602756e90a62f278adca4891fd Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Thu, 30 May 2019 13:55:32 +0300 Subject: [PATCH 15/31] add a decoding method to transform bytes to string --- openmc/lattice.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index 281cc3add..429d8b691 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -142,7 +142,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): pitch = group['pitch'][()] outer = group['outer'][()] if ('orientation' in group.keys()): - orientation = group['orientation'][()] + orientation = group['orientation'][()].decode() else: orientation = "oy" From fb7a895ce9fa5212f0ef9be161241d3571b92687 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Thu, 30 May 2019 14:13:27 +0300 Subject: [PATCH 16/31] removed unused marks --- openmc/lattice.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index 429d8b691..cec9a6725 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -940,7 +940,7 @@ class HexLattice(Lattice): string = 'HexLattice\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) - string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', "OX" if self._hextype else "OY" )#//DR + string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', "OX" if self._hextype else "OY" ) string += '{0: <16}{1}{2}\n'.format('\t# Rings', '=\t', self._num_rings) string += '{0: <16}{1}{2}\n'.format('\t# Axial', '=\t', self._num_axial) string += '{0: <16}{1}{2}\n'.format('\tCenter', '=\t', @@ -1802,7 +1802,7 @@ class HexLattice(Lattice): # Join the rows together and return the string. return '\n'.join(rows) @staticmethod - def show_indices_ox(num_rings):#//DR + def show_indices_ox(num_rings): """Return a diagram of the hexagonal lattice with OX orientation layout with indices. This method can be used to show the proper indices to be used when From 4e0c52a62826770325e9484b09f956fd262e0ce8 Mon Sep 17 00:00:00 2001 From: dryuri92 <39188804+dryuri92@users.noreply.github.com> Date: Thu, 30 May 2019 22:06:36 +0300 Subject: [PATCH 17/31] rename HexOrientation to Orientation --- include/openmc/lattice.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index 6b30486b8..b79e77660 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -27,9 +27,6 @@ enum class LatticeType { rect, hex }; -enum class HexOrientation { - oy, ox -}; //============================================================================== // Global variables @@ -269,6 +266,9 @@ public: void to_hdf5_inner(hid_t group_id) const; private: + enum class Orientation { + y, x + }; //! Fill universes_ vector for OY orientation void fill_lattice_oy(std::vector univ_words); @@ -277,7 +277,7 @@ private: int n_rings_; //!< Number of radial tile positions int n_axial_; //!< Number of axial tile positions - HexOrientation orientation_; //!< Flat side up (oy) vs. sharp side up (ox) + Orientation orientation_; //!< Flat side up (oy) vs. sharp side up (ox) Position center_; //!< Global center of lattice std::array pitch_; //!< Lattice tile width and height }; From 0333c5d500d5e205ad48aafdc2a448be91717c6e Mon Sep 17 00:00:00 2001 From: dryuri92 <39188804+dryuri92@users.noreply.github.com> Date: Thu, 30 May 2019 22:10:07 +0300 Subject: [PATCH 18/31] call function by ref --- include/openmc/lattice.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index b79e77660..195f62785 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -270,14 +270,14 @@ private: y, x }; //! Fill universes_ vector for OY orientation - void fill_lattice_oy(std::vector univ_words); + void fill_lattice_oy(const std::vector& univ_words); //! Fill universes_ vector for OX orientation - void fill_lattice_ox(std::vector univ_words); + void fill_lattice_ox(const std::vector& univ_words); int n_rings_; //!< Number of radial tile positions int n_axial_; //!< Number of axial tile positions - Orientation orientation_; //!< Flat side up (oy) vs. sharp side up (ox) + Orientation orientation_; //!< Flat side up (oy) vs. sharp side up (ox) Position center_; //!< Global center of lattice std::array pitch_; //!< Lattice tile width and height }; From 40b5c23aabdf9bb8968d1d6afb846d908d71199e Mon Sep 17 00:00:00 2001 From: dryuri92 <39188804+dryuri92@users.noreply.github.com> Date: Thu, 30 May 2019 22:35:35 +0300 Subject: [PATCH 19/31] add call fucntion by reference --- src/lattice.cpp | 34 +++++++++++++++------------------- 1 file changed, 15 insertions(+), 19 deletions(-) diff --git a/src/lattice.cpp b/src/lattice.cpp index b98ba6366..ec9fa908c 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -444,15 +444,15 @@ HexLattice::HexLattice(pugi::xml_node lat_node) if (check_for_node(lat_node, "orient")) { std::string orientation = get_node_value(lat_node, "orient"); if (orientation == "OY") { - orientation_ = HexOrientation::oy; + orientation_ = Orientation::oy; } else if (orientation == "OX") { - orientation_ = HexOrientation::ox; + orientation_ = Orientation::ox; } else { fatal_error("Unrecognized orientation '" + orientation + "' for lattice " + std::to_string(id_)); } } else { - orientation_ = HexOrientation::oy; + orientation_ = Orientation::oy; } // Read the lattice center. @@ -504,22 +504,18 @@ HexLattice::HexLattice(pugi::xml_node lat_node) // the following code walks a set of index values across the skewed array // in a manner that matches the input order. Note that i_x = 0, i_a = 0 // or i_a = 0, i_y = 0 corresponds to the center of the hexagonal lattice. - universes_.resize((2*n_rings_-1) * (2*n_rings_-1) * n_axial_, C_NONE); - - if (orientation_ == HexOrientation::oy) { + if (orientation_ == Orientation::oy) { fill_lattice_oy(univ_words); } else { fill_lattice_ox(univ_words); } - - } - + //============================================================================== - + void -HexLattice::fill_lattice_ox(std::vector univ_words) +HexLattice::fill_lattice_ox(const std::vector& univ_words) { int input_index = 0; for (int m = 0; m < n_axial_; m++) { @@ -572,7 +568,7 @@ HexLattice::fill_lattice_ox(std::vector univ_words) //============================================================================== void -HexLattice::fill_lattice_oy(std::vector univ_words) +HexLattice::fill_lattice_oy(const std::vector& univ_words) { int input_index = 0; for (int m = 0; m < n_axial_; m++) { @@ -712,7 +708,7 @@ const double beta_dir; double gamma_dir; double delta_dir; - if (orientation_ == HexOrientation::oy) { + if (orientation_ == Orientation::oy) { beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0; delta_dir = u.y; @@ -741,7 +737,7 @@ const r_t = get_local_position(r, i_xyz_t); } double beta; - if (orientation_ == HexOrientation::oy) { + if (orientation_ == Orientation::oy) { beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0; } else { beta = r_t.x; @@ -765,7 +761,7 @@ const r_t = get_local_position(r, i_xyz_t); } double gamma; - if (orientation_ == HexOrientation::oy) { + if (orientation_ == Orientation::oy) { gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0; } else { gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0; @@ -792,7 +788,7 @@ const r_t = get_local_position(r, i_xyz_t); } double delta; - if (orientation_ == HexOrientation::oy) { + if (orientation_ == Orientation::oy) { delta = r_t.y; } else { delta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0; @@ -852,7 +848,7 @@ HexLattice::get_indices(Position r, Direction u) const } int i1, i2; - if (orientation_ == HexOrientation::oy) { + if (orientation_ == Orientation::oy) { // Convert coordinates into skewed bases. The (x, alpha) basis is used to // find the index of the global coordinates to within 4 cells. double alpha = r_o.y - r_o.x / std::sqrt(3.0); @@ -930,7 +926,7 @@ Position HexLattice::get_local_position(Position r, const std::array i_xyz) const { - if (orientation_ == HexOrientation::oy) { + if (orientation_ == Orientation::oy) { // x_l = x_g - (center + pitch_x*cos(30)*index_x) r.x -= center_.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0]; @@ -1010,7 +1006,7 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const write_string(lat_group, "type", "hexagonal", false); write_dataset(lat_group, "n_rings", n_rings_); write_dataset(lat_group, "n_axial", n_axial_); - if (orientation_ == HexOrientation::oy) { + if (orientation_ == Orientation::oy) { write_string(lat_group, "orientation", "oy", false); } else { write_string(lat_group, "orientation", "ox", false); From feac9c83a8ff251a96652ec93cb90b6679beca9b Mon Sep 17 00:00:00 2001 From: dryuri92 <39188804+dryuri92@users.noreply.github.com> Date: Thu, 30 May 2019 22:38:54 +0300 Subject: [PATCH 20/31] is3d separately --- src/lattice.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/lattice.cpp b/src/lattice.cpp index ec9fa908c..d2e29e106 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -933,10 +933,6 @@ const // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0] + (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0); - if (is_3d_) { - r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; - } - } else { // x_l = x_g - (center + pitch_x*index_a + pitch_y*sin(30)*index_y) r.x -= (center_.x + (i_xyz[0] - n_rings_ + 1) * pitch_[0] @@ -944,10 +940,11 @@ const // y_l = y_g - (center + pitch_y*cos(30)*index_y) r.y -= center_.y + std::sqrt(3.0)/2.0 * (i_xyz[1] - n_rings_ + 1) * pitch_[0]; - if (is_3d_) { + } + + if (is_3d_) { r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; } - } return r; } From 40e8f3c22cfbb830c42470a43ca651d584ec7271 Mon Sep 17 00:00:00 2001 From: dryuri92 <39188804+dryuri92@users.noreply.github.com> Date: Thu, 30 May 2019 22:40:09 +0300 Subject: [PATCH 21/31] simple change --- src/initialize.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/initialize.cpp b/src/initialize.cpp index 58d96bf40..061025745 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -36,11 +36,9 @@ int openmc_init(int argc, char* argv[], const void* intracomm) { using namespace openmc; - - + char cCurrentPath[256]; - - + #ifdef OPENMC_MPI // Check if intracomm was passed MPI_Comm comm; From de0ffecdd517e6ca77de3d2a1ee0bcb6d2843857 Mon Sep 17 00:00:00 2001 From: dryuri92 <39188804+dryuri92@users.noreply.github.com> Date: Fri, 31 May 2019 00:59:54 +0300 Subject: [PATCH 22/31] Accept changes by P.Romano --- openmc/lattice.py | 179 ++++++++++++++++++---------------------------- 1 file changed, 70 insertions(+), 109 deletions(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index cec9a6725..8d2d82ed0 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -141,7 +141,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): center = group['center'][()] pitch = group['pitch'][()] outer = group['outer'][()] - if ('orientation' in group.keys()): + if ('orientation' in group): orientation = group['orientation'][()].decode() else: orientation = "oy" @@ -157,7 +157,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): # If the Universe specified outer the Lattice is not void if outer >= 0: lattice.outer = universes[outer] - if (orientation.lower() is "oy"): + if (orientation.lower() == "oy"): # Build array of Universe pointers for the Lattice. Note that # we need to convert between the HDF5's square array of # (x, alpha, z) to the Python API's format of a ragged nested @@ -883,6 +883,7 @@ class HexLattice(Lattice): that when universes are assigned to lattice elements using the :attr:`HexLattice.universes` property, the array indices do not correspond to natural indices. + Parameters ---------- lattice_id : int, optional @@ -918,8 +919,9 @@ class HexLattice(Lattice): possible, where z is the axial index, r is in the ring index (starting from the outermost ring), and i is the index with a ring starting from the top and proceeding clockwise. - orientation : str by default 'OY' orientation of main lattice diagonal another option - - 'OX' + orientation : {'x', 'y'} + str by default 'y' orientation of main lattice diagonal another option + - 'x' num_rings : int Number of radial ring positions in the xy-plane num_axial : int @@ -934,13 +936,13 @@ class HexLattice(Lattice): self._num_rings = None self._num_axial = None self._center = None - self._hextype = 0 # by Default orientaion is OY + self._orientation = 'y' def __repr__(self): string = 'HexLattice\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) - string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', "OX" if self._hextype else "OY" ) + string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', "OX" if (self._orientation== 'x') else "OY" ) string += '{0: <16}{1}{2}\n'.format('\t# Rings', '=\t', self._num_rings) string += '{0: <16}{1}{2}\n'.format('\t# Axial', '=\t', self._num_axial) string += '{0: <16}{1}{2}\n'.format('\tCenter', '=\t', @@ -971,10 +973,7 @@ class HexLattice(Lattice): @property def orientation(self): - if self._hextype: - return "OX" - else: - return "OY" + return self._orientation @property def num_axial(self): @@ -1032,8 +1031,7 @@ class HexLattice(Lattice): @orientation.setter def orientation(self, orientation): cv.check_value('orientation', orientation.lower(), ('ox', 'oy')) - if orientation.lower() == "ox": - self._hextype = 1 + self._orientation = orientation.lower() @Lattice.pitch.setter def pitch(self, pitch): @@ -1125,7 +1123,7 @@ class HexLattice(Lattice): ------- 3-tuple of int Indices of corresponding lattice element in :math:`(x,\alpha,z)` - or :math:`(\alpha,y,z)`bases + or :math:`(\alpha,y,z)` bases numpy.ndarray Carestian coordinates of the point in the corresponding lattice element coordinate system @@ -1139,7 +1137,7 @@ class HexLattice(Lattice): else: z = point[2] - self.center[2] iz = floor(z/self.pitch[1] + 0.5*self.num_axial) - if self._hextype: + if self._orientation == 'x': alpha = y - x*sqrt(3.) i1 = floor(-alpha/(sqrt(3.0) * self.pitch[0])) i2 = floor(y/(sqrt(0.75) * self.pitch[0])) @@ -1149,10 +1147,10 @@ class HexLattice(Lattice): i2 = floor(alpha/self.pitch[0]) # Check four lattice elements to see which one is closest based on local # coordinates - numbersaxy=[(i1, i2, iz), (i1 + 1, i2, iz), (i1, i2 + 1, iz), - (i1 + 1, i2 + 1, iz)] + indices = [(i1, i2, iz), (i1 + 1, i2, iz), (i1, i2 + 1, iz), + (i1 + 1, i2 + 1, iz)] d_min = np.inf - for idx in numbersaxy: + for idx in indices: p = self.get_local_coordinates(point, idx) d = p[0]**2 + p[1]**2 if d < d_min: @@ -1171,7 +1169,7 @@ class HexLattice(Lattice): Cartesian coordinates of point idx : Iterable of int Indices of lattice element in :math:`(x,\alpha,z)` - or :math:`(\alpha,y,z)` basesis + or :math:`(\alpha,y,z)` bases Returns ------- @@ -1180,17 +1178,17 @@ class HexLattice(Lattice): system """ - if self._hextype: + if self._orientation == 'x': x = point[0] - (self.center[0] + self.pitch[0]*idx[0] + 0.5*self.pitch[0]*idx[1]) y = point[1] - (self.center[1] + - (sqrt(0.75)*self.pitch[0]*idx[1])) + sqrt(0.75)*self.pitch[0]*idx[1]) else: x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0]) y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0]) - + if self._num_axial is None: z = point[2] else: @@ -1210,31 +1208,11 @@ class HexLattice(Lattice): Returns ------- - 2- or 3-tuple of int + 2- or 3-tuple of int Indices used when setting the :attr:`HexLattice.universes` property """ - if self._hextype: - return self._get_universe_index_ox(idx) - else: - return self._get_universe_index_oy(idx) - - def _get_universe_index_oy(self, idx): - r"""Return index in the universes array corresponding to a lattice element index - - Parameters - ---------- - idx : Iterable of int - Lattice element indices in the :math:`(x,\alpha,z)` coordinate - system - - Returns - ------- - 2- or 3-tuple of int - Indices used when setting the :attr:`HexLattice.universes` property - - """ - + # First we determine which ring the index corresponds to. x = idx[0] a = idx[1] @@ -1253,54 +1231,15 @@ class HexLattice(Lattice): i_within = 3*g - x else: i_within = 5*g - z - + + if (self._orientation == 'x') and (i_within > 0): + i_within = 6*(self._num_rings - i_ring - 1) - i_within_ + if self.num_axial is None: return (i_ring, i_within) else: return (idx[2], i_ring, i_within) - - - def _get_universe_index_ox(self, idx): - r"""Return index in the universes array corresponding to a lattice element index - numeration opposite clockwise - - Parameters - ---------- - idx : Iterable of int - Lattice element indices in the :math:`(\alpha,y,z)` coordinate - system - - Returns - ------- - 2- or 3-tuple of int - Indices used when setting the :attr:`HexLattice.universes` property - - """ - - # First we determine which ring the index corresponds to. - a = idx[0] - y = idx[1] - z = -a - y - g = max(abs(y), abs(a), abs(z)) - - # Next we use a clever method to figure out where along the ring we are. - i_ring = self._num_rings - 1 - g - if y >= 0: - if a >= 0: - i_within = y - else: - i_within = 2*g + z - else: - if a <= 0: - i_within = 3*g - y - else: - i_within = 5*g - z - - if self.num_axial is None: - return (i_ring, i_within) - else: - return (idx[2], i_ring, i_within) - + def is_valid_index(self, idx): r"""Determine whether lattice element index is within defined range @@ -1352,7 +1291,7 @@ class HexLattice(Lattice): lattice_subelement.set("n_rings", str(self._num_rings)) # If orientation is "OX" export it to XML - if self._hextype: + if self._orientation== 'x': lattice_subelement.set("orient", "OX") if self._num_axial is not None: @@ -1483,10 +1422,10 @@ class HexLattice(Lattice): each sub-list represents a single ring. The first list should be the outer ring. """ - if self._hextype: - return self._repr_axial_slice_ox(universes) + if self._orientation== 'x': + return self._repr_axial_slice_ox(universes) else: - return self._repr_axial_slice_oy(universes) + return self._repr_axial_slice_oy(universes) def _repr_axial_slice_ox(self, universes): """Return string representation for the given 2D group of universes in 'OX' orientation case. @@ -1533,7 +1472,7 @@ class HexLattice(Lattice): for i in range(r): # Add the universe. universe = universes[r_prime][theta] - rows[y].insert(0,id_form.format(universe._id)) + rows[y].insert(0, id_form.format(universe._id)) # Translate the indices. theta += 1 @@ -1542,7 +1481,7 @@ class HexLattice(Lattice): for i in range(r): # Add the universe. universe = universes[r_prime][theta] - rows[y].insert(0,id_form.format(universe._id)) + rows[y].insert(0, id_form.format(universe._id)) # Translate the indices. y -= 1 @@ -1584,12 +1523,11 @@ class HexLattice(Lattice): for y in range(self._num_rings - 1): rows[y] = (self._num_rings - 1 - y)*pad + rows[y] rows[-1 - y] = (self._num_rings - 1 - y)*pad + rows[-1 - y] - + # Join the rows together and return the string. universe_ids = '\n'.join(rows) return universe_ids - - + def _repr_axial_slice_oy(self, universes): """Return string representation for the given 2D group of universes in 'OY' orientation case.. @@ -1697,9 +1635,8 @@ class HexLattice(Lattice): # Join the rows together and return the string. universe_ids = '\n'.join(rows) return universe_ids - - @staticmethod - def show_indices(num_rings): + + def _show_indices_y(num_rings): """Return a diagram of the hexagonal lattice layout with indices. This method can be used to show the proper indices to be used when @@ -1801,24 +1738,24 @@ class HexLattice(Lattice): # Join the rows together and return the string. return '\n'.join(rows) - @staticmethod - def show_indices_ox(num_rings): + + def _show_indices_x(num_rings): """Return a diagram of the hexagonal lattice with OX orientation layout with indices. This method can be used to show the proper indices to be used when setting the :attr:`HexLattice.universes` property. For example, running this method with num_rings=3 will return the similar diagram:: - (0, 4) (0, 3) (0, 2) + (0, 4) (0, 3) (0, 2) + + (0, 5) (1, 2) (1, 1) (0, 1) + + (0, 6) (1, 3) (2, 0) (1, 0) (0, 0) - (0, 5) (1, 2) (1, 1) (0, 1) - - (0, 6) (1, 3) (2, 0) (1, 0) (0, 0) - - (0, 7) (1, 4) (1, 5) (0,11) - - (0, 8) (0, 9) (0,10) + (0, 7) (1, 4) (1, 5) (0,11) + (0, 8) (0, 9) (0,10) + Parameters ---------- num_rings : int @@ -1897,4 +1834,28 @@ class HexLattice(Lattice): # Join the rows together and return the string. return '\n\n'.join(rows) + + @staticmethod + def show_indices(num_rings, orientation="y"): + """Return a diagram of the hexagonal lattice layout with indices. + + Parameters + ---------- + num_rings : int + Number of rings in the hexagonal lattice + orientation : {"x", "y"} + str type of hex lattice orientation + + Returns + ------- + str + Diagram of the hexagonal lattice showing indices + + """ + + if self._orientation == 'x': + return self._show_indices_x(num_rings) + else: + return self._show_indices_y(num_rings) + From d86026689254fb1f36e1f80a9d5eaa13c4fac87b Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Fri, 31 May 2019 02:38:06 +0300 Subject: [PATCH 23/31] test lattice OX according to PEP8 --- include/openmc/lattice.h | 4 +- openmc/lattice.py | 309 +++++++++--------- tests/regression_tests/lattice_hex_OX/test.py | 157 +++++---- 3 files changed, 255 insertions(+), 215 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index 195f62785..c06920759 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -266,8 +266,8 @@ public: void to_hdf5_inner(hid_t group_id) const; private: - enum class Orientation { - y, x + enum class Orientation { + oy, ox }; //! Fill universes_ vector for OY orientation void fill_lattice_oy(const std::vector& univ_words); diff --git a/openmc/lattice.py b/openmc/lattice.py index 8d2d82ed0..fed12d192 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -36,7 +36,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): outer : openmc.Universe A universe to fill all space outside the lattice universes : Iterable of Iterable of openmc.Universe - A two- or three-dimensional list/array of universes filling each element + A two-or three-dimensional list/array of universes filling each element of the lattice """ @@ -142,10 +142,9 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): pitch = group['pitch'][()] outer = group['outer'][()] if ('orientation' in group): - orientation = group['orientation'][()].decode() + orientation = group['orientation'][()].decode() else: - orientation = "oy" - + orientation = "y" universe_ids = group['universes'][()] @@ -158,120 +157,119 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): if outer >= 0: lattice.outer = universes[outer] if (orientation.lower() == "oy"): - # Build array of Universe pointers for the Lattice. Note that - # we need to convert between the HDF5's square array of - # (x, alpha, z) to the Python API's format of a ragged nested - # list of (z, ring, theta). - uarray = [] - for z in range(n_axial): - # Add a list for this axial level. - uarray.append([]) - x = n_rings - 1 - a = 2*n_rings - 2 - for r in range(n_rings - 1, 0, -1): - # Add a list for this ring. - uarray[-1].append([]) + # Build array of Universe pointers for the Lattice. Note that + # we need to convert between the HDF5's square array of + # (x, alpha, z) to the Python API's format of a ragged nested + # list of (z, ring, theta). + uarray = [] + for z in range(n_axial): + # Add a list for this axial level. + uarray.append([]) + x = n_rings - 1 + a = 2*n_rings - 2 + for r in range(n_rings - 1, 0, -1): + # Add a list for this ring. + uarray[-1].append([]) - # Climb down the top-right. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - x += 1 - a -= 1 + # Climb down the top-right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + x += 1 + a -= 1 - # Climb down the right. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - a -= 1 + # Climb down the right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + a -= 1 - # Climb down the bottom-right. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - x -= 1 + # Climb down the bottom-right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + x -= 1 - # Climb up the bottom-left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - x -= 1 - a += 1 + # Climb up the bottom-left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + x -= 1 + a += 1 - # Climb up the left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - a += 1 + # Climb up the left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + a += 1 - # Climb up the top-left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, a, x]) - x += 1 + # Climb up the top-left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, a, x]) + x += 1 - # Move down to the next ring. - a -= 1 + # Move down to the next ring. + a -= 1 - # Convert the ids into Universe objects. - uarray[-1][-1] = [universes[u_id] - for u_id in uarray[-1][-1]] + # Convert the ids into Universe objects. + uarray[-1][-1] = [universes[u_id] + for u_id in uarray[-1][-1]] - # Handle the degenerate center ring separately. - u_id = universe_ids[z, a, x] - uarray[-1].append([universes[u_id]]) + # Handle the degenerate center ring separately. + u_id = universe_ids[z, a, x] + uarray[-1].append([universes[u_id]]) else: - # Build array of Universe pointers for the Lattice. Note that - # we need to convert between the HDF5's square array of - # (alpha, y, z) to the Python API's format of a ragged nested - # list of (z, ring, theta). - uarray = [] - for z in range(n_axial): - # Add a list for this axial level. - uarray.append([]) - a = 2*n_rings - 2 - y = n_rings - 1 - for r in range(n_rings - 1, 0, -1): - # Add a list for this ring. - uarray[-1].append([]) + # Build array of Universe pointers for the Lattice. Note that + # we need to convert between the HDF5's square array of + # (alpha, y, z) to the Python API's format of a ragged nested + # list of (z, ring, theta). + uarray = [] + for z in range(n_axial): + # Add a list for this axial level. + uarray.append([]) + a = 2*n_rings - 2 + y = n_rings - 1 + for r in range(n_rings - 1, 0, -1): + # Add a list for this ring. + uarray[-1].append([]) - # Climb up the top-left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, y, a]) - a -= 1 - y += 1 - + # Climb up the top-left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a -= 1 + y += 1 - # Climb the left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, y, a]) - a -= 1 + # Climb the left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a -= 1 - # Climb down the bottom-left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, y, a]) - y -= 1 + # Climb down the bottom-left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + y -= 1 - # Climb down the bottom-right. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, y, a]) - a += 1 - y -= 1 + # Climb down the bottom-right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a += 1 + y -= 1 - # Climb up the right. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, y, a]) - a += 1 + # Climb up the right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a += 1 - # Climb up the top-right. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, y, a]) - y += 1 + # Climb up the top-right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + y += 1 - # Move down to the next ring. - a -= 1 + # Move down to the next ring. + a -= 1 - # Convert the ids into Universe objects. - uarray[-1][-1] = [universes[u_id] - for u_id in uarray[-1][-1]] + # Convert the ids into Universe objects. + uarray[-1][-1] = [universes[u_id] + for u_id in uarray[-1][-1]] - # Handle the degenerate center ring separately. - u_id = universe_ids[z, y, a] - uarray[-1].append([universes[u_id]]) + # Handle the degenerate center ring separately. + u_id = universe_ids[z, y, a] + uarray[-1].append([universes[u_id]]) # Add the universes to the lattice. if len(pitch) == 2: @@ -407,7 +405,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): ---------- idx : Iterable of int Lattice element indices. For a rectangular lattice, the indices are - given in the :math:`(x,y)` or :math:`(x,y,z)` coordinate system. For + given in the :math:`(x,y)` or :math:`(x,y,z)` coordinate system.For hexagonal lattices, they are given in the :math:`x,\alpha` or :math:`x,\alpha,z` coordinate systems. @@ -709,7 +707,8 @@ class RectLattice(Lattice): return (x, y, z) def get_universe_index(self, idx): - """Return index in the universes array corresponding to a lattice element index + """Return index in the universes array corresponding + to a lattice element index Parameters ---------- @@ -853,7 +852,8 @@ class RectLattice(Lattice): lat_id = int(get_text(elem, 'id')) name = get_text(elem, 'name') lat = cls(lat_id, name) - lat.lower_left = [float(i) for i in get_text(elem, 'lower_left').split()] + lat.lower_left = [float(i) + for i in get_text(elem, 'lower_left').split()] lat.pitch = [float(i) for i in get_text(elem, 'pitch').split()] outer = get_text(elem, 'outer') if outer is not None: @@ -863,7 +863,7 @@ class RectLattice(Lattice): dimension = get_text(elem, 'dimension').split() shape = np.array(dimension, dtype=int)[::-1] uarray = np.array([get_universe(int(i)) for i in - get_text(elem, 'universes').split()]) + get_text(elem, 'universes').split()]) uarray.shape = shape lat.universes = uarray return lat @@ -883,7 +883,7 @@ class HexLattice(Lattice): that when universes are assigned to lattice elements using the :attr:`HexLattice.universes` property, the array indices do not correspond to natural indices. - + Parameters ---------- lattice_id : int, optional @@ -942,7 +942,9 @@ class HexLattice(Lattice): string = 'HexLattice\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) - string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', "OX" if (self._orientation== 'x') else "OY" ) + string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', + "OX" if (self._orientation == 'x') + else "OY") string += '{0: <16}{1}{2}\n'.format('\t# Rings', '=\t', self._num_rings) string += '{0: <16}{1}{2}\n'.format('\t# Axial', '=\t', self._num_axial) string += '{0: <16}{1}{2}\n'.format('\tCenter', '=\t', @@ -970,11 +972,11 @@ class HexLattice(Lattice): @property def num_rings(self): return self._num_rings - + @property def orientation(self): return self._orientation - + @property def num_axial(self): return self._num_axial @@ -1031,7 +1033,7 @@ class HexLattice(Lattice): @orientation.setter def orientation(self, orientation): cv.check_value('orientation', orientation.lower(), ('ox', 'oy')) - self._orientation = orientation.lower() + self._orientation = orientation.lower() @Lattice.pitch.setter def pitch(self, pitch): @@ -1078,7 +1080,7 @@ class HexLattice(Lattice): # Check the center ring. if len(axial_slice[-1]) != 1: msg = 'HexLattice ID={0:d} has the wrong number of ' \ - 'elements in the innermost ring. Only 1 element is ' \ + 'elements in the innermost ring.Only 1 element is ' \ 'allowed in the innermost ring.'.format(self._id) raise ValueError(msg) @@ -1086,8 +1088,8 @@ class HexLattice(Lattice): for r in range(self._num_rings-1): if len(axial_slice[r]) != 6*(self._num_rings - 1 - r): msg = 'HexLattice ID={0:d} has the wrong number of ' \ - 'elements in ring number {1:d} (counting from the '\ - 'outermost ring). This ring should have {2:d} ' \ + 'elements in ring number{1:d}(counting from the '\ + 'outermost ring). This ring should have {2:d} ' \ 'elements.'.format(self._id, r, 6*(self._num_rings - 1 - r)) raise ValueError(msg) @@ -1148,9 +1150,10 @@ class HexLattice(Lattice): # Check four lattice elements to see which one is closest based on local # coordinates indices = [(i1, i2, iz), (i1 + 1, i2, iz), (i1, i2 + 1, iz), - (i1 + 1, i2 + 1, iz)] + (i1 + 1, i2 + 1, iz)] d_min = np.inf - for idx in indices: + + for idx in indices: p = self.get_local_coordinates(point, idx) d = p[0]**2 + p[1]**2 if d < d_min: @@ -1168,7 +1171,7 @@ class HexLattice(Lattice): point : Iterable of float Cartesian coordinates of point idx : Iterable of int - Indices of lattice element in :math:`(x,\alpha,z)` + Indices of lattice element in :math:`(x,\alpha,z)` or :math:`(\alpha,y,z)` bases Returns @@ -1179,40 +1182,42 @@ class HexLattice(Lattice): """ if self._orientation == 'x': - x = point[0] - (self.center[0] + self.pitch[0]*idx[0] + + x = point[0] - (self.center[0] + self.pitch[0]*idx[0] + 0.5*self.pitch[0]*idx[1]) - y = point[1] - (self.center[1] + + y = point[1] - (self.center[1] + sqrt(0.75)*self.pitch[0]*idx[1]) else: - x = point[0] - (self.center[0] + x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0]) - y = point[1] - (self.center[1] + y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0]) - + if self._num_axial is None: z = point[2] else: z = point[2] - (self.center[2] + (idx[2] + 0.5 - 0.5*self.num_axial)* self.pitch[1]) return (x, y, z) - + def get_universe_index(self, idx): - r"""Return index in the universes array corresponding to a lattice element index + r"""Return index in the universes array corresponding + to a lattice element index Parameters ---------- idx : Iterable of int Lattice element indices in the :math:`(x,\alpha,z)` coordinate - system in 'OY' orientation case, or indices in the :math:`(\alpha,y,z)` coordinate - system in 'OX' one + system in 'OY' orientation case, or indices in the + :math:`(\alpha,y,z)` coordinate system in 'OX' one Returns ------- 2- or 3-tuple of int - Indices used when setting the :attr:`HexLattice.universes` property + Indices used when setting the :attr:`HexLattice.universes` + property """ - + # First we determine which ring the index corresponds to. x = idx[0] a = idx[1] @@ -1231,23 +1236,23 @@ class HexLattice(Lattice): i_within = 3*g - x else: i_within = 5*g - z - + if (self._orientation == 'x') and (i_within > 0): i_within = 6*(self._num_rings - i_ring - 1) - i_within_ - + if self.num_axial is None: return (i_ring, i_within) else: return (idx[2], i_ring, i_within) - + def is_valid_index(self, idx): r"""Determine whether lattice element index is within defined range Parameters ---------- idx : Iterable of int - Lattice element indices in the both :math:`(x,\alpha,z)` and :math:`(\alpha,y,z)` coordinate - system + Lattice element indices in the both :math:`(x,\alpha,z)` + and :math:`(\alpha,y,z)` coordinate system Returns ------- @@ -1291,7 +1296,7 @@ class HexLattice(Lattice): lattice_subelement.set("n_rings", str(self._num_rings)) # If orientation is "OX" export it to XML - if self._orientation== 'x': + if self._orientation == 'x': lattice_subelement.set("orient", "OX") if self._num_axial is not None: @@ -1377,13 +1382,13 @@ class HexLattice(Lattice): # Create empty nested lists for one axial level univs = [[None for _ in range(max(6*(n_rings - 1 - r), 1))] - for r in range(n_rings)] + for r in range(n_rings)] if n_axial > 1: univs = [deepcopy(univs) for i in range(n_axial)] # Get flat array of universes numbers uarray = np.array([get_universe(int(i)) for i in - get_text(elem, 'universes').split()]) + get_text(elem, 'universes').split()]) # Fill nested lists j = 0 @@ -1415,6 +1420,7 @@ class HexLattice(Lattice): j += 1 lat.universes = univs return lat + def _repr_axial_slice(self, universes): """Return string representation for the given 2D group of universes. @@ -1422,13 +1428,14 @@ class HexLattice(Lattice): each sub-list represents a single ring. The first list should be the outer ring. """ - if self._orientation== 'x': + if self._orientation == 'x': return self._repr_axial_slice_ox(universes) else: return self._repr_axial_slice_oy(universes) def _repr_axial_slice_ox(self, universes): - """Return string representation for the given 2D group of universes in 'OX' orientation case. + """Return string representation for the given 2D group of universes + in 'OX' orientation case. The 'universes' argument should be a list of lists of universes where each sub-list represents a single ring. The first list should be the @@ -1523,13 +1530,14 @@ class HexLattice(Lattice): for y in range(self._num_rings - 1): rows[y] = (self._num_rings - 1 - y)*pad + rows[y] rows[-1 - y] = (self._num_rings - 1 - y)*pad + rows[-1 - y] - + # Join the rows together and return the string. universe_ids = '\n'.join(rows) return universe_ids - + def _repr_axial_slice_oy(self, universes): - """Return string representation for the given 2D group of universes in 'OY' orientation case.. + """Return string representation for the given 2D group of universes in + 'OY' orientation case.. The 'universes' argument should be a list of lists of universes where each sub-list represents a single ring. The first list should be the @@ -1635,7 +1643,7 @@ class HexLattice(Lattice): # Join the rows together and return the string. universe_ids = '\n'.join(rows) return universe_ids - + def _show_indices_y(num_rings): """Return a diagram of the hexagonal lattice layout with indices. @@ -1738,24 +1746,25 @@ class HexLattice(Lattice): # Join the rows together and return the string. return '\n'.join(rows) - + def _show_indices_x(num_rings): - """Return a diagram of the hexagonal lattice with OX orientation layout with indices. + """Return a diagram of the hexagonal lattice with OX orientation + layout with indices. This method can be used to show the proper indices to be used when - setting the :attr:`HexLattice.universes` property. For example, running + setting the :attr:`HexLattice.universes` property. For example,running this method with num_rings=3 will return the similar diagram:: (0, 4) (0, 3) (0, 2) - + (0, 5) (1, 2) (1, 1) (0, 1) - + (0, 6) (1, 3) (2, 0) (1, 0) (0, 0) (0, 7) (1, 4) (1, 5) (0,11) (0, 8) (0, 9) (0,10) - + Parameters ---------- num_rings : int @@ -1798,12 +1807,12 @@ class HexLattice(Lattice): for i in range(r): # Climb left the top-left. - rows[y].insert(0,str_form.format(r_prime, theta)) + rows[y].insert(0, str_form.format(r_prime, theta)) theta += 1 for i in range(r): # Climb down the middle left. - rows[y].insert(0,str_form.format(r_prime, theta)) + rows[y].insert(0, str_form.format(r_prime, theta)) y -= 1 theta += 1 @@ -1834,11 +1843,11 @@ class HexLattice(Lattice): # Join the rows together and return the string. return '\n\n'.join(rows) - + @staticmethod def show_indices(num_rings, orientation="y"): """Return a diagram of the hexagonal lattice layout with indices. - + Parameters ---------- num_rings : int @@ -1857,5 +1866,3 @@ class HexLattice(Lattice): return self._show_indices_x(num_rings) else: return self._show_indices_y(num_rings) - - diff --git a/tests/regression_tests/lattice_hex_OX/test.py b/tests/regression_tests/lattice_hex_OX/test.py index 0af36e6f7..bcef2357d 100644 --- a/tests/regression_tests/lattice_hex_OX/test.py +++ b/tests/regression_tests/lattice_hex_OX/test.py @@ -2,48 +2,56 @@ from tests.testing_harness import PyAPITestHarness import openmc import numpy as np + class HexLatticeOXTestHarness(PyAPITestHarness): + def _build_inputs(self): + materials = openmc.Materials() - fuel_mat = openmc.Material(material_id=1,name="UO2") + fuel_mat = openmc.Material(material_id=1, name="UO2") fuel_mat.set_density('sum') - fuel_mat.add_nuclide('U235',0.87370e-03) - fuel_mat.add_nuclide('U238',1.87440e-02) - fuel_mat.add_nuclide('O16',3.92350e-02) + fuel_mat.add_nuclide('U235', 0.87370e-03) + fuel_mat.add_nuclide('U238', 1.87440e-02) + fuel_mat.add_nuclide('O16', 3.92350e-02) materials.append(fuel_mat) - coolant = openmc.Material(material_id=2,name="borated H2O") + coolant = openmc.Material(material_id=2, name="borated H2O") coolant.set_density('sum') - coolant.add_nuclide('H1',0.06694) - coolant.add_nuclide('O16',0.03347) - coolant.add_nuclide('B10',6.6262e-6) - coolant.add_nuclide('B11',2.6839e-5) + coolant.add_nuclide('H1', 0.06694) + coolant.add_nuclide('O16', 0.03347) + coolant.add_nuclide('B10', 6.6262e-6) + coolant.add_nuclide('B11', 2.6839e-5) materials.append(coolant) - absorber = openmc.Material(material_id=3,name="pellet B4C") + absorber = openmc.Material(material_id=3, name="pellet B4C") absorber.set_density('sum') - absorber.add_nuclide('C0',0.01966) - absorber.add_nuclide('B11',4.7344e-6) - absorber.add_nuclide('B10',1.9177e-5) + absorber.add_nuclide('C0', 0.01966) + absorber.add_nuclide('B11', 4.7344e-6) + absorber.add_nuclide('B10', 1.9177e-5) materials.append(absorber) - zirc = openmc.Material(material_id=4,name="Zirc4") + zirc = openmc.Material(material_id=4, name="Zirc4") zirc.set_density('sum') - zirc.add_element('Zr',4.23e-2) + zirc.add_element('Zr', 4.23e-2) materials.append(zirc) materials.export_to_xml() - ### Geometry ### + # Geometry # + pin_rad = 0.7 # cm assembly_pitch = 1.235 # cm hexagonal_pitch = 23.6 # cm length = 10.0 # cm - #Fuel pin surfaces - cylfuelin = openmc.ZCylinder(surface_id=1,r=0.386) - cylfuelout = openmc.ZCylinder(surface_id=2,r=0.4582) - #Fuel cells + + # Fuel pin surfaces + + cylfuelin = openmc.ZCylinder(surface_id=1, r=0.386) + cylfuelout = openmc.ZCylinder(surface_id=2, r=0.4582) + + # Fuel cells + infcell = openmc.Cell(cell_id=1) infcell.region = -cylfuelin infcell.fill = fuel_mat @@ -55,14 +63,19 @@ class HexLatticeOXTestHarness(PyAPITestHarness): outfcell = openmc.Cell(cell_id=3) outfcell.region = +cylfuelout outfcell.fill = coolant - #Fuel universe - fuel_ch_univ = openmc.Universe(universe_id=1,name="Fuel channel", - cells=[infcell, clfcell,outfcell]) - - #Central tube surfaces - cyltubein = openmc.ZCylinder(surface_id=3,r=0.45) - cyltubeout = openmc.ZCylinder(surface_id=4,r=0.5177) - #Central tube cells + + # Fuel universe + + fuel_ch_univ = openmc.Universe(universe_id=1, name="Fuel channel", + cells=[infcell, clfcell, outfcell]) + + # Central tube surfaces + + cyltubein = openmc.ZCylinder(surface_id=3, r=0.45) + cyltubeout = openmc.ZCylinder(surface_id=4, r=0.5177) + + # Central tube cells + inctcell = openmc.Cell(cell_id=4) inctcell.region = -cyltubein inctcell.fill = coolant @@ -74,16 +87,22 @@ class HexLatticeOXTestHarness(PyAPITestHarness): outctcell = openmc.Cell(cell_id=6) outctcell.region = +cyltubeout outctcell.fill = coolant - #Central tubel universe - tube_ch_univ = openmc.Universe(universe_id=2,name="Central tube channel", - cells=[inctcell, clctcell,outctcell]) - #Absorber tube surfaces - cylabsin = openmc.ZCylinder(surface_id=5,r=0.35) - cylabsout = openmc.ZCylinder(surface_id=6,r=0.41) - cylabsclin = openmc.ZCylinder(surface_id=7,r=0.545) - cylabsclout = openmc.ZCylinder(surface_id=8,r=0.6323) - #Absorber tube cells + # Central tubel universe + + tube_ch_univ = openmc.Universe(universe_id=2, + name="Central tube channel", + cells=[inctcell, clctcell, outctcell]) + + # Absorber tube surfaces + + cylabsin = openmc.ZCylinder(surface_id=5, r=0.35) + cylabsout = openmc.ZCylinder(surface_id=6, r=0.41) + cylabsclin = openmc.ZCylinder(surface_id=7, r=0.545) + cylabsclout = openmc.ZCylinder(surface_id=8, r=0.6323) + + # Absorber tube cells + inabscell = openmc.Cell(cell_id=7) inabscell.region = -cylabsin inabscell.fill = absorber @@ -95,7 +114,7 @@ class HexLatticeOXTestHarness(PyAPITestHarness): interabscell = openmc.Cell(cell_id=9) interabscell.region = -cylabsclin & +cylabsout interabscell.fill = coolant - + clatcell = openmc.Cell(cell_id=10) clatcell.region = -cylabsclout & +cylabsclin clatcell.fill = zirc @@ -103,57 +122,69 @@ class HexLatticeOXTestHarness(PyAPITestHarness): outabscell = openmc.Cell(cell_id=11) outabscell.region = +cylabsclout outabscell.fill = coolant - #Absorber tube universe - abs_ch_univ = openmc.Universe(universe_id=3,name="Central tube channel",\ - cells=[inabscell,clabscell,interabscell,clatcell,outabscell]) - #Assembly surfaces - edge_length = (1./np.sqrt(3.0)) * hexagonal_pitch - fuel_bottom = openmc.ZPlane(surface_id=9,z0=0.0, - boundary_type = 'reflective') - fuel_top = openmc.ZPlane(surface_id=10,z0=length, - boundary_type = 'reflective') + + # Absorber tube universe + + abs_ch_univ = openmc.Universe(universe_id=3, + name="Central tube channel", + cells=[inabscell, clabscell, + interabscell, + clatcell, outabscell]) + # Assembly surfaces + + edge_length = (1./np.sqrt(3.0)) * hexagonal_pitch + fuel_bottom = openmc.ZPlane(surface_id=9, z0=0.0, + boundary_type='reflective') + fuel_top = openmc.ZPlane(surface_id=10, z0=length, + boundary_type='reflective') + # a hex surface for the core to go inside of + hexprism = openmc.model.get_hexagonal_prism(edge_length=edge_length, - origin=(0.0, 0.0), - boundary_type = 'reflective', + origin=(0.0, 0.0), + boundary_type='reflective', orientation='x') region = hexprism & +fuel_bottom & -fuel_top inf_mat = openmc.Cell(cell_id=12) inf_mat.fill = coolant - inf_mat_univ = openmc.Universe(universe_id=4,cells=[inf_mat,]) + inf_mat_univ = openmc.Universe(universe_id=4, cells=[inf_mat, ]) + # Fill lattice by channels + nring = 11 - universes = [] # array of rings - for ring in range(nring - 1,-1,-1): + universes = [] + for ring in range(nring - 1, -1, -1): arr = [] arr.append(fuel_ch_univ) for cell in range(ring * 6 - 1): arr.append(fuel_ch_univ) universes.append(arr) universes[-1] = [tube_ch_univ] - channels = [(7,16),(7,13),(7,10),(7,7),(7,4),(7,1),(5,0),\ - (4,33),(5,25),(4,27),(5,20),(4,21),(5,15),\ - (4,15),(5,10),(4,9),(5,5),(4,3)] - for i,j in channels: + channels = [(7, 16), (7, 13), (7, 10), (7, 7), (7, 4), (7, 1), (5, 0), + (4, 33), (5, 25), (4, 27), (5, 20), (4, 21), (5, 15), + (4, 15), (5, 10), (4, 9), (5, 5), (4, 3)] + for i, j in channels: universes[i][j] = abs_ch_univ lattice = openmc.HexLattice(name="regular fuel assembly") lattice.orientation = "ox" lattice.center = (0., 0., length/2.0) - lattice.pitch = (assembly_pitch,length/2.0) + lattice.pitch = (assembly_pitch, length/2.0) lattice.universes = 2*[universes] lattice.outer = inf_mat_univ - - assembly_cell = openmc.Cell(cell_id=13,name="container assembly cell") + + assembly_cell = openmc.Cell(cell_id=13, + name="container assembly cell") assembly_cell.region = region assembly_cell.fill = lattice - root_univ = openmc.Universe(universe_id=5,name="root universe", cells=[assembly_cell,]) + root_univ = openmc.Universe(universe_id=5, name="root universe", + cells=[assembly_cell, ]) geom = openmc.Geometry(root_univ) geom.export_to_xml() - ### Settings ### + # Settings # settings = openmc.Settings() settings.run_mode = 'eigenvalue' @@ -164,13 +195,15 @@ class HexLatticeOXTestHarness(PyAPITestHarness): source.space = openmc.stats.Box(ll, ur) source.strength = 1.0 settings.source = source - settings.output = {'summary' : False} + settings.output = {'summary': False} settings.batches = 10 settings.inactive = 5 settings.particles = 1000 settings.seed = 22 settings.export_to_xml() + def test_lattice_hex_ox_surf(): + harness = HexLatticeOXTestHarness('statepoint.10.h5') harness.main() From 98a95d048a39f0d666a8524e38a6928753e77c66 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Fri, 31 May 2019 10:30:40 +0300 Subject: [PATCH 24/31] test works with x --- openmc/lattice.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index fed12d192..43d7eb525 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -707,7 +707,7 @@ class RectLattice(Lattice): return (x, y, z) def get_universe_index(self, idx): - """Return index in the universes array corresponding + """Return index in the universes array corresponding to a lattice element index Parameters @@ -852,7 +852,7 @@ class RectLattice(Lattice): lat_id = int(get_text(elem, 'id')) name = get_text(elem, 'name') lat = cls(lat_id, name) - lat.lower_left = [float(i) + lat.lower_left = [float(i) for i in get_text(elem, 'lower_left').split()] lat.pitch = [float(i) for i in get_text(elem, 'pitch').split()] outer = get_text(elem, 'outer') @@ -943,8 +943,8 @@ class HexLattice(Lattice): string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', - "OX" if (self._orientation == 'x') - else "OY") + "x" if (self._orientation == 'x') + else "y") string += '{0: <16}{1}{2}\n'.format('\t# Rings', '=\t', self._num_rings) string += '{0: <16}{1}{2}\n'.format('\t# Axial', '=\t', self._num_axial) string += '{0: <16}{1}{2}\n'.format('\tCenter', '=\t', @@ -1032,7 +1032,7 @@ class HexLattice(Lattice): @orientation.setter def orientation(self, orientation): - cv.check_value('orientation', orientation.lower(), ('ox', 'oy')) + cv.check_value('orientation', orientation.lower(), ('x', 'y')) self._orientation = orientation.lower() @Lattice.pitch.setter @@ -1382,7 +1382,7 @@ class HexLattice(Lattice): # Create empty nested lists for one axial level univs = [[None for _ in range(max(6*(n_rings - 1 - r), 1))] - for r in range(n_rings)] + for r in range(n_rings)] if n_axial > 1: univs = [deepcopy(univs) for i in range(n_axial)] @@ -1434,7 +1434,7 @@ class HexLattice(Lattice): return self._repr_axial_slice_oy(universes) def _repr_axial_slice_ox(self, universes): - """Return string representation for the given 2D group of universes + """Return string representation for the given 2D group of universes in 'OX' orientation case. The 'universes' argument should be a list of lists of universes where @@ -1536,7 +1536,7 @@ class HexLattice(Lattice): return universe_ids def _repr_axial_slice_oy(self, universes): - """Return string representation for the given 2D group of universes in + """Return string representation for the given 2D group of universes in 'OY' orientation case.. The 'universes' argument should be a list of lists of universes where @@ -1644,6 +1644,7 @@ class HexLattice(Lattice): universe_ids = '\n'.join(rows) return universe_ids + @staticmethod def _show_indices_y(num_rings): """Return a diagram of the hexagonal lattice layout with indices. @@ -1747,8 +1748,9 @@ class HexLattice(Lattice): # Join the rows together and return the string. return '\n'.join(rows) + @staticmethod def _show_indices_x(num_rings): - """Return a diagram of the hexagonal lattice with OX orientation + """Return a diagram of the hexagonal lattice with OX orientation layout with indices. This method can be used to show the proper indices to be used when @@ -1862,7 +1864,7 @@ class HexLattice(Lattice): """ - if self._orientation == 'x': - return self._show_indices_x(num_rings) + if orientation == 'x': + return HexLattice._show_indices_x(num_rings) else: - return self._show_indices_y(num_rings) + return HexLattice._show_indices_y(num_rings) From a465dcdb1666bf15c0a6c9f859e06b800c0738cd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 4 Jun 2019 14:15:12 -0500 Subject: [PATCH 25/31] Consistently use 'x' / 'y' for hex lattice orientation --- include/openmc/lattice.h | 16 ++--- openmc/lattice.py | 36 ++++++------ src/lattice.cpp | 58 +++++++++---------- .../lattice_hex_OX/inputs_true.dat | 2 +- tests/regression_tests/lattice_hex_OX/test.py | 6 +- 5 files changed, 60 insertions(+), 58 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index c06920759..8f2b6c4a1 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -267,17 +267,19 @@ public: private: enum class Orientation { - oy, ox - }; - //! Fill universes_ vector for OY orientation - void fill_lattice_oy(const std::vector& univ_words); + y, //!< Flat side perpendicular to y-axis + x //!< Flat side perpendicular to x-axis + }; - //! Fill universes_ vector for OX orientation - void fill_lattice_ox(const std::vector& univ_words); + //! Fill universes_ vector for 'y' orientation + void fill_lattice_y(const std::vector& univ_words); + + //! Fill universes_ vector for 'x' orientation + void fill_lattice_x(const std::vector& univ_words); int n_rings_; //!< Number of radial tile positions int n_axial_; //!< Number of axial tile positions - Orientation orientation_; //!< Flat side up (oy) vs. sharp side up (ox) + Orientation orientation_; //!< Orientation of lattice Position center_; //!< Global center of lattice std::array pitch_; //!< Lattice tile width and height }; diff --git a/openmc/lattice.py b/openmc/lattice.py index 43d7eb525..4f7caa20a 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -141,7 +141,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): center = group['center'][()] pitch = group['pitch'][()] outer = group['outer'][()] - if ('orientation' in group): + if 'orientation' in group: orientation = group['orientation'][()].decode() else: orientation = "y" @@ -156,7 +156,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): # If the Universe specified outer the Lattice is not void if outer >= 0: lattice.outer = universes[outer] - if (orientation.lower() == "oy"): + if orientation == "y": # Build array of Universe pointers for the Lattice. Note that # we need to convert between the HDF5's square array of # (x, alpha, z) to the Python API's format of a ragged nested @@ -1207,8 +1207,8 @@ class HexLattice(Lattice): ---------- idx : Iterable of int Lattice element indices in the :math:`(x,\alpha,z)` coordinate - system in 'OY' orientation case, or indices in the - :math:`(\alpha,y,z)` coordinate system in 'OX' one + system in 'y' orientation case, or indices in the + :math:`(\alpha,y,z)` coordinate system in 'x' one Returns ------- @@ -1297,7 +1297,7 @@ class HexLattice(Lattice): lattice_subelement.set("n_rings", str(self._num_rings)) # If orientation is "OX" export it to XML if self._orientation == 'x': - lattice_subelement.set("orient", "OX") + lattice_subelement.set("orientation", "x") if self._num_axial is not None: lattice_subelement.set("n_axial", str(self._num_axial)) @@ -1429,13 +1429,13 @@ class HexLattice(Lattice): outer ring. """ if self._orientation == 'x': - return self._repr_axial_slice_ox(universes) + return self._repr_axial_slice_x(universes) else: - return self._repr_axial_slice_oy(universes) + return self._repr_axial_slice_y(universes) - def _repr_axial_slice_ox(self, universes): + def _repr_axial_slice_x(self, universes): """Return string representation for the given 2D group of universes - in 'OX' orientation case. + in 'x' orientation case. The 'universes' argument should be a list of lists of universes where each sub-list represents a single ring. The first list should be the @@ -1535,9 +1535,9 @@ class HexLattice(Lattice): universe_ids = '\n'.join(rows) return universe_ids - def _repr_axial_slice_oy(self, universes): + def _repr_axial_slice_y(self, universes): """Return string representation for the given 2D group of universes in - 'OY' orientation case.. + 'y' orientation case.. The 'universes' argument should be a list of lists of universes where each sub-list represents a single ring. The first list should be the @@ -1750,22 +1750,22 @@ class HexLattice(Lattice): @staticmethod def _show_indices_x(num_rings): - """Return a diagram of the hexagonal lattice with OX orientation + """Return a diagram of the hexagonal lattice with x orientation layout with indices. This method can be used to show the proper indices to be used when setting the :attr:`HexLattice.universes` property. For example,running this method with num_rings=3 will return the similar diagram:: - (0, 4) (0, 3) (0, 2) + (0, 4) (0, 3) (0, 2) - (0, 5) (1, 2) (1, 1) (0, 1) + (0, 5) (1, 2) (1, 1) (0, 1) - (0, 6) (1, 3) (2, 0) (1, 0) (0, 0) + (0, 6) (1, 3) (2, 0) (1, 0) (0, 0) - (0, 7) (1, 4) (1, 5) (0,11) + (0, 7) (1, 4) (1, 5) (0,11) - (0, 8) (0, 9) (0,10) + (0, 8) (0, 9) (0,10) Parameters ---------- @@ -1855,7 +1855,7 @@ class HexLattice(Lattice): num_rings : int Number of rings in the hexagonal lattice orientation : {"x", "y"} - str type of hex lattice orientation + Orientation of the hexagonal lattice Returns ------- diff --git a/src/lattice.cpp b/src/lattice.cpp index d2e29e106..1838df133 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -440,19 +440,19 @@ HexLattice::HexLattice(pugi::xml_node lat_node) is_3d_ = false; } - // Read the lattice orientation. Default to OY. - if (check_for_node(lat_node, "orient")) { - std::string orientation = get_node_value(lat_node, "orient"); - if (orientation == "OY") { - orientation_ = Orientation::oy; - } else if (orientation == "OX") { - orientation_ = Orientation::ox; + // Read the lattice orientation. Default to 'y'. + if (check_for_node(lat_node, "orientation")) { + std::string orientation = get_node_value(lat_node, "orientation"); + if (orientation == "y") { + orientation_ = Orientation::y; + } else if (orientation == "x") { + orientation_ = Orientation::x; } else { fatal_error("Unrecognized orientation '" + orientation + "' for lattice " + std::to_string(id_)); } } else { - orientation_ = Orientation::oy; + orientation_ = Orientation::y; } // Read the lattice center. @@ -497,25 +497,25 @@ HexLattice::HexLattice(pugi::xml_node lat_node) // Parse the universes. // Universes in hexagonal lattices are stored in a manner that represents - // a skewed coordinate system: (x, alpha) in case of OY orientation - // and (alpha,y) in OX one rather than (x, y). There is + // a skewed coordinate system: (x, alpha) in case of 'y' orientation + // and (alpha,y) in 'x' one rather than (x, y). There is // no obvious, direct relationship between the order of universes in the // input and the order that they will be stored in the skewed array so // the following code walks a set of index values across the skewed array // in a manner that matches the input order. Note that i_x = 0, i_a = 0 // or i_a = 0, i_y = 0 corresponds to the center of the hexagonal lattice. universes_.resize((2*n_rings_-1) * (2*n_rings_-1) * n_axial_, C_NONE); - if (orientation_ == Orientation::oy) { - fill_lattice_oy(univ_words); + if (orientation_ == Orientation::y) { + fill_lattice_y(univ_words); } else { - fill_lattice_ox(univ_words); + fill_lattice_x(univ_words); } } - + //============================================================================== - + void -HexLattice::fill_lattice_ox(const std::vector& univ_words) +HexLattice::fill_lattice_x(const std::vector& univ_words) { int input_index = 0; for (int m = 0; m < n_axial_; m++) { @@ -568,7 +568,7 @@ HexLattice::fill_lattice_ox(const std::vector& univ_words) //============================================================================== void -HexLattice::fill_lattice_oy(const std::vector& univ_words) +HexLattice::fill_lattice_y(const std::vector& univ_words) { int input_index = 0; for (int m = 0; m < n_axial_; m++) { @@ -692,13 +692,13 @@ const { // Short description of the direction vectors used here. The beta, gamma, and // delta vectors point towards the flat sides of each hexagonal tile. - // OY - orientation: + // Y - orientation: // basis0 = (1, 0) // basis1 = (-1/sqrt(3), 1) = +120 degrees from basis0 // beta = (sqrt(3)/2, 1/2) = +30 degrees from basis0 // gamma = (sqrt(3)/2, -1/2) = -60 degrees from beta // delta = (0, 1) = +60 degrees from beta - // OX - orientation: + // X - orientation: // basis0 = (1/sqrt(3), -1) // basis1 = (0, 1) = +120 degrees from basis0 // beta = (1, 0) = +30 degrees from basis0 @@ -708,7 +708,7 @@ const double beta_dir; double gamma_dir; double delta_dir; - if (orientation_ == Orientation::oy) { + if (orientation_ == Orientation::y) { beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0; delta_dir = u.y; @@ -737,7 +737,7 @@ const r_t = get_local_position(r, i_xyz_t); } double beta; - if (orientation_ == Orientation::oy) { + if (orientation_ == Orientation::y) { beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0; } else { beta = r_t.x; @@ -761,7 +761,7 @@ const r_t = get_local_position(r, i_xyz_t); } double gamma; - if (orientation_ == Orientation::oy) { + if (orientation_ == Orientation::y) { gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0; } else { gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0; @@ -788,7 +788,7 @@ const r_t = get_local_position(r, i_xyz_t); } double delta; - if (orientation_ == Orientation::oy) { + if (orientation_ == Orientation::y) { delta = r_t.y; } else { delta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0; @@ -848,7 +848,7 @@ HexLattice::get_indices(Position r, Direction u) const } int i1, i2; - if (orientation_ == Orientation::oy) { + if (orientation_ == Orientation::y) { // Convert coordinates into skewed bases. The (x, alpha) basis is used to // find the index of the global coordinates to within 4 cells. double alpha = r_o.y - r_o.x / std::sqrt(3.0); @@ -926,7 +926,7 @@ Position HexLattice::get_local_position(Position r, const std::array i_xyz) const { - if (orientation_ == Orientation::oy) { + if (orientation_ == Orientation::y) { // x_l = x_g - (center + pitch_x*cos(30)*index_x) r.x -= center_.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0]; @@ -941,7 +941,7 @@ const r.y -= center_.y + std::sqrt(3.0)/2.0 * (i_xyz[1] - n_rings_ + 1) * pitch_[0]; } - + if (is_3d_) { r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1]; } @@ -1003,10 +1003,10 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const write_string(lat_group, "type", "hexagonal", false); write_dataset(lat_group, "n_rings", n_rings_); write_dataset(lat_group, "n_axial", n_axial_); - if (orientation_ == Orientation::oy) { - write_string(lat_group, "orientation", "oy", false); + if (orientation_ == Orientation::y) { + write_string(lat_group, "orientation", "y", false); } else { - write_string(lat_group, "orientation", "ox", false); + write_string(lat_group, "orientation", "x", false); } if (is_3d_) { write_dataset(lat_group, "pitch", pitch_); diff --git a/tests/regression_tests/lattice_hex_OX/inputs_true.dat b/tests/regression_tests/lattice_hex_OX/inputs_true.dat index 9130cf284..9c0b102a9 100644 --- a/tests/regression_tests/lattice_hex_OX/inputs_true.dat +++ b/tests/regression_tests/lattice_hex_OX/inputs_true.dat @@ -13,7 +13,7 @@ - + 1.235 5.0 4
0.0 0.0 5.0
diff --git a/tests/regression_tests/lattice_hex_OX/test.py b/tests/regression_tests/lattice_hex_OX/test.py index bcef2357d..c958387dd 100644 --- a/tests/regression_tests/lattice_hex_OX/test.py +++ b/tests/regression_tests/lattice_hex_OX/test.py @@ -148,7 +148,7 @@ class HexLatticeOXTestHarness(PyAPITestHarness): inf_mat = openmc.Cell(cell_id=12) inf_mat.fill = coolant - inf_mat_univ = openmc.Universe(universe_id=4, cells=[inf_mat, ]) + inf_mat_univ = openmc.Universe(universe_id=4, cells=[inf_mat]) # Fill lattice by channels @@ -167,7 +167,7 @@ class HexLatticeOXTestHarness(PyAPITestHarness): for i, j in channels: universes[i][j] = abs_ch_univ lattice = openmc.HexLattice(name="regular fuel assembly") - lattice.orientation = "ox" + lattice.orientation = "x" lattice.center = (0., 0., length/2.0) lattice.pitch = (assembly_pitch, length/2.0) lattice.universes = 2*[universes] @@ -179,7 +179,7 @@ class HexLatticeOXTestHarness(PyAPITestHarness): assembly_cell.fill = lattice root_univ = openmc.Universe(universe_id=5, name="root universe", - cells=[assembly_cell, ]) + cells=[assembly_cell]) geom = openmc.Geometry(root_univ) geom.export_to_xml() From 3794c3eec44ba69a67b0028cb76954aa3f7c20a6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 4 Jun 2019 14:42:57 -0500 Subject: [PATCH 26/31] Update documentation for hex lattice --- docs/source/io_formats/geometry.rst | 7 +++++++ openmc/lattice.py | 6 +++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/docs/source/io_formats/geometry.rst b/docs/source/io_formats/geometry.rst index 94511842d..a89bfb1e3 100644 --- a/docs/source/io_formats/geometry.rst +++ b/docs/source/io_formats/geometry.rst @@ -318,6 +318,13 @@ the following attributes or sub-elements: *Default*: None + :orientation: + The orientation of the hexagonal lattice. The string "x" indicates that two + sides of the lattice are perpendicular to the x-axis, whereas the string "y" + indicates that two sides are perpendicular to the y-axis. + + *Default*: "y" + :center: The coordinates of the center of the lattice. If the lattice does not have axial sections then only the x- and y-coordinates are specified. diff --git a/openmc/lattice.py b/openmc/lattice.py index 4f7caa20a..b2bdc8a56 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -884,6 +884,10 @@ class HexLattice(Lattice): :attr:`HexLattice.universes` property, the array indices do not correspond to natural indices. + .. versionchanged:: 0.11 + The orientation of the lattice can now be changed with the + :attr:`orientation` attribute. + Parameters ---------- lattice_id : int, optional @@ -1295,7 +1299,7 @@ class HexLattice(Lattice): self._outer.create_xml_subelement(xml_element) lattice_subelement.set("n_rings", str(self._num_rings)) - # If orientation is "OX" export it to XML + # If orientation is "x" export it to XML if self._orientation == 'x': lattice_subelement.set("orientation", "x") From d1a2640ac350cc4acde81193c6ef9362e9237a92 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 4 Jun 2019 16:00:50 -0500 Subject: [PATCH 27/31] Rename lattice_hex_OX -> lattice_hex_x --- .../{lattice_hex_OX => lattice_hex_x}/__init__.py | 0 .../{lattice_hex_OX => lattice_hex_x}/inputs_true.dat | 0 .../{lattice_hex_OX => lattice_hex_x}/results_true.dat | 0 tests/regression_tests/{lattice_hex_OX => lattice_hex_x}/test.py | 0 4 files changed, 0 insertions(+), 0 deletions(-) rename tests/regression_tests/{lattice_hex_OX => lattice_hex_x}/__init__.py (100%) rename tests/regression_tests/{lattice_hex_OX => lattice_hex_x}/inputs_true.dat (100%) rename tests/regression_tests/{lattice_hex_OX => lattice_hex_x}/results_true.dat (100%) rename tests/regression_tests/{lattice_hex_OX => lattice_hex_x}/test.py (100%) diff --git a/tests/regression_tests/lattice_hex_OX/__init__.py b/tests/regression_tests/lattice_hex_x/__init__.py similarity index 100% rename from tests/regression_tests/lattice_hex_OX/__init__.py rename to tests/regression_tests/lattice_hex_x/__init__.py diff --git a/tests/regression_tests/lattice_hex_OX/inputs_true.dat b/tests/regression_tests/lattice_hex_x/inputs_true.dat similarity index 100% rename from tests/regression_tests/lattice_hex_OX/inputs_true.dat rename to tests/regression_tests/lattice_hex_x/inputs_true.dat diff --git a/tests/regression_tests/lattice_hex_OX/results_true.dat b/tests/regression_tests/lattice_hex_x/results_true.dat similarity index 100% rename from tests/regression_tests/lattice_hex_OX/results_true.dat rename to tests/regression_tests/lattice_hex_x/results_true.dat diff --git a/tests/regression_tests/lattice_hex_OX/test.py b/tests/regression_tests/lattice_hex_x/test.py similarity index 100% rename from tests/regression_tests/lattice_hex_OX/test.py rename to tests/regression_tests/lattice_hex_x/test.py From 853438083ae01fae6b18c5860e8e7157943c75ee Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 7 Jun 2019 15:25:09 -0500 Subject: [PATCH 28/31] Change order of x-oriented hex lattice universes to clockwise --- openmc/lattice.py | 96 +++++++++++++------- tests/regression_tests/lattice_hex_x/test.py | 6 +- 2 files changed, 64 insertions(+), 38 deletions(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index b2bdc8a56..30d9ad6de 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -1241,8 +1241,8 @@ class HexLattice(Lattice): else: i_within = 5*g - z - if (self._orientation == 'x') and (i_within > 0): - i_within = 6*(self._num_rings - i_ring - 1) - i_within_ + if self._orientation == 'x' and g > 0: + i_within = (i_within + 5*g) % (6*g) if self.num_axial is None: return (i_ring, i_within) @@ -1376,6 +1376,7 @@ class HexLattice(Lattice): lat = cls(lat_id, name) lat.center = [float(i) for i in get_text(elem, 'center').split()] lat.pitch = [float(i) for i in get_text(elem, 'pitch').split()] + lat.orientation = get_text(elem, 'orientation') outer = get_text(elem, 'outer') if outer is not None: lat.outer = get_universe(int(outer)) @@ -1390,7 +1391,7 @@ class HexLattice(Lattice): if n_axial > 1: univs = [deepcopy(univs) for i in range(n_axial)] - # Get flat array of universes numbers + # Get flat array of universes uarray = np.array([get_universe(int(i)) for i in get_text(elem, 'universes').split()]) @@ -1400,28 +1401,53 @@ class HexLattice(Lattice): # Get list for a single axial level axial_level = univs[z] if n_axial > 1 else univs - # Start iterating from top - x, alpha = 0, n_rings - 1 - while True: - # Set entry in list based on (x,alpha,z) coordinates - _, i_ring, i_within = lat.get_universe_index((x, alpha, z)) - axial_level[i_ring][i_within] = uarray[j] + if lat.orientation == 'y': + # Start iterating from top + x, alpha = 0, n_rings - 1 + while True: + # Set entry in list based on (x,alpha,z) coordinates + _, i_ring, i_within = lat.get_universe_index((x, alpha, z)) + axial_level[i_ring][i_within] = uarray[j] - # Move to the right - x += 2 - alpha -= 1 - if not lat.is_valid_index((x, alpha, z)): - # Move down in y direction - alpha += x - 1 - x = 1 - x + # Move to the right + x += 2 + alpha -= 1 if not lat.is_valid_index((x, alpha, z)): - # Move to the right - x += 2 - alpha -= 1 + # Move down in y direction + alpha += x - 1 + x = 1 - x if not lat.is_valid_index((x, alpha, z)): - # Reached the bottom + # Move to the right + x += 2 + alpha -= 1 + if not lat.is_valid_index((x, alpha, z)): + # Reached the bottom + break + j += 1 + else: + # Start iterating from top + alpha, y = 1 - n_rings, n_rings - 1 + while True: + # Set entry in list based on (alpha,y,z) coordinates + _, i_ring, i_within = lat.get_universe_index((alpha, y, z)) + axial_level[i_ring][i_within] = uarray[j] + + # Move to the right + alpha += 1 + if not lat.is_valid_index((alpha, y, z)): + # Move down to next row + alpha = 1 - n_rings + y -= 1 + + # Check if we've reached the bottom + if y == -n_rings: break - j += 1 + + while not lat.is_valid_index((alpha, y, z)): + # Move to the right + alpha += 1 + j += 1 + lat.universes = univs return lat @@ -1469,7 +1495,7 @@ class HexLattice(Lattice): theta = 0 y = middle - # Climb up the top-right. + # Climb down the bottom-right for i in range(r): # Add the universe. universe = universes[r_prime][theta] @@ -1479,7 +1505,7 @@ class HexLattice(Lattice): y += 1 theta += 1 - # Climb left the top-left. + # Climb left across the bottom for i in range(r): # Add the universe. universe = universes[r_prime][theta] @@ -1488,7 +1514,7 @@ class HexLattice(Lattice): # Translate the indices. theta += 1 - # Climb down the middle left. + # Climb up the bottom-left for i in range(r): # Add the universe. universe = universes[r_prime][theta] @@ -1498,7 +1524,7 @@ class HexLattice(Lattice): y -= 1 theta += 1 - # Climb down the down left. + # Climb up the top-left for i in range(r): # Add the universe. universe = universes[r_prime][theta] @@ -1508,7 +1534,7 @@ class HexLattice(Lattice): y -= 1 theta += 1 - # Climb right the down right. + # Climb right across the top for i in range(r): # Add the universe. universe = universes[r_prime][theta] @@ -1517,7 +1543,7 @@ class HexLattice(Lattice): # Translate the indices. theta += 1 - # Climb up the middle-right. + # Climb down the top-right for i in range(r): # Add the universe. universe = universes[r_prime][theta] @@ -1528,7 +1554,7 @@ class HexLattice(Lattice): theta += 1 # Flip the rows and join each row into a single string. - rows = [pad.join(x) for x in rows[::-1]] + rows = [pad.join(x) for x in rows] # Pad the beginning of the rows so they line up properly. for y in range(self._num_rings - 1): @@ -1806,41 +1832,41 @@ class HexLattice(Lattice): y = middle for i in range(r): - # Climb up the top-right. + # Climb down the bottom-right rows[y].append(str_form.format(r_prime, theta)) y += 1 theta += 1 for i in range(r): - # Climb left the top-left. + # Climb left across the bottom rows[y].insert(0, str_form.format(r_prime, theta)) theta += 1 for i in range(r): - # Climb down the middle left. + # Climb up the bottom-left rows[y].insert(0, str_form.format(r_prime, theta)) y -= 1 theta += 1 for i in range(r): - # Climb down the down left. + # Climb up the top-left rows[y].insert(0, str_form.format(r_prime, theta)) y -= 1 theta += 1 for i in range(r): - # Climb right the down right. + # Climb right across the top rows[y].append(str_form.format(r_prime, theta)) theta += 1 for i in range(r): - # Climb up the middle-right. + # Climb down the top-right rows[y].append(str_form.format(r_prime, theta)) y += 1 theta += 1 # Flip the rows and join each row into a single string. - rows = [pad.join(x) for x in rows[::-1]] + rows = [pad.join(x) for x in rows] # Pad the beginning of the rows so they line up properly. for y in range(num_rings - 1): diff --git a/tests/regression_tests/lattice_hex_x/test.py b/tests/regression_tests/lattice_hex_x/test.py index c958387dd..405e9d1e9 100644 --- a/tests/regression_tests/lattice_hex_x/test.py +++ b/tests/regression_tests/lattice_hex_x/test.py @@ -161,9 +161,9 @@ class HexLatticeOXTestHarness(PyAPITestHarness): arr.append(fuel_ch_univ) universes.append(arr) universes[-1] = [tube_ch_univ] - channels = [(7, 16), (7, 13), (7, 10), (7, 7), (7, 4), (7, 1), (5, 0), - (4, 33), (5, 25), (4, 27), (5, 20), (4, 21), (5, 15), - (4, 15), (5, 10), (4, 9), (5, 5), (4, 3)] + channels = [(7, 2), (7, 5), (7, 8), (7, 11), (7, 14), (7, 17), (5, 0), + (4, 3), (5, 5), (4, 9), (5, 10), (4, 15), (5, 15), + (4, 21), (5, 20), (4, 27), (5, 25), (4, 33)] for i, j in channels: universes[i][j] = abs_ch_univ lattice = openmc.HexLattice(name="regular fuel assembly") From d36adef2806947a39a1d25655f7af21683348a3e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 8 Jun 2019 06:57:57 -0500 Subject: [PATCH 29/31] Fix unit test using Geometry.from_xml on a hex lattice --- openmc/lattice.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index 30d9ad6de..0ac47b906 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -1376,7 +1376,7 @@ class HexLattice(Lattice): lat = cls(lat_id, name) lat.center = [float(i) for i in get_text(elem, 'center').split()] lat.pitch = [float(i) for i in get_text(elem, 'pitch').split()] - lat.orientation = get_text(elem, 'orientation') + lat.orientation = get_text(elem, 'orientation', 'y') outer = get_text(elem, 'outer') if outer is not None: lat.outer = get_universe(int(outer)) From e7f3e34088496b2c43d9e7bbb0f3da81d1644e3c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 10 Jun 2019 13:13:52 -0500 Subject: [PATCH 30/31] Fix from_hdf5 for hex lattice with x orientation. Assorted other fixes. --- docs/source/io_formats/geometry.rst | 4 +- include/openmc/lattice.h | 4 +- openmc/lattice.py | 49 ++++++++++--------- src/initialize.cpp | 4 +- src/relaxng/geometry.rnc | 2 +- src/relaxng/geometry.rng | 14 ++++-- .../lattice_hex_x/inputs_true.dat | 3 -- tests/regression_tests/lattice_hex_x/test.py | 3 -- tests/unit_tests/test_lattice.py | 9 ++++ 9 files changed, 50 insertions(+), 42 deletions(-) diff --git a/docs/source/io_formats/geometry.rst b/docs/source/io_formats/geometry.rst index a89bfb1e3..617ac6e87 100644 --- a/docs/source/io_formats/geometry.rst +++ b/docs/source/io_formats/geometry.rst @@ -320,8 +320,8 @@ the following attributes or sub-elements: :orientation: The orientation of the hexagonal lattice. The string "x" indicates that two - sides of the lattice are perpendicular to the x-axis, whereas the string "y" - indicates that two sides are perpendicular to the y-axis. + sides of the lattice are parallel to the x-axis, whereas the string "y" + indicates that two sides are parallel to the y-axis. *Default*: "y" diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index 8f2b6c4a1..cd73b9640 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -267,8 +267,8 @@ public: private: enum class Orientation { - y, //!< Flat side perpendicular to y-axis - x //!< Flat side perpendicular to x-axis + y, //!< Flat side of lattice parallel to y-axis + x //!< Flat side of lattice parallel to x-axis }; //! Fill universes_ vector for 'y' orientation diff --git a/openmc/lattice.py b/openmc/lattice.py index 0ac47b906..663d051e3 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -228,38 +228,38 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): # Add a list for this ring. uarray[-1].append([]) - # Climb up the top-left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, y, a]) - a -= 1 - y += 1 - - # Climb the left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, y, a]) - a -= 1 - - # Climb down the bottom-left. - for i in range(r): - uarray[-1][-1].append(universe_ids[z, y, a]) - y -= 1 - # Climb down the bottom-right. for i in range(r): uarray[-1][-1].append(universe_ids[z, y, a]) - a += 1 y -= 1 - # Climb up the right. + # Climb across the bottom. for i in range(r): uarray[-1][-1].append(universe_ids[z, y, a]) - a += 1 + a -= 1 - # Climb up the top-right. + # Climb up the bottom-left. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a -= 1 + y +=1 + + # Climb up the top-left. for i in range(r): uarray[-1][-1].append(universe_ids[z, y, a]) y += 1 + # Climb across the top. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a += 1 + + # Climb down the top-right. + for i in range(r): + uarray[-1][-1].append(universe_ids[z, y, a]) + a += 1 + y -= 1 + # Move down to the next ring. a -= 1 @@ -405,9 +405,11 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): ---------- idx : Iterable of int Lattice element indices. For a rectangular lattice, the indices are - given in the :math:`(x,y)` or :math:`(x,y,z)` coordinate system.For + given in the :math:`(x,y)` or :math:`(x,y,z)` coordinate system. For hexagonal lattices, they are given in the :math:`x,\alpha` or - :math:`x,\alpha,z` coordinate systems. + :math:`x,\alpha,z` coordinate systems for "y" orientations and + :math:`\alpha,y` or :math:`\alpha,y,z` coordinate systems for "x" + orientations. Returns ------- @@ -947,8 +949,7 @@ class HexLattice(Lattice): string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', - "x" if (self._orientation == 'x') - else "y") + self._orientation) string += '{0: <16}{1}{2}\n'.format('\t# Rings', '=\t', self._num_rings) string += '{0: <16}{1}{2}\n'.format('\t# Axial', '=\t', self._num_axial) string += '{0: <16}{1}{2}\n'.format('\tCenter', '=\t', diff --git a/src/initialize.cpp b/src/initialize.cpp index 061025745..f8a4f4a39 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -36,9 +36,7 @@ int openmc_init(int argc, char* argv[], const void* intracomm) { using namespace openmc; - - char cCurrentPath[256]; - + #ifdef OPENMC_MPI // Check if intracomm was passed MPI_Comm comm; diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index 30fdf2c8a..c34934071 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -36,7 +36,7 @@ element geometry { attribute dimension { list { xsd:positiveInteger+ } }) & (element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) & (element pitch { list { xsd:double+ } } | attribute pitch { list { xsd:double+ } }) & - (element orient { ( "OX" ) } | attribute orient { ( "OX" ) })? & + (element orientation { ( "x" | "y" ) } | attribute orientation { ( "x" | "y" ) })? & (element universes { list { xsd:int+ } } | attribute universes { list { xsd:int+ } }) & (element outer { xsd:int } | attribute outer { xsd:int })? }* diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng index e6c5f6ca1..52b9825be 100644 --- a/src/relaxng/geometry.rng +++ b/src/relaxng/geometry.rng @@ -296,11 +296,17 @@ - - OX + + + x + y + - - OX + + + x + y + diff --git a/tests/regression_tests/lattice_hex_x/inputs_true.dat b/tests/regression_tests/lattice_hex_x/inputs_true.dat index 9c0b102a9..b528a97b7 100644 --- a/tests/regression_tests/lattice_hex_x/inputs_true.dat +++ b/tests/regression_tests/lattice_hex_x/inputs_true.dat @@ -119,8 +119,5 @@ -13.62546635287517 -13.62546635287517 0.0 13.62546635287517 13.62546635287517 10.0 - - false - 22 diff --git a/tests/regression_tests/lattice_hex_x/test.py b/tests/regression_tests/lattice_hex_x/test.py index 405e9d1e9..538f790da 100644 --- a/tests/regression_tests/lattice_hex_x/test.py +++ b/tests/regression_tests/lattice_hex_x/test.py @@ -6,7 +6,6 @@ import numpy as np class HexLatticeOXTestHarness(PyAPITestHarness): def _build_inputs(self): - materials = openmc.Materials() fuel_mat = openmc.Material(material_id=1, name="UO2") @@ -195,7 +194,6 @@ class HexLatticeOXTestHarness(PyAPITestHarness): source.space = openmc.stats.Box(ll, ur) source.strength = 1.0 settings.source = source - settings.output = {'summary': False} settings.batches = 10 settings.inactive = 5 settings.particles = 1000 @@ -204,6 +202,5 @@ class HexLatticeOXTestHarness(PyAPITestHarness): def test_lattice_hex_ox_surf(): - harness = HexLatticeOXTestHarness('statepoint.10.h5') harness.main() diff --git a/tests/unit_tests/test_lattice.py b/tests/unit_tests/test_lattice.py index fe86e2a78..c4a351776 100644 --- a/tests/unit_tests/test_lattice.py +++ b/tests/unit_tests/test_lattice.py @@ -202,6 +202,13 @@ def test_get_universe(rlat2, rlat3, hlat2, hlat3): assert hlat2.get_universe((1, 0)) == u1 assert hlat2.get_universe((-2, 2)) == u1 + hlat2.orientation = 'x' + assert hlat2.get_universe((2, 0)) == u2 + assert hlat2.get_universe((1, 0)) == u2 + assert hlat2.get_universe((1, 1)) == u1 + assert hlat2.get_universe((-1, 1)) == u1 + hlat2.orientation = 'y' + u1, u2, u3, outer = hlat3.univs assert hlat3.get_universe((0, 0, 0)) == u2 assert hlat3.get_universe((0, 0, 1)) == u3 @@ -342,3 +349,5 @@ def test_show_indices(): for i in range(1, 11): lines = openmc.HexLattice.show_indices(i).split('\n') assert len(lines) == 4*i - 3 + lines_x = openmc.HexLattice.show_indices(i, 'x').split('\n') + assert len(lines) == 4*i - 3 From a6c099e29ae5e9c3b1795aeea3dd5030706dcb5d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 17 Jun 2019 06:24:45 -0500 Subject: [PATCH 31/31] A few more fixes pointed out by @smharper --- openmc/lattice.py | 16 ++++++++-------- src/relaxng/geometry.rnc | 2 +- src/relaxng/geometry.rng | 32 ++++++++++++++++---------------- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index 663d051e3..4776222da 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -1085,7 +1085,7 @@ class HexLattice(Lattice): # Check the center ring. if len(axial_slice[-1]) != 1: msg = 'HexLattice ID={0:d} has the wrong number of ' \ - 'elements in the innermost ring.Only 1 element is ' \ + 'elements in the innermost ring. Only 1 element is ' \ 'allowed in the innermost ring.'.format(self._id) raise ValueError(msg) @@ -1093,7 +1093,7 @@ class HexLattice(Lattice): for r in range(self._num_rings-1): if len(axial_slice[r]) != 6*(self._num_rings - 1 - r): msg = 'HexLattice ID={0:d} has the wrong number of ' \ - 'elements in ring number{1:d}(counting from the '\ + 'elements in ring number {1:d} (counting from the '\ 'outermost ring). This ring should have {2:d} ' \ 'elements.'.format(self._id, r, 6*(self._num_rings - 1 - r)) @@ -1788,15 +1788,15 @@ class HexLattice(Lattice): setting the :attr:`HexLattice.universes` property. For example,running this method with num_rings=3 will return the similar diagram:: - (0, 4) (0, 3) (0, 2) - - (0, 5) (1, 2) (1, 1) (0, 1) - - (0, 6) (1, 3) (2, 0) (1, 0) (0, 0) + (0, 8) (0, 9) (0,10) (0, 7) (1, 4) (1, 5) (0,11) - (0, 8) (0, 9) (0,10) + (0, 6) (1, 3) (2, 0) (1, 0) (0, 0) + + (0, 5) (1, 2) (1, 1) (0, 1) + + (0, 4) (0, 3) (0, 2) Parameters ---------- diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index c34934071..e3c88c445 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -36,7 +36,6 @@ element geometry { attribute dimension { list { xsd:positiveInteger+ } }) & (element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) & (element pitch { list { xsd:double+ } } | attribute pitch { list { xsd:double+ } }) & - (element orientation { ( "x" | "y" ) } | attribute orientation { ( "x" | "y" ) })? & (element universes { list { xsd:int+ } } | attribute universes { list { xsd:int+ } }) & (element outer { xsd:int } | attribute outer { xsd:int })? }* @@ -49,6 +48,7 @@ element geometry { (element n_axial { xsd:int } | attribute n_axial { xsd:int })? & (element center { list { xsd:double+ } } | attribute center { list { xsd:double+ } }) & (element pitch { list { xsd:double+ } } | attribute pitch { list { xsd:double+ } }) & + (element orientation { ( "x" | "y" ) } | attribute orientation { ( "x" | "y" ) })? & (element universes { list { xsd:int+ } } | attribute universes { list { xsd:int+ } }) & (element outer { xsd:int } | attribute outer { xsd:int })? }* diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng index 52b9825be..56bf38580 100644 --- a/src/relaxng/geometry.rng +++ b/src/relaxng/geometry.rng @@ -294,22 +294,6 @@ - - - - - x - y - - - - - x - y - - - - @@ -414,6 +398,22 @@ + + + + + x + y + + + + + x + y + + + +