Merge pull request #59 from tboudreaux/feature/mixedPolytrope

Static 3D FEM Polytropic Model
This commit is contained in:
2025-06-12 04:11:03 -04:00
committed by GitHub
5 changed files with 107 additions and 33 deletions

View File

@@ -1,3 +1,23 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include <cmath> #include <cmath>
#include <stdexcept> #include <stdexcept>
#include <array> #include <array>
@@ -52,7 +72,7 @@ The coefficients to the fit are from reaclib.jinaweb.org .
*/ */
namespace nnApprox8{ namespace serif::network::approx8{
// using namespace std; // using namespace std;
using namespace boost::numeric::odeint; using namespace boost::numeric::odeint;
@@ -225,7 +245,7 @@ namespace nnApprox8{
// a Jacobian matrix for implicit solvers // a Jacobian matrix for implicit solvers
void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) { void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) {
Constants& constants = Constants::getInstance(); serif::constant::Constants& constants = serif::constant::Constants::getInstance();
const double avo = constants.get("N_a").value; const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value; const double clight = constants.get("c").value;
// EOS // EOS
@@ -330,7 +350,7 @@ namespace nnApprox8{
} }
void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) { void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) {
Constants& constants = Constants::getInstance(); serif::constant::Constants& constants = serif::constant::Constants::getInstance();
const double avo = constants.get("N_a").value; const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value; const double clight = constants.get("c").value;
@@ -424,7 +444,7 @@ namespace nnApprox8{
dydt[Net::iener] = -enuc*avo*clight*clight; dydt[Net::iener] = -enuc*avo*clight*clight;
} }
nuclearNetwork::NetOut Approx8Network::evaluate(const nuclearNetwork::NetIn &netIn) { NetOut Approx8Network::evaluate(const NetIn &netIn) {
m_y = convert_netIn(netIn); m_y = convert_netIn(netIn);
m_tmax = netIn.tmax; m_tmax = netIn.tmax;
m_dt0 = netIn.dt0; m_dt0 = netIn.dt0;
@@ -468,7 +488,7 @@ namespace nnApprox8{
m_y[i] /= ysum; m_y[i] /= ysum;
} }
nuclearNetwork::NetOut netOut; NetOut netOut;
std::vector<double> outComposition; std::vector<double> outComposition;
outComposition.reserve(Net::nvar); outComposition.reserve(Net::nvar);
@@ -486,7 +506,7 @@ namespace nnApprox8{
m_stiff = stiff; m_stiff = stiff;
} }
vector_type Approx8Network::convert_netIn(const nuclearNetwork::NetIn &netIn) { vector_type Approx8Network::convert_netIn(const NetIn &netIn) {
if (netIn.composition.size() != Net::niso) { if (netIn.composition.size() != Net::niso) {
LOG_ERROR(m_logger, "Error: composition size mismatch in convert_netIn"); LOG_ERROR(m_logger, "Error: composition size mismatch in convert_netIn");
throw std::runtime_error("Error: composition size mismatch in convert_netIn"); throw std::runtime_error("Error: composition size mismatch in convert_netIn");

View File

@@ -1,14 +1,34 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Authors: Aaron Dotter, Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "network.h" #include "network.h"
#include "probe.h" #include "probe.h"
#include "quill/LogMacros.h" #include "quill/LogMacros.h"
namespace nuclearNetwork { namespace serif::network {
Network::Network() : Network::Network() :
m_config(Config::getInstance()), m_config(serif::config::Config::getInstance()),
m_logManager(Probe::LogManager::getInstance()), m_logManager(serif::probe::LogManager::getInstance()),
m_logger(m_logManager.getLogger("log")) { m_logger(m_logManager.getLogger("log")) {
} }
nuclearNetwork::NetOut nuclearNetwork::Network::evaluate(const NetIn &netIn) { NetOut Network::evaluate(const NetIn &netIn) {
// You can throw an exception here or log a warning if it should never be used. // You can throw an exception here or log a warning if it should never be used.
LOG_ERROR(m_logger, "nuclearNetwork::Network::evaluate() is not implemented"); LOG_ERROR(m_logger, "nuclearNetwork::Network::evaluate() is not implemented");
throw std::runtime_error("nuclearNetwork::Network::evaluate() is not implemented"); throw std::runtime_error("nuclearNetwork::Network::evaluate() is not implemented");

View File

@@ -1,5 +1,24 @@
#ifndef APPROX8_H /* ***********************************************************************
#define APPROX8_H //
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <array> #include <array>
@@ -34,7 +53,7 @@ typedef boost::numeric::ublas::matrix< double > matrix_type;
*/ */
typedef std::array<double,7> vec7; typedef std::array<double,7> vec7;
namespace nnApprox8{ namespace serif::network::approx8{
using namespace boost::numeric::odeint; using namespace boost::numeric::odeint;
@@ -275,14 +294,14 @@ namespace nnApprox8{
* @class Approx8Network * @class Approx8Network
* @brief Class for the Approx8 nuclear reaction network. * @brief Class for the Approx8 nuclear reaction network.
*/ */
class Approx8Network : public nuclearNetwork::Network { class Approx8Network : public Network {
public: public:
/** /**
* @brief Evaluates the nuclear network. * @brief Evaluates the nuclear network.
* @param netIn Input parameters for the network. * @param netIn Input parameters for the network.
* @return Output results from the network. * @return Output results from the network.
*/ */
virtual nuclearNetwork::NetOut evaluate(const nuclearNetwork::NetIn &netIn); virtual NetOut evaluate(const NetIn &netIn);
/** /**
* @brief Sets whether the solver should use a stiff method. * @brief Sets whether the solver should use a stiff method.
@@ -306,9 +325,7 @@ namespace nnApprox8{
* @param netIn Input parameters for the network. * @param netIn Input parameters for the network.
* @return Internal state vector. * @return Internal state vector.
*/ */
vector_type convert_netIn(const nuclearNetwork::NetIn &netIn); vector_type convert_netIn(const NetIn &netIn);
}; };
} // namespace nnApprox8 } // namespace nnApprox8
#endif

View File

@@ -1,5 +1,24 @@
#ifndef NETWORK_H /* ***********************************************************************
#define NETWORK_H //
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Authors: Emily Boudreaux, Aaron Dotter
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <vector> #include <vector>
@@ -7,7 +26,7 @@
#include "config.h" #include "config.h"
#include "quill/Logger.h" #include "quill/Logger.h"
namespace nuclearNetwork { namespace serif::network {
/** /**
* @struct NetIn * @struct NetIn
@@ -83,11 +102,9 @@ namespace nuclearNetwork {
virtual NetOut evaluate(const NetIn &netIn); virtual NetOut evaluate(const NetIn &netIn);
protected: protected:
Config& m_config; ///< Configuration instance serif::config::Config& m_config; ///< Configuration instance
Probe::LogManager& m_logManager; ///< Log manager instance serif::probe::LogManager& m_logManager; ///< Log manager instance
quill::Logger* m_logger; ///< Logger instance quill::Logger* m_logger; ///< Logger instance
}; };
} // namespace nuclearNetwork } // namespace nuclearNetwork
#endif // NETWORK_H

View File

@@ -14,13 +14,13 @@ class approx8Test : public ::testing::Test {};
* @brief Test the constructor of the Config class. * @brief Test the constructor of the Config class.
*/ */
TEST_F(approx8Test, constructor) { TEST_F(approx8Test, constructor) {
Config& config = Config::getInstance(); serif::config::Config& config = serif::config::Config::getInstance();
config.loadConfig(TEST_CONFIG); config.loadConfig(TEST_CONFIG);
EXPECT_NO_THROW(nnApprox8::Approx8Network()); EXPECT_NO_THROW(serif::network::approx8::Approx8Network());
} }
TEST_F(approx8Test, setStiff) { TEST_F(approx8Test, setStiff) {
nnApprox8::Approx8Network network; serif::network::approx8::Approx8Network network;
EXPECT_NO_THROW(network.setStiff(true)); EXPECT_NO_THROW(network.setStiff(true));
EXPECT_TRUE(network.isStiff()); EXPECT_TRUE(network.isStiff());
EXPECT_NO_THROW(network.setStiff(false)); EXPECT_NO_THROW(network.setStiff(false));
@@ -28,8 +28,8 @@ TEST_F(approx8Test, setStiff) {
} }
TEST_F(approx8Test, evaluate) { TEST_F(approx8Test, evaluate) {
nnApprox8::Approx8Network network; serif::network::approx8::Approx8Network network;
nuclearNetwork::NetIn netIn; serif::network::NetIn netIn;
std::vector<double> comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4}; std::vector<double> comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
@@ -41,10 +41,10 @@ TEST_F(approx8Test, evaluate) {
netIn.tmax = 3.15e17; netIn.tmax = 3.15e17;
netIn.dt0 = 1e12; netIn.dt0 = 1e12;
nuclearNetwork::NetOut netOut; serif::network::NetOut netOut;
EXPECT_NO_THROW(netOut = network.evaluate(netIn)); EXPECT_NO_THROW(netOut = network.evaluate(netIn));
EXPECT_DOUBLE_EQ(netOut.composition[nnApprox8::Net::ih1], 0.50166260916650918); EXPECT_DOUBLE_EQ(netOut.composition[serif::network::approx8::Net::ih1], 0.50166260916650918);
EXPECT_DOUBLE_EQ(netOut.composition[nnApprox8::Net::ihe4],0.48172270591286032); EXPECT_DOUBLE_EQ(netOut.composition[serif::network::approx8::Net::ihe4],0.48172270591286032);
EXPECT_DOUBLE_EQ(netOut.energy, 1.6433049870528356e+18); EXPECT_DOUBLE_EQ(netOut.energy, 1.6433049870528356e+18);
} }