Skip to content

Commit 1fa9cdc

Browse files
van4elottivan4elotti
committed
refactored matrixreal and sigp methods
1 parent e87e822 commit 1fa9cdc

2 files changed

Lines changed: 56 additions & 116 deletions

File tree

‎include/openbps/matrix.h‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -262,6 +262,14 @@ class DecayMatrix : public BaseMatrix<double> {
262262
//! Form a real decay nuclide matrix
263263
//!
264264
//!\param[in] chain with decay information and fill the data storage
265+
//!\param[in] material with nuclear concentration
266+
//!\return result a matrix with all transition for the material
267+
xt::xarray<double> get_matrix_real(Chain& chain,
268+
const Materials& mat, std::string matrix_kind);
269+
//! Methods
270+
//! Form a real decay nuclide matrix
271+
//!
272+
//!\param[in] chain with decay information and fill the data storage
265273
void form_matrixreal(Chain& chain);
266274
//! Form a deviation decay nuclide matrix for unceratanties analysis
267275
//!
@@ -315,6 +323,16 @@ class IterMatrix : public DecayMatrix {
315323
//!\return result a vector with all transition for every nuclide:Dev
316324
xt::xarray<double> dsigp(Chain& chain,
317325
const Materials& mat);
326+
//! Utility for functions sigp(...) and dsigp(...)
327+
//!
328+
//!\param[in] chain with decay information and fill the data storage
329+
//!\param[in] material with nuclear concentration
330+
//!\param[in] boolean var (true for dsigp, false for sigps)
331+
//!\return result a vector with all transition for every
332+
//!nuclide:Dev in case of d == true nuclide:Real in case d == false
333+
xt::xarray<double> formSigp(Chain& chain,
334+
const Materials& mat,
335+
bool d);
318336

319337
}; //class IterMatrix
320338

‎src/matrix.cpp‎

Lines changed: 38 additions & 116 deletions
Original file line numberDiff line numberDiff line change
@@ -75,14 +75,8 @@ void DecayMatrix::form_matrixdev(Chain& chain) {
7575
}
7676
}
7777

78-
//==============================================================================
79-
// IterativeMatrix implementation
80-
//==============================================================================
81-
82-
//! Methods
83-
//! Form a real decay nuclide matrix
84-
xt::xarray<double> IterMatrix::matrixreal(Chain& chain,
85-
const Materials& mat) {
78+
xt::xarray<double> DecayMatrix::get_matrix_real(Chain& chain,
79+
const Materials& mat, std::string matrix_kind) {
8680
// Variables for calculation fissiop product yields by spectrum
8781
std::pair<std::vector<double>, std::vector<double>> pair2;
8882
std::vector<double> weight;
@@ -94,7 +88,8 @@ xt::xarray<double> IterMatrix::matrixreal(Chain& chain,
9488
for (size_t i = 0; i < NN; i++) {
9589
std::copy(&this->data_[i][0], &this->data_[i][NN],
9690
(result.begin() + i * NN));
97-
result(i, i) = 0.0;
91+
if (matrix_kind == "iter")
92+
result(i, i) = 0.0;
9893
}
9994
int icompos {mat.numcomposition};
10095
if (icompos > -1) {
@@ -147,14 +142,27 @@ xt::xarray<double> IterMatrix::matrixreal(Chain& chain,
147142
} // fission yields
148143
} // if reaction = chain.reaction
149144
} // run over reaction
145+
if (matrix_kind == "cram")
146+
result(i, i) -= rr * PWD * mat.normpower;
150147
} // if nuclide is in crossection data and chain
151148
} // for xslib in composition
152149
} // if composition is presented
153150
} // for all nuclides
154-
155151
return result;
156152
}
157153

154+
//==============================================================================
155+
// IterativeMatrix implementation
156+
//==============================================================================
157+
158+
//! Methods
159+
//! Form a real decay nuclide matrix
160+
xt::xarray<double> IterMatrix::matrixreal(Chain& chain,
161+
const Materials& mat) {
162+
// Variables for calculation fissiop product yields by spectrum
163+
return DecayMatrix::get_matrix_real(chain, mat, "iter");
164+
}
165+
158166
//! Form a deviation decay nuclide matrix for unceratanties analysis
159167
xt::xarray<double> IterMatrix::matrixdev(Chain& chain,
160168
const Materials& mat) {
@@ -169,7 +177,7 @@ xt::xarray<double> IterMatrix::matrixdev(Chain& chain,
169177
for (size_t i = 0; i < NN; i++) {
170178
std::copy(&this->data_[i][0], &this->data_[i][NN],
171179
(result.begin() + i * NN));
172-
result(i, i) = 0.0;
180+
result(i, i) = 0.0; // ------------------
173181
}
174182
int icompos {mat.numcomposition};
175183
if (icompos > -1) {
@@ -186,7 +194,7 @@ xt::xarray<double> IterMatrix::matrixdev(Chain& chain,
186194
double rr {0.0};
187195
// Sum reaction-rates over energy group
188196
for (auto& r: obj.rxs)
189-
rr += r.Dev();
197+
rr += r.Dev();//------------------
190198
// Iterate over chain nuclides transition
191199
for (auto& r : nuclides[inucl]->get_reactions()) {
192200
// If match
@@ -230,22 +238,23 @@ xt::xarray<double> IterMatrix::matrixdev(Chain& chain,
230238
return result;
231239
}
232240

233-
//! Form a dev decay nuclide vector for all transition from every nuclide
234-
xt::xarray<double> IterMatrix::sigp(Chain& chain,
235-
const Materials& mat) {
241+
242+
xt::xarray<double> IterMatrix::formSigp(Chain& chain,
243+
const Materials& mat,
244+
bool d) {
236245
std::vector<size_t> shape {chain.name_idx.size()};
237246
xt::xarray<double> result(shape, 0.0);
238247
udouble decay_ {0.0, 0.0};
239248
//Run over all nuclides
240-
if (configure::verbose)
249+
if (configure::verbose && !d)
241250
std::cout << "SIGP: "<<std::endl;
242251
for (auto it = chain.name_idx.begin(); it != chain.name_idx.end(); it++) {
243252
size_t inucl = it->second; // from nuclide
244253

245254
size_t i = chain.get_nuclide_index(it->first);
246255
if (nuclides[inucl]->half_life.Real() > 0) {
247256
decay_ = (log(2.0) / nuclides[inucl]->half_life);
248-
result(i) = decay_.Real();
257+
result(i) = d ? decay_.Dev() : decay_.Real();
249258
}
250259
int icompos {mat.numcomposition};
251260
if (icompos > -1) {
@@ -256,51 +265,31 @@ xt::xarray<double> IterMatrix::sigp(Chain& chain,
256265
double rr {0.0};
257266
// Sum reaction-rates over energy group
258267
for (auto& r: obj.rxs)
259-
rr += r.Real();
268+
rr += d ? r.Dev() : r.Real() ;
260269
result(i) += rr * PWD * mat.normpower;
261270
}
262271
}
263272
}
264-
if (configure::verbose)
273+
if (configure::verbose && !d)
265274
std::cout << "sigp: " << it->first << " " << result(i)<<std::endl;
266275
}
267276

268277
return result;
269278
}
270279

271280
//! Form a dev decay nuclide vector for all transition from every nuclide
272-
xt::xarray<double> IterMatrix::dsigp(Chain& chain,
273-
const Materials& mat) {
274-
std::vector<size_t> shape {chain.name_idx.size()};
275-
xt::xarray<double> result(shape, 0.0);
276-
udouble decay_ {0.0, 0.0};
277-
//Run over all nuclides
278-
for (auto it = chain.name_idx.begin(); it != chain.name_idx.end(); it++) {
279-
size_t inucl = it->second; // from nuclide
281+
xt::xarray<double> IterMatrix::sigp(Chain& chain,
282+
const Materials& mat) {
283+
return formSigp(chain, mat, false);
284+
}
280285

281-
size_t i = chain.get_nuclide_index(it->first);
282-
if (nuclides[inucl]->half_life.Real() > 0) {
283-
decay_ = (log(2.0) / nuclides[inucl]->half_life);
284-
result(i) = decay_.Dev();
285-
}
286-
int icompos {mat.numcomposition};
287-
if (icompos > -1) {
288-
// For xslib
289-
for (auto& obj : compositions[icompos]->xslib) {
290-
// If cross section presented for nuclides
291-
if (obj.xsname == it->first) {
292-
double rr {0.0};
293-
// Sum reaction-rates over energy group
294-
for (auto& r: obj.rxs)
295-
rr += r.Dev();
296-
result(i) += rr * PWD * mat.normpower;
297-
}
298-
}
299-
}
300-
}
301-
return result;
286+
//! Form a dev decay nuclide vector for all transition from every nuclide
287+
xt::xarray<double> IterMatrix::dsigp(Chain& chain,
288+
const Materials& mat) {
289+
return formSigp(chain, mat, true);
302290
}
303291

292+
304293
//==============================================================================
305294
// ChebyshevMatrix implementation
306295
//==============================================================================
@@ -309,74 +298,7 @@ xt::xarray<double> IterMatrix::dsigp(Chain& chain,
309298
xt::xarray<double> CramMatrix::matrixreal(Chain& chain,
310299
const Materials& mat) {
311300
// Variables for calculation fissiop product yields by spectrum
312-
std::pair<std::vector<double>, std::vector<double>> pair2;
313-
std::vector<double> weight;
314-
size_t k {0};
315-
//
316-
size_t NN {chain.name_idx.size()}; //!< Nuclide number
317-
std::vector<std::size_t> shape = { NN, NN };
318-
xt::xarray<double> result(shape, 0.0);
319-
for (size_t i = 0; i < NN; i++)
320-
std::copy(&this->data_[i][0], &this->data_[i][NN],
321-
(result.begin() + i * NN));
322-
int icompos {mat.numcomposition};
323-
if (icompos > -1) {
324-
// Get neutron flux - energy value descretization
325-
pair2 = compositions[icompos]->get_fluxenergy();
326-
for (auto it = chain.name_idx.begin();
327-
it != chain.name_idx.end(); it++) {
328-
size_t inucl = it->second; // from nuclide
329-
size_t i = chain.get_nuclide_index(it->first);
330-
// For xslib
331-
for (auto& obj : compositions[icompos]->xslib) {
332-
// If cross section presented for nuclides
333-
if (obj.xsname == it->first) {
334-
double rr {0.0};
335-
// Sum reaction-rates over energy group
336-
for (auto& r: obj.rxs)
337-
rr += r.Real();
338-
// Iterate over chain nuclides transition
339-
for (auto& r : nuclides[inucl]->get_reactions()) {
340-
// If match
341-
if (r.first == obj.xstype) {
342-
if (obj.xstype != "fission") {
343-
// And non-fission then add
344-
size_t k = chain.get_nuclide_index(r.second);
345-
result(k, i) += rr * PWD * mat.normpower;
346-
} else {
347-
// Considering fission reaction by energy separately
348-
std::vector<double> energies =
349-
nuclides[inucl]->get_nfy_energies();
350-
// Fission yields by product
351-
for (auto& item: nuclides[inucl]->
352-
get_yield_product()) {
353-
k = chain.get_nuclide_index(item.first);
354-
double br {0.0};
355-
double norm {0.0};
356-
if (weight.empty())
357-
weight = transition(pair2.first,
358-
pair2.second,
359-
energies);
360-
for (int l = 0; l < weight.size(); l++) {
361-
br += weight[l] * item.second[l];
362-
norm += weight[l];
363-
} // for weight
364-
365-
result(k, i) += br / norm * rr * PWD *
366-
mat.normpower;
367-
norm = 1.0;
368-
} // for product
369-
weight.clear();
370-
} // fission yields
371-
} // if reaction = chain.reaction
372-
} // run over reaction
373-
result(i, i) -= rr * PWD * mat.normpower;
374-
} // if nuclide is in crossection data and chain
375-
} // for xslib in composition
376-
} // if composition is presented
377-
} // for all nuclides
378-
379-
return result;
301+
return DecayMatrix::get_matrix_real(chain, mat, "cram");
380302
}
381303

382304
//==============================================================================

0 commit comments

Comments
 (0)