refactored matrixreal and sigp methods

CMakeLists.txt

@@ -85,6 +85,8 @@ message(STATUS "C++ flags: ${cxxflags}")
 #===============================================================================
 # Setup process
 #===============================================================================
+#Linking with thread library
+SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -pthread -lrt")
 
 list(APPEND lib_SOURCES
       "src/chain.cpp"

configure.xml

@@ -2,23 +2,16 @@
 <configure>
   <chain>./examples/chain_endfb71.xml</chain>
   <nuclides>./examples/nuclides.xml</nuclides>
-  <impxslibs>
-    <impxslib>./examples/XMAS172.xml</impxslib>
-  </impxslibs>
-  <inpmaterials>./examples/U235thn_1sec/materialsi.xml</inpmaterials>
-  <reaction>./examples/U235thn_1sec/reactions.xml</reaction>
-  <outmaterials>./examples/U235thn_1sec/materialso.xml</outmaterials>
-  <filters>
-    <filter type="exnuclide">U238 U239 U235 Pu238 Pu239 Pu240 Pu241 Pu242 Am241 Am242_m1</filter>
-  </filters>
-  <numbers>10</numbers>
-  <timestep>1</timestep>
-  <!-- IT IS  equivalent time description
-   timerecord year="0" month="0" day="0.0" hour="0" minute="0" second="1"/-->
+  <inpmaterials>/mnt/e/code/os/inpmaterialsOS.xml</inpmaterials>
+  <reaction>/mnt/e/code/os/reactionsOS.xml</reaction>
+  <outmaterials>/mnt/e/code/outmaterialsOS.xml</outmaterials>
+  <numbers>1</numbers>
+  <!--timestep>1</timestep-->
+  <timerecord year="0" month="0" day="592.0" hour="0" minute="0" second="0"></timerecord> 
   <method>chebyshev</method>  
   <epb>1.e-3</epb>
   <decaykey>false</decaykey>
-  <uncertanties>false</uncertanties>
-  <decay_print>false</decay_print>
+  <uncertanties>true</uncertanties>
+  <decay_print>true</decay_print>
   <cram_order>16</cram_order>
 </configure>

configureold.xml

@@ -0,0 +1,24 @@
+<?xml version="1.0"?>
+<configure>
+  <chain>./examples/chain_endfb71.xml</chain>
+  <nuclides>./examples/nuclides.xml</nuclides>
+  <impxslibs>
+    <impxslib>./examples/XMAS172.xml</impxslib>
+  </impxslibs>
+  <inpmaterials>./examples/U235thn_1sec/materialso.xml</inpmaterials>
+  <reaction>./examples/U235thn_1sec/reactions.xml</reaction>
+  <outmaterials>./examples/U235thn_1sec/materialso1.xml</outmaterials>
+  <filters>
+    <filter type="exnuclide">U238 U239 U235 Pu238 Pu239 Pu240 Pu241 Pu242 Am241 Am242_m1</filter>
+  </filters>
+  <numbers>10</numbers>
+  <timestep>1</timestep>
+  <!-- IT IS  equivalent time description
+   timerecord year="0" month="0" day="0.0" hour="0" minute="0" second="1"/-->
+  <method>chebyshev</method>  
+  <epb>1.e-3</epb>
+  <decaykey>false</decaykey>
+  <uncertanties>true</uncertanties>
+  <decay_print>true</decay_print>
+  <cram_order>16</cram_order>
+</configure>

include/openbps/executer.h

@@ -6,6 +6,7 @@
 #include "xtensor/xcomplex.hpp"
 #include "xtensor/xbuilder.hpp"
 #include "configure.h"
+#include "matrix.h"
 #include "chain.h"
 
 namespace openbps {
@@ -43,6 +44,7 @@ void run_solver();
 //!
 void init_solver();
 
+void executethr(int imat, Chain& chain, DecayMatrix& dm, DecayMatrix& ddm);
 //! Iterative method implementation v.1.
 //! Method based on algorythm described at paper Seleznev E.F., Belov A.A.,
 //! Belousov V.I., Chernova I.S. "BPSD CODE UPGRADE FOR SOLVING

include/openbps/matrix.h

@@ -117,9 +117,8 @@ class BaseMatrix
     //!
     //! The indexing operator [][]
     T* operator [] (int i) {
-        if (i >= nrows_)
-            return nullptr;
-        return data_[i];
+        if (i < nrows_)
+            return data_[i];
     }
     //! M1 | M2
     const BaseMatrix& operator |(const BaseMatrix& M2) {
@@ -162,7 +161,6 @@ class BaseMatrix
                     << m1.data_[i][j] << "   ";
             out << std::endl;
         }
-        return out;
     }
 
     //==========================================================================
@@ -342,6 +340,14 @@ class CramMatrix : public DecayMatrix {
     //!\return result a matrix with all transition for the material
     xt::xarray<double> matrixreal(Chain& chain,
                                   const Materials& mat);
+
+    //! Form a deviation decay nuclide matrix for unceratanties analysis
+    //!
+    //!\param[in] chain with decay information and fill the data storage
+    //!\param[in] material with nuclear concentration
+    //!\return result a matrix with all transition for the material
+    xt::xarray<double> matrixdev(Chain& chain,
+                                 const Materials& mat);
 }; //class ChebyshevMatrix
 
 //==============================================================================

outlog.csv

@@ -0,0 +1,22 @@
+dt;Act, sec-1;Q, Mev;dAct, sec-1;dQ, Mev;
+u235;
+dt;Act, sec-1;Q, Mev;dAct, sec-1;dQ, Mev;
+u235;
+0.1;4.80338e+15;1.85452e+16;1.40574e+14;5.32352e+14;
+0.2;5.12429e+15;1.97197e+16;1.50363e+14;5.6691e+14;
+0.3;5.43404e+15;2.08474e+16;1.59851e+14;6.00173e+14;
+0.4;5.73345e+15;2.19322e+16;1.69057e+14;6.32248e+14;
+0.5;6.02325e+15;2.29772e+16;1.78002e+14;6.63229e+14;
+0.6;6.30408e+15;2.39854e+16;1.86703e+14;6.93196e+14;
+0.7;6.57652e+15;2.49592e+16;1.95173e+14;7.22221e+14;
+0.8;6.84109e+15;2.59011e+16;2.03427e+14;7.50366e+14;
+0.9;7.09826e+15;2.6813e+16;2.11477e+14;7.77688e+14;
+1;7.34845e+15;2.76969e+16;2.19332e+14;8.04236e+14;
+dt;Act, sec-1;Q, Mev;dAct, sec-1;dQ, Mev;
+os16;
+dt;Act, sec-1;Q, Mev;dAct, sec-1;dQ, Mev;
+os16;
+5.11488e+07;1.31841e+14;2.19058e+14;1.71959e+12;4.45731e+12;
+dt;Act, sec-1;Q, Mev;dAct, sec-1;dQ, Mev;
+os16;
+5.11488e+07;1.31841e+14;2.19058e+14;1.71959e+12;4.45731e+12;

scripts/scriptprepro.py

@@ -268,8 +268,10 @@ def exec(nuclidnames, filepaths, execdir, nuclidlist):
     targetdir = sys.argv[-1]
     execdir = os.path.join(targetdir, "proceed")
     if (_MOD == "endfb"):
+    #Parsing ENDFB file structure lib
         nucls, paths = parse_endfblib(libdir)
     else:
+    #Parsing TENDL file structure lib
         nucls, paths = parse_tendlib(libdir)
     if (len(nuclidlist) == 0):
         nuclidlist = nucls

src/configure.cpp

@@ -71,6 +71,9 @@ int parse_command_line(int argc, char* argv[])
             if (arg == "-i" || arg == "--input") {
                 path_input = "";
             }
+            if (arg == "-v" || arg == "--verbose") {
+                verbose = true;
+            }
             if (arg == "-o" || arg == "--output") {
                 configure::outwrite = true;
             }

src/executer.cpp

@@ -1,5 +1,6 @@
 #include "openbps/executer.h"
 #include <vector>
+#include <thread>
 #include <cmath>
 #include <iostream>
 #include <fstream>
@@ -138,22 +139,13 @@ void init_solver() {
     }
 }
 
-//! Run a calculation process
-void run_solver() {
-    bool isMaterial {false};
-    if (configure::verbose)
-        std::cout << "We start execution\n";
-    // Read a chain from xml file
-    Chain chain = read_chain_xml(configure::chain_file);
-    // Form a calculation matrix
-    DecayMatrix dm(chain.name_idx.size());
-    dm.form_matrixreal(chain);
-    DecayMatrix ddm(chain.name_idx.size());
-    ddm.form_matrixdev(chain);
+void executethr(int imat, Chain& chain, DecayMatrix& dm, DecayMatrix& ddm) {
+    //std::vector<std::unique_ptr<Materials>>::iterator mat = materials.begin() + imat;
+//    auto mat = materials[imat];
     xt::xarray<double> dy;
-    for (auto& mat : materials) {
-        xt::xarray<double> y =
-                make_concentration(chain, mat->namenuclides, mat->conc);
+    bool isMaterial {false};
+    xt::xarray<double> y =
+                make_concentration(chain, materials[imat]->namenuclides, materials[imat]->conc);
         // Prepare dump to store every step of calculation results
         if (configure::outwrite) {
             configure::dumpoutput.clear();
@@ -171,15 +163,15 @@ void run_solver() {
         case Mode::iteration:
         {
             IterMatrix im(dm);
-            auto matrix = im.matrixreal(chain, *mat);
-            auto sigp = im.sigp(chain, *mat);
+            auto matrix = im.matrixreal(chain, *materials[imat]);
+            auto sigp = im.sigp(chain, *materials[imat]);
             // Perform calculation with uncertanties analysis
             if (configure::uncertantie_mod) {
-                dy = make_concentration(chain, mat->namenuclides,
-                                        mat->conc, true);
+                dy = make_concentration(chain,materials[imat]->namenuclides,
+                                        materials[imat]->conc, true);
                 IterMatrix imim(ddm);
-                auto dmatrix = imim.matrixdev(chain, *mat);
-                auto dsigp = imim.dsigp(chain, *mat);
+                auto dmatrix = imim.matrixdev(chain, *materials[imat]);
+                auto dsigp = imim.dsigp(chain, *materials[imat]);
                 iterative(matrix, sigp, y,
                           dmatrix, dsigp, dy);
             } else {
@@ -191,10 +183,30 @@ void run_solver() {
         case Mode::chebyshev:
         {
             CramMatrix cm(dm);
-            auto matrix = cm.matrixreal(chain, *mat);
+            auto matrix = cm.matrixreal(chain, *materials[imat]);
             if (configure::order == 8) {
+                std::vector<std::size_t> shape = {y.size()};
+                xt::xarray<double> dy2  = xt::zeros<double>(shape);
+
                 cram(matrix, y, alpha16, theta16,
                      configure::order, alpha160);
+                if (configure::uncertantie_mod) {
+                    configure::outwrite = false;
+
+                    CramMatrix devcm(ddm);
+                    dy = make_concentration(chain, materials[imat]->namenuclides,
+                                        materials[imat]->conc, true);
+                    auto dmatrix = devcm.matrixdev(chain, *materials[imat]);
+                    cram(matrix, dy, alpha16, theta16,
+                     configure::order, alpha160);
+
+                    std::copy(y.begin(), y.end(), dy2.begin());
+                    cram(dmatrix, dy2, alpha16, theta16,
+                     configure::order, alpha160);
+                    dy = dy + xt::abs(y - dy2);
+                    configure::outwrite = true;
+
+                }
             } else {
                 cram(matrix, y, alpha48, theta48,
                      configure::order, alpha480);
@@ -209,9 +221,9 @@ void run_solver() {
                              dy[j]<< std::endl;
             if (configure::rewrite) {
                 if (configure::uncertantie_mod) {
-                    mat->add_nuclide(item.first, udouble(y[j], dy[j]));
+                    materials[imat]->add_nuclide(item.first, udouble(y[j], dy[j]));
                 } else {
-                    mat->add_nuclide(item.first, y[j]);
+                    materials[imat]->add_nuclide(item.first, y[j]);
                 }
             }
             j++;
@@ -221,14 +233,37 @@ void run_solver() {
             // If materials fitler is present
             // then applying it
             if (materialfilter != nullptr) {
-                materialfilter->apply(mat->Name(), isMaterial);
+                materialfilter->apply(materials[imat]->Name(), isMaterial);
             } else {
                 isMaterial = true;
             }
             if (isMaterial)
-                apply_filters(chain, mat->Name());
+                apply_filters(chain, materials[imat]->Name());
          }
-    }//for materials
+}
+
+//! Run a calculation process
+void run_solver() {
+    bool isMaterial {false};
+    if (configure::verbose)
+        std::cout << "We start execution\n";
+    // Read a chain from xml file
+    Chain chain = read_chain_xml(configure::chain_file);
+    // Form a calculation matrix
+    DecayMatrix dm(chain.name_idx.size());
+    dm.form_matrixreal(chain);
+    DecayMatrix ddm(chain.name_idx.size());
+    ddm.form_matrixdev(chain);
+    std::vector<std::thread> threads;
+    for (int i = 0; i < materials.size(); i++) {
+        std::thread thr=std::thread(executethr, i, std::ref(chain), std::ref(dm), std::ref(ddm));
+        threads.emplace_back(std::move(thr));
+    }
+    for(auto& thr : threads) {
+        thr.join();
+    }
+
+    std::cout << "Done!" << std::endl;
     // Write down the getting nuclear concentration for every material
     // into *.xml file
     if (configure::rewrite)

src/matrix.cpp

@@ -379,6 +379,81 @@ xt::xarray<double> CramMatrix::matrixreal(Chain& chain,
     return result;
 }
 
+//! Form a deviation decay nuclide matrix for unceratanties analysis
+xt::xarray<double> CramMatrix::matrixdev(Chain& chain,
+                                           const Materials& mat) {
+    // Variables for calculation fissiop product yields by spectrum
+    std::pair<std::vector<double>, std::vector<double>> pair2;
+    std::vector<double> weight;
+    size_t k {0};
+    //
+    size_t NN {chain.name_idx.size()}; //!< Nuclide number
+    std::vector<std::size_t> shape = { NN, NN };
+    xt::xarray<double> result(shape, 0.0);
+    for (size_t i = 0; i < NN; i++) {
+        std::copy(&this->data_[i][0], &this->data_[i][NN],
+                  (result.begin() + i * NN));
+    }
+    int icompos {mat.numcomposition};
+    if (icompos > -1) {
+        // Get neutron flux - energy value descretization
+        pair2 = compositions[icompos]->get_fluxenergy();
+        for (auto it = chain.name_idx.begin();
+             it != chain.name_idx.end(); it++) {
+            size_t inucl = it->second; // from nuclide
+            size_t i = chain.get_nuclide_index(it->first);
+            // For xslib
+            for (auto& obj : compositions[icompos]->xslib) {
+            // If cross section presented for nuclides
+            if (obj.xsname == it->first) {
+                double rr {0.0};
+                // Sum reaction-rates over energy group
+                for (auto& r: obj.rxs)
+                    rr += r.Dev();
+                // Iterate over chain nuclides transition
+                for (auto& r : nuclides[inucl]->get_reactions()) {
+                    // If match
+                    if (r.first == obj.xstype) {
+                        if (obj.xstype != "fission") {
+                            // And non-fission then add
+                            size_t k = chain.get_nuclide_index(r.second);
+                            result(k, i) += rr * PWD * mat.normpower;
+                        } else {
+                            // Considering fission reaction by energy separately
+                            std::vector<double> energies =
+                                    nuclides[inucl]->get_nfy_energies();
+                            // Fission yields by product
+                            for (auto& item: nuclides[inucl]->
+                                     get_yield_product()) {
+                                k = chain.get_nuclide_index(item.first);
+                                double br {0.0};
+                                double norm {0.0};
+                                if (weight.empty())
+                                    weight = transition(pair2.first,
+                                                        pair2.second,
+                                                        energies);
+                                    for (int l = 0; l < weight.size(); l++) {
+                                         br += weight[l] * item.second[l];
+                                         norm += weight[l];
+                                    } // for weight
+
+                                    result(k, i) += br / norm * rr * PWD *
+                                                   mat.normpower;
+                                    norm = 1.0;
+                            } // for product
+                            weight.clear();
+                        } // fission yields
+                    } // if reaction = chain.reaction
+                } // run over reaction
+                result(i, i) -= rr * PWD * mat.normpower;
+            } // if nuclide is in crossection data and chain
+        } // for xslib in composition
+    } // if composition is presented
+    } // for all nuclides
+
+    return result;
+}
+
 //==============================================================================
 // Non class methods
 //==============================================================================