diff --git a/models/hymod/include/Hymod.h b/models/hymod/include/Hymod.h index fbdc75df24..916bdb065b 100644 --- a/models/hymod/include/Hymod.h +++ b/models/hymod/include/Hymod.h @@ -120,11 +120,11 @@ class hymod_kernel for ( unsigned long i = 0; i < nash_cascade.size(); ++i ) { //construct a single outlet nonlinear reservoir - nash_cascade[i] = Nonlinear_Reservoir(0, params.max_storage_meters, state.Sr[i], params.Kq, 1, 0); + nash_cascade[i] = Nonlinear_Reservoir(0, params.max_storage_meters, state.Sr[i], params.Kq, 1, 0, 100); } // initalize groundwater reservoir - Nonlinear_Reservoir groundwater(0, params.max_storage_meters, state.groundwater_storage_meters, params.Ks, 1, 0); + Nonlinear_Reservoir groundwater(0, params.max_storage_meters, state.groundwater_storage_meters, params.Ks, 1, 0, 100); // add flux to the current state state.storage_meters += input_flux_meters; diff --git a/models/kernels/Nonlinear_Reservoir.hpp b/models/kernels/Nonlinear_Reservoir.hpp index b17040ebb2..a7c9e9e6f9 100644 --- a/models/kernels/Nonlinear_Reservoir.hpp +++ b/models/kernels/Nonlinear_Reservoir.hpp @@ -26,13 +26,13 @@ class Reservoir_Outlet public: //Default Constructor - Reservoir_Outlet(): activation_threshold_meters(0.0), a(0.0), b(0.0) + Reservoir_Outlet(): a(0.0), b(0.0), activation_threshold_meters(0.0), max_velocity_meters_per_second(0.0) { } //Parameterized Constructor - Reservoir_Outlet(double a, double b, double activation_threshold_meters): activation_threshold_meters(activation_threshold_meters), a(a), b(b) + Reservoir_Outlet(double a, double b, double activation_threshold_meters, double max_velocity_meters_per_second): a(a), b(b), activation_threshold_meters(activation_threshold_meters), max_velocity_meters_per_second(max_velocity_meters_per_second) { } @@ -40,13 +40,27 @@ class Reservoir_Outlet //Function to return the velocity in meters per second of the discharge through the outlet double velocity_meters_per_second(reservoir_parameters ¶meters_struct, reservoir_state &storage_struct) { + double velocity_meters_per_second_local; + //Return velocity of 0.0 if the storage passed in is less than the activation threshold if (storage_struct.current_storage_height_meters <= activation_threshold_meters) return 0.0; + //Calculate the velocity in meters per second of the discharge through the outlet + velocity_meters_per_second_local = a * std::pow((storage_struct.current_storage_height_meters - activation_threshold_meters) + / (parameters_struct.maximum_storage_meters - activation_threshold_meters), b); + + //If calculated oulet velocity is greater than max velocity, then set to max velocity and return a warning. + if (velocity_meters_per_second_local > max_velocity_meters_per_second) + { + velocity_meters_per_second_local = max_velocity_meters_per_second; + + //TODO: Return appropriate warning + cout << "WARNING: Nonlinear reservoir calculated an outlet velocity over max velocity, and therefore set the outlet velocity to max velocity." << endl; + } + //Return the velocity in meters per second of the discharge through the outlet - return a * std::pow((storage_struct.current_storage_height_meters - activation_threshold_meters) - / (parameters_struct.maximum_storage_meters - activation_threshold_meters), b); + return velocity_meters_per_second_local; }; //Accessor to return activation_threshold_meters @@ -59,6 +73,7 @@ class Reservoir_Outlet double a; double b; double activation_threshold_meters; + double max_velocity_meters_per_second; }; @@ -82,13 +97,13 @@ class Nonlinear_Reservoir //Constructor for a reservoir with only one outlet. Nonlinear_Reservoir(double minimum_storage_meters, double maximum_storage_meters, double current_storage_height_meters, double a, - double b, double activation_threshold_meters) : Nonlinear_Reservoir(minimum_storage_meters, maximum_storage_meters, current_storage_height_meters) + double b, double activation_threshold_meters, double max_velocity_meters_per_second) : Nonlinear_Reservoir(minimum_storage_meters, maximum_storage_meters, current_storage_height_meters) { //Ensure that the activation threshold is less than the maximum storage //if (activation_threshold_meters > maximum_storage_meters) //TODO: Return appropriate error - this->outlets.push_back(Reservoir_Outlet(a, b, activation_threshold_meters)); + this->outlets.push_back(Reservoir_Outlet(a, b, activation_threshold_meters, max_velocity_meters_per_second)); } //Constructor for a reservoir with multiple outlets. @@ -130,7 +145,7 @@ class Nonlinear_Reservoir if (state.current_storage_height_meters < parameters.minimum_storage_meters) { //TODO: Return appropriate warning - cout << "WARNING: Nonlinear reservoir calculated a storage below the minimum storage." << endl; + //cout << "WARNING: Nonlinear reservoir calculated a storage below the minimum storage." << endl; //Return to storage before falling below minimum storage. state.current_storage_height_meters += outlet_velocity_meters_per_second * delta_time_seconds; diff --git a/test/models/hymod/include/Nonlinear_Reservoir_Test.cpp b/test/models/hymod/include/Nonlinear_Reservoir_Test.cpp index 38229725cd..ce70fbc6cc 100644 --- a/test/models/hymod/include/Nonlinear_Reservoir_Test.cpp +++ b/test/models/hymod/include/Nonlinear_Reservoir_Test.cpp @@ -16,6 +16,8 @@ class NonlinearReservoirKernelTest : public ::testing::Test { void setupOneOutletNonlinearReservoir(); + void setupOneOutletHighStorageNonlinearReservoir(); + void setupMultipleOutletNonlinearReservoir(); void setupMultipleOutletOutOfOrderNonlinearReservoir(); @@ -24,6 +26,8 @@ class NonlinearReservoirKernelTest : public ::testing::Test { std::shared_ptr OneOutletReservoir; //smart pointer to a Nonlinear_Reservoir with one outlet + std::shared_ptr OneOutletHighStorageReservoir; //smart pointer to a Nonlinear_Reservoir with one outlet and high storage + std::shared_ptr MultipleOutletReservoir; //smart pointer to a Nonlinear_Reservoir with multiple outlets std::shared_ptr MultipleOutletOutOfOrderNonlinearReservoir; //smart pointer to a Nonlinear_Reservoir with multiple outlets out of order @@ -46,6 +50,8 @@ void NonlinearReservoirKernelTest::SetUp() { setupOneOutletNonlinearReservoir(); + setupOneOutletHighStorageNonlinearReservoir(); + setupMultipleOutletNonlinearReservoir(); setupMultipleOutletOutOfOrderNonlinearReservoir(); @@ -64,17 +70,23 @@ void NonlinearReservoirKernelTest::setupNoOutletNonlinearReservoir() //Construct a reservoir with one outlet void NonlinearReservoirKernelTest::setupOneOutletNonlinearReservoir() { - OneOutletReservoir = std::make_shared(0.0, 8.0, 3.5, 0.5, 0.7, 4.0); + OneOutletReservoir = std::make_shared(0.0, 8.0, 3.5, 0.5, 0.7, 4.0, 100.0); +} + +//Construct a reservoir with one outlet and high storage +void NonlinearReservoirKernelTest::setupOneOutletHighStorageNonlinearReservoir() +{ + OneOutletHighStorageReservoir = std::make_shared(0.0, 8000.0, 3.5, 1.1, 1.2, 4.0, 0.005); } //Construct a reservoir with multiple outlets void NonlinearReservoirKernelTest::setupMultipleOutletNonlinearReservoir() { - ReservoirOutlet1 = std::make_shared(0.2, 0.4, 4.0); + ReservoirOutlet1 = std::make_shared(0.2, 0.4, 4.0, 100.0); - ReservoirOutlet2 = std::make_shared(0.3, 0.5, 10.0); + ReservoirOutlet2 = std::make_shared(0.3, 0.5, 10.0, 100.0); - ReservoirOutlet3 = std::make_shared(0.4, 0.6, 15.0); + ReservoirOutlet3 = std::make_shared(0.4, 0.6, 15.0, 100.0); ReservoirOutletsVector.push_back(*ReservoirOutlet1); @@ -88,11 +100,11 @@ void NonlinearReservoirKernelTest::setupMultipleOutletNonlinearReservoir() //Construct a reservoir with multiple outlets that are not ordered from lowest to highest activation threshold void NonlinearReservoirKernelTest::setupMultipleOutletOutOfOrderNonlinearReservoir() { - ReservoirOutlet1 = std::make_shared(0.2, 0.4, 4.0); + ReservoirOutlet1 = std::make_shared(0.2, 0.4, 4.0, 100.0); - ReservoirOutlet2 = std::make_shared(0.3, 0.5, 10.0); + ReservoirOutlet2 = std::make_shared(0.3, 0.5, 10.0, 100.0); - ReservoirOutlet3 = std::make_shared(0.4, 0.6, 15.0); + ReservoirOutlet3 = std::make_shared(0.4, 0.6, 15.0, 100.0); ReservoirOutletsVectorOutOfOrder.push_back(*ReservoirOutlet3); @@ -172,7 +184,7 @@ TEST_F(NonlinearReservoirKernelTest, TestRunOneOutletReservoirExceedMaxStorage) double excess; double final_storage; - in_flux_meters_per_second = 3.0; + in_flux_meters_per_second = 5.0; OneOutletReservoir->response_meters_per_second(in_flux_meters_per_second, 10, excess); @@ -207,6 +219,27 @@ TEST_F(NonlinearReservoirKernelTest, TestRunOneOutletReservoirFallsBelowMinimumS } +//Test Nonlinear Reservoir with one outlet where the calculated outlet velocity exceeds max velocity. +TEST_F(NonlinearReservoirKernelTest, TestRunOneOutletReservoirExceedMaxVelocity) +{ + double in_flux_meters_per_second; + double excess; + double final_storage; + + in_flux_meters_per_second = 80.0; + + OneOutletHighStorageReservoir->response_meters_per_second(in_flux_meters_per_second, 10, excess); + + final_storage = OneOutletHighStorageReservoir->get_storage_height_meters(); + + final_storage = round( final_storage * 100.0 ) / 100.0; + + EXPECT_DOUBLE_EQ (803.450, final_storage); + + ASSERT_TRUE(true); +} + + //Test Nonlinear Reservoir with multiple outlets TEST_F(NonlinearReservoirKernelTest, TestRunMultipleOutletReservoir) {