Skip to content

Commit cd4687b

Browse files
GiuseppeOrlando878776GiuseppeOrlando878776
committed
Included an auxiliary function to check data
1 parent 562d9a5 commit cd4687b

File tree

3 files changed

+172
-0
lines changed

3 files changed

+172
-0
lines changed

‎1D_RIEMANN/test_configuration/All_topology_solver.hpp‎

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,9 @@ namespace fs = std::filesystem;
2121
#include "schemes/Rusanov_flux.hpp"
2222
#include "schemes/non_conservative_flux.hpp"
2323

24+
/*--- Define preprocessor to check whether to control data or not ---*/
25+
#define VERBOSE
26+
2427
// This is the class for the simulation of the all-topology model
2528
//
2629
template<std::size_t dim>
@@ -118,6 +121,8 @@ class All_Topology_Solver {
118121

119122
Number get_max_lambda() const; /*--- Compute the estimate of the maximum eigenvalue ---*/
120123

124+
void check_data(); /*--- Auxiliary routine to check if (small) spurious negative values are present ---*/
125+
121126
void update_auxiliary_fields(); /*--- Routine to update auxiliary fields for output and time step update ---*/
122127
};
123128

@@ -424,6 +429,57 @@ All_Topology_Solver<dim>::get_max_lambda() const {
424429
return global_res;
425430
}
426431

432+
// Auxiliary routine to check if spurious negative values arise
433+
//
434+
template<std::size_t dim>
435+
void All_Topology_Solver<dim>::check_data() {
436+
samurai::for_each_cell(mesh,
437+
[&](const auto& cell)
438+
{
439+
// Start with the volume fraction
440+
if(conserved_variables[cell][Indices::ALPHA1_INDEX] < static_cast<Number>(0.0)) {
441+
std::cerr << "Negative volume fraction for phase 1" << std::endl;
442+
save("_diverged", conserved_variables);
443+
exit(1);
444+
}
445+
else if(conserved_variables[cell][Indices::ALPHA1_INDEX] > static_cast<Number>(1.0)) {
446+
std::cerr << "Exceeding volume fraction for phase 1" << std::endl;
447+
save("_diverged", conserved_variables);
448+
exit(1);
449+
}
450+
else if(std::isnan(conserved_variables[cell][Indices::ALPHA1_INDEX])) {
451+
std::cerr << "NaN volume fraction for phase 1" << std::endl;
452+
save("_diverged", conserved_variables);
453+
exit(1);
454+
}
455+
456+
// Sanity check for m1
457+
if(conserved_variables[cell][Indices::ALPHA1_RHO1_INDEX] < static_cast<Number>(0.0)) {
458+
std::cerr << "Negative mass for phase 1" << std::endl;
459+
save("_diverged", conserved_variables);
460+
exit(1);
461+
}
462+
else if(std::isnan(conserved_variables[cell][Indices::ALPHA1_RHO1_INDEX])) {
463+
std::cerr << "NaN mass for phase 1" << std::endl;
464+
save("_diverged", conserved_variables);
465+
exit(1);
466+
}
467+
468+
// Sanity check for m2
469+
if(conserved_variables[cell][Indices::ALPHA2_RHO2_INDEX] < static_cast<Number>(0.0)) {
470+
std::cerr << "Negative mass for phase 2" << std::endl;
471+
save("_diverged", conserved_variables);
472+
exit(1);
473+
}
474+
else if(std::isnan(conserved_variables[cell][Indices::ALPHA2_RHO2_INDEX])){
475+
std::cerr << "NaN mass for phase 2" << std::endl;
476+
save("_diverged", conserved_variables);
477+
exit(1);
478+
}
479+
}
480+
);
481+
}
482+
427483
// Update auxiliary fields after solution of the system
428484
//
429485
template<std::size_t dim>
@@ -593,6 +649,9 @@ void All_Topology_Solver<dim>::run(const std::size_t nfiles) {
593649
- dt*Rusanov_flux(conserved_variables)
594650
- dt*NonConservative_flux(conserved_variables);
595651
samurai::swap(conserved_variables, conserved_variables_np1);
652+
#ifdef VERBOSE
653+
check_data();
654+
#endif
596655

597656
// Save the results
598657
if(t >= static_cast<Number>(nsave + 1)*dt_save || t == Tf) {

‎2D_RIEMANN/test_configuration/All_topology_solver.hpp‎

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,9 @@ namespace fs = std::filesystem;
2424
#include "schemes/Rusanov_flux.hpp"
2525
#include "schemes/non_conservative_flux.hpp"
2626

27+
/*--- Define preprocessor to check whether to control data or not ---*/
28+
#define VERBOSE
29+
2730
// This is the class for the simulation of the all-topology model
2831
//
2932
template<std::size_t dim>
@@ -118,6 +121,8 @@ class All_Topology_Solver {
118121

119122
Number get_max_lambda() const; /*--- Compute the estimate of the maximum eigenvalue ---*/
120123

124+
void check_data(); /*--- Auxiliary routine to check if (small) spurious negative values are present ---*/
125+
121126
void update_auxiliary_fields(); /*--- Routine to update auxiliary fields for output and time step update ---*/
122127
};
123128

@@ -433,6 +438,57 @@ All_Topology_Solver<dim>::get_max_lambda() const {
433438
return global_res;
434439
}
435440

441+
// Auxiliary routine to check if spurious negative values arise
442+
//
443+
template<std::size_t dim>
444+
void All_Topology_Solver<dim>::check_data() {
445+
samurai::for_each_cell(mesh,
446+
[&](const auto& cell)
447+
{
448+
// Start with the volume fraction
449+
if(conserved_variables[cell][Indices::ALPHA1_INDEX] < static_cast<Number>(0.0)) {
450+
std::cerr << "Negative volume fraction for phase 1" << std::endl;
451+
save("_diverged", conserved_variables);
452+
exit(1);
453+
}
454+
else if(conserved_variables[cell][Indices::ALPHA1_INDEX] > static_cast<Number>(1.0)) {
455+
std::cerr << "Exceeding volume fraction for phase 1" << std::endl;
456+
save("_diverged", conserved_variables);
457+
exit(1);
458+
}
459+
else if(std::isnan(conserved_variables[cell][Indices::ALPHA1_INDEX])) {
460+
std::cerr << "NaN volume fraction for phase 1" << std::endl;
461+
save("_diverged", conserved_variables);
462+
exit(1);
463+
}
464+
465+
// Sanity check for m1
466+
if(conserved_variables[cell][Indices::ALPHA1_RHO1_INDEX] < static_cast<Number>(0.0)) {
467+
std::cerr << "Negative mass for phase 1" << std::endl;
468+
save("_diverged", conserved_variables);
469+
exit(1);
470+
}
471+
else if(std::isnan(conserved_variables[cell][Indices::ALPHA1_RHO1_INDEX])) {
472+
std::cerr << "NaN mass for phase 1" << std::endl;
473+
save("_diverged", conserved_variables);
474+
exit(1);
475+
}
476+
477+
// Sanity check for m2
478+
if(conserved_variables[cell][Indices::ALPHA2_RHO2_INDEX] < static_cast<Number>(0.0)) {
479+
std::cerr << "Negative mass for phase 2" << std::endl;
480+
save("_diverged", conserved_variables);
481+
exit(1);
482+
}
483+
else if(std::isnan(conserved_variables[cell][Indices::ALPHA2_RHO2_INDEX])){
484+
std::cerr << "NaN mass for phase 2" << std::endl;
485+
save("_diverged", conserved_variables);
486+
exit(1);
487+
}
488+
}
489+
);
490+
}
491+
436492
// Update auxiliary fields after solution of the system
437493
//
438494
template<std::size_t dim>
@@ -599,6 +655,9 @@ void All_Topology_Solver<dim>::run(const std::size_t nfiles) {
599655
- dt*Rusanov_flux(conserved_variables)
600656
- dt*NonConservative_flux(conserved_variables);
601657
samurai::swap(conserved_variables, conserved_variables_np1);
658+
#ifdef VERBOSE
659+
check_data();
660+
#endif
602661

603662
// Save the results
604663
if(t >= static_cast<Number>(nsave + 1)*dt_save || t == Tf) {

‎include/schemes/Rusanov_flux.hpp‎

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,8 @@
88

99
#include "flux_base.hpp"
1010

11+
#define VERBOSE_FLUX
12+
1113
namespace samurai {
1214
/**
1315
* Implementation of a Rusanov flux
@@ -56,6 +58,22 @@ namespace samurai {
5658
const auto m2_L = qL(Indices::ALPHA2_RHO2_INDEX);
5759
const auto m2E2_L = qL(Indices::ALPHA2_RHO2_E2_INDEX);
5860

61+
// Verify if it is admissible
62+
#ifdef VERBOSE_FLUX
63+
if(m1_L < static_cast<Number>(0.0)) {
64+
throw std::runtime_error(std::string("Negative mass phase 1 left state: " + std::to_string(m1_L)));
65+
}
66+
if(m2_L < static_cast<Number>(0.0)) {
67+
throw std::runtime_error(std::string("Negative mass phase 2 left state: " + std::to_string(m2_L)));
68+
}
69+
if(alpha1_L < static_cast<Number>(0.0)) {
70+
throw std::runtime_error(std::string("Negative volume fraction phase 1 left state: " + std::to_string(alpha1_L)));
71+
}
72+
else if(alpha1_L > static_cast<Number>(1.0)) {
73+
throw std::runtime_error(std::string("Exceeding volume fraction phase 1 left state: " + std::to_string(alpha1_L)));
74+
}
75+
#endif
76+
5977
// Phase 1
6078
const auto inv_m1_L = static_cast<Number>(1.0)/m1_L; /*--- TODO: Add treatment for vanishing volume fraction ---*/
6179
const auto vel1_L_d = qL(Indices::ALPHA1_RHO1_U1_INDEX + curr_d)*inv_m1_L; /*--- TODO: Add treatment for vanishing volume fraction ---*/
@@ -68,6 +86,11 @@ namespace samurai {
6886
}
6987
const auto p1_L = this->EOS_phase1.pres_value_Rhoe(rho1_L, e1_L);
7088
const auto c1_L = this->EOS_phase1.c_value_RhoP(rho1_L, p1_L);
89+
#ifdef VERBOSE_FLUX
90+
if(std::isnan(c1_L)) {
91+
throw std::runtime_error(std::string("NaN speed of sound phase 1 left state"));
92+
}
93+
#endif
7194

7295
// Phase 2
7396
const auto inv_m2_L = static_cast<Number>(1.0)/m2_L; /*--- TODO: Add treatment for vanishing volume fraction ---*/
@@ -81,6 +104,11 @@ namespace samurai {
81104
}
82105
const auto p2_L = this->EOS_phase2.pres_value_Rhoe(rho2_L, e2_L);
83106
const auto c2_L = this->EOS_phase2.c_value_RhoP(rho2_L, p2_L);
107+
#ifdef VERBOSE_FLUX
108+
if(std::isnan(c2_L)) {
109+
throw std::runtime_error(std::string("NaN speed of sound phase 2 left state"));
110+
}
111+
#endif
84112

85113
/*--- Right state ---*/
86114
// Pre-fetch variables that will be used several times so as to exploit possible vectorization
@@ -91,6 +119,22 @@ namespace samurai {
91119
const auto m2_R = qR(Indices::ALPHA2_RHO2_INDEX);
92120
const auto m2E2_R = qR(Indices::ALPHA2_RHO2_E2_INDEX);
93121

122+
// Verify if it is admissible
123+
#ifdef VERBOSE_FLUX
124+
if(m1_R < static_cast<Number>(0.0)) {
125+
throw std::runtime_error(std::string("Negative mass phase 1 right state: " + std::to_string(m1_R)));
126+
}
127+
if(m2_R < static_cast<Number>(0.0)) {
128+
throw std::runtime_error(std::string("Negative mass phase 2 right state: " + std::to_string(m2_R)));
129+
}
130+
if(alpha1_R < static_cast<Number>(0.0)) {
131+
throw std::runtime_error(std::string("Negative volume fraction phase 1 right state: " + std::to_string(alpha1_R)));
132+
}
133+
else if(alpha1_R > static_cast<Number>(1.0)) {
134+
throw std::runtime_error(std::string("Exceeding volume fraction phase 1 right state: " + std::to_string(alpha1_R)));
135+
}
136+
#endif
137+
94138
// Phase 1
95139
const auto inv_m1_R = static_cast<Number>(1.0)/m1_R; /*--- TODO: Add treatment for vanishing volume fraction ---*/
96140
const auto vel1_R_d = qR(Indices::ALPHA1_RHO1_U1_INDEX + curr_d)*inv_m1_R; /*--- TODO: Add treatment for vanishing volume fraction ---*/
@@ -103,6 +147,11 @@ namespace samurai {
103147
}
104148
const auto p1_R = this->EOS_phase1.pres_value_Rhoe(rho1_R, e1_R);
105149
const auto c1_R = this->EOS_phase1.c_value_RhoP(rho1_R, p1_R);
150+
#ifdef VERBOSE_FLUX
151+
if(std::isnan(c1_R)) {
152+
throw std::runtime_error(std::string("NaN speed of sound phase 1 right state"));
153+
}
154+
#endif
106155

107156
// Phase 2
108157
const auto inv_m2_R = static_cast<Number>(1.0)/m2_R; /*--- TODO: Add treatment for vanishing volume fraction ---*/
@@ -116,6 +165,11 @@ namespace samurai {
116165
}
117166
const auto p2_R = this->EOS_phase2.pres_value_Rhoe(rho2_R, e2_R);
118167
const auto c2_R = this->EOS_phase2.c_value_RhoP(rho2_R, p2_R);
168+
#ifdef VERBOSE_FLUX
169+
if(std::isnan(c2_R)) {
170+
throw std::runtime_error(std::string("NaN speed of sound phase 2 right state"));
171+
}
172+
#endif
119173

120174
/*--- Compute the flux ---*/
121175
const auto lambda = std::max(std::max(std::abs(vel1_L_d) + c1_L, std::abs(vel1_R_d) + c1_R),

0 commit comments

Comments
 (0)