Loading...
Searching...
No Matches
4DSSE: A 4D Stellar Structure and Evolution Code

# Introduction

Welcome to the documentation for 4DSSE (4D Stellar Structure and Evolution), a new code designed for simulating stellar phenomena in three spatial dimensions plus time. This project is currently under active development.

The primary goal of 4DSSE is to provide a flexible and extensible framework for advanced stellar modeling, incorporating modern numerical techniques and physics modules.

# Building 4DSSE

The project uses Meson as its build system. MFEM is a core dependency and will be automatically downloaded and built if not found.

Prerequisites:

- A C++ compiler (supporting C++17 or later)
- Meson (Install via pip: `pip install meson`)
- Python 3

Build Steps:

1. Clone the repository (if you haven't already): `bash git clone <repository-url> cd 4dsse`
2. Using the `mk` script (recommended for ease of use):
   - To build the project: `bash ./mk`
   - To build without tests: `bash ./mk --noTest`
3. Using the `4DSSEConsole.sh` script: This script provides a simple interface for building and debugging. `bash ./4DSSEConsole.sh` Follow the on-screen prompts.
4. Manual Meson build:
   - Setup the build directory (e.g., `build`): `bash meson setup build`
   - Compile the project: `bash meson compile -C build`
   - To run tests (if built with tests enabled): `bash meson test -C build`

The compiled executables and libraries will typically be found in the `build` directory.

# High-Level Usage Examples

Below are some high-level examples of how to use key components of 4DSSE.

## Solving for a Polytrope

The `PolySolver` class handles the setup and solution of the Lane-Emden equation for polytropic stellar models.

#include "polySolver.h"
#include "config.h" // For global configuration
#include "probe.h" // For logging
int main() {
// Initialize configuration and logging
Config::getInstance().loadConfig("path/to/your/config.yaml");
Probe::LogManager::getInstance().getLogger("main_log"); // Initialize a logger
try {
// Create a PolySolver for polytropic index n=1.5 and FE order 2
PolySolver solver(1.5, 2);
solver.solve(); // Solve the system
// Access the solution
mfem::GridFunction& theta_solution = solver.getSolution();
// ... process or visualize theta_solution ...
} catch (const std::exception& e) {
std::cerr << "An error occurred: " << e.what() << std::endl;
return 1;
}
return 0;
}
main
int main(int argv, char *argc[])
Definition comp.cpp:5
config.h
polySolver.h
probe.h

## Managing Chemical Compositions

The `Composition` class allows for defining and managing chemical compositions.

#include "composition.h"
#include <iostream>
#include <vector>
int main() {
try {
// Define symbols and their mass fractions
std::vector<std::string> symbols = {"H-1", "He-4"}; // Use specific isotopes
std::vector<double> mass_fractions = {0.75, 0.25};
// Create and finalize the composition
composition::Composition comp(symbols, mass_fractions);
// Get mass fraction of a specific element
std::cout << "Mass fraction of H-1: " << comp.getMassFraction("H-1") << std::endl;
// Get global properties
auto global_props = comp.getComposition().second;
std::cout << "Mean particle mass: " << global_props.meanParticleMass << std::endl;
} catch (const std::exception& e) {
std::cerr << "Composition error: " << e.what() << std::endl;
}
return 0;
}
composition.h

## Nuclear Reaction Networks

The `Network` and `Approx8Network` classes provide interfaces for nuclear reaction network calculations.

#include "network.h"
#include "approx8.h" // Specific network implementation
#include <iostream>
#include <vector>
int main() {
nnApprox8::Approx8Network approx8_net; // Using the Approx8 network
approx8_net.setStiff(true); // Example: use stiff solver
nuclearNetwork::NetIn input;
input.composition = {0.7, 0.0, 0.28, 0.01, 0.005, 0.004, 0.0005, 0.0005}; // H1, He3, He4, C12, N14, O16, Ne20, Mg24
input.temperature = 1.5e7; // K
input.density = 150.0; // g/cm^3
input.tmax = 1.0e10; // s
input.dt0 = 1.0e6; // s
// input.energy can also be set if needed
try {
nuclearNetwork::NetOut output = approx8_net.evaluate(input);
std::cout << "Number of steps: " << output.num_steps << std::endl;
std::cout << "Final H-1 mass fraction (approx): " << output.composition[nnApprox8::Net::ih1] << std::endl;
} catch (const std::exception& e) {
std::cerr << "Network evaluation error: " << e.what() << std::endl;
}
return 0;
}
approx8.h
Header file for the Approx8 nuclear reaction network.
network.h

## Accessing Physical Constants

The `Constants` singleton provides access to a database of physical constants.

#include "const.h"
#include <iostream>
int main() {
Constants& consts = Constants::getInstance();
if (consts.isLoaded()) {
Constant G = consts.get("Gravitational constant");
std::cout << G.name << ": " << G.value << " " << G.unit << std::endl;
Constant c = consts["Speed of light in vacuum"]; // Can also use operator[]
std::cout << c.name << ": " << c.value << " " << c.unit << std::endl;
} else {
std::cerr << "Failed to load constants." << std::endl;
}
return 0;
}
const.h

## Configuration Management

The `Config` singleton manages settings from a YAML configuration file.

#include "config.h"
#include <iostream>
int main() {
Config& config = Config::getInstance();
if (config.loadConfig("path/to/your/config.yaml")) {
// Get a string value, with a default if not found
std::string outputPath = config.get<std::string>("Output:Path", "./output/");
std::cout << "Output path: " << outputPath << std::endl;
// Get an integer value
int maxIter = config.get<int>("Solver:MaxIterations", 100);
std::cout << "Max iterations: " << maxIter << std::endl;
} else {
std::cerr << "Failed to load configuration." << std::endl;
}
return 0;
}

## Logging

The `Probe::LogManager` provides a way to manage and use loggers.

#include "probe.h"
#include "config.h" // Often used to configure logging
int main() {
// Assuming config is loaded and might define log file, level etc.
// Config::getInstance().loadConfig("config.yaml");
Probe::LogManager& logManager = Probe::LogManager::getInstance();
quill::Logger* mainLogger = logManager.getLogger("main_app_log"); // Get or create logger
// Example: Create a new file logger if not configured through a central mechanism
// quill::Logger* fileLogger = logManager.newFileLogger("app_trace.log", "trace_log");
LOG_INFO(mainLogger, "Application started. Version: {}", "1.0.0");
// ... application logic ...
LOG_ERROR(mainLogger, "An unexpected error occurred in module X.");
return 0;
}

## Equation of State (EOS)

The `EosIO` class loads EOS tables, and the `helmholtz` namespace provides functions to use them, for example, the Helmholtz EOS.

#include "eosIO.h"
#include "helm.h"
#include <iostream>
int main() {
try {
// Load the Helmholtz EOS table
EosIO helm_eos_io("path/to/helm_table.dat"); // Replace with actual path
EOSTable& table_variant = helm_eos_io.getTable();
// Assuming it's a HELMTable, get it (add error checking in real code)
auto* helm_table_ptr = std::get_if<std::unique_ptr<helmholtz::HELMTable>>(&table_variant);
if (!helm_table_ptr || !(*helm_table_ptr) || !(*helm_table_ptr)->loaded) {
std::cerr << "Failed to load or access HELM table." << std::endl;
return 1;
}
const helmholtz::HELMTable& table = **helm_table_ptr;
// Define input conditions
helmholtz::EOSInput input;
input.T = 1.0e7; // Temperature in K
input.rho = 100.0; // Density in g/cm^3
input.abar = 1.0; // Mean atomic mass (e.g., for pure hydrogen)
input.zbar = 1.0; // Mean atomic number (e.g., for pure hydrogen)
// Get EOS results
helmholtz::EOS results = helmholtz::get_helm_EOS(input, table);
std::cout << "Total Pressure (Ptot): " << results.ptot << " dyne/cm^2" << std::endl;
std::cout << "Total Energy (Etot): " << results.etot << " erg/g" << std::endl;
} catch (const std::exception& e) {
std::cerr << "EOS error: " << e.what() << std::endl;
return 1;
}
return 0;
}
EOSio.h
helm.h

## Mesh Handling

The `MeshIO` class facilitates loading and managing computational meshes.

#include "meshIO.h"
#include "mfem.hpp" // For mfem::Mesh
#include <iostream>
int main() {
try {
// Initialize MeshIO with a mesh file and a scale factor
MeshIO mesh_handler("path/to/your/mesh.msh", 1.0); // Replace with actual path
if (mesh_handler.IsLoaded()) {
mfem::Mesh& mesh = mesh_handler.GetMesh();
std::cout << "Mesh loaded successfully with " << mesh.GetNE() << " elements." << std::endl;
// Optionally, rescale the mesh
// mesh_handler.LinearRescale(2.0);
// std::cout << "Mesh rescaled. New bounding box: ";
// mfem::Vector min, max;
// mesh.GetBoundingBox(min, max);
// min.Print(std::cout); max.Print(std::cout);
} else {
std::cerr << "Failed to load mesh." << std::endl;
}
} catch (const std::exception& e) {
std::cerr << "MeshIO error: " << e.what() << std::endl;
}
return 0;
}
meshIO.h

# Key Modules and Components

4DSSE is organized into several key modules:

- Polytrope Solver (`polySolver.h`, `polytropeOperator.h`): Provides tools to solve the Lane-Emden equation for polytropic stellar structures using a mixed finite element method. `PolytropeOperator` defines the nonlinear system and its Jacobian, while `PolySolver` orchestrates the solution process. The `SchurCompliment` and `GMRESInverter` classes are helper components for the linear algebra involved.
- Equation of State (EOS) (`helm.h`, `eosIO.h`): Manages Equation of State data. `helm.h` provides an implementation of the Helmholtz EOS (Timmes & Swesty 2000), including structures for table data (`HELMTable`), input parameters (`EOSInput`), and output results (`EOS`). It also defines functions for reading tables and calculating EOS quantities. `eosIO.h` provides the `EosIO` class for loading EOS tables from files, currently supporting the HELM table format.
- Chemical Composition (`composition.h`, `atomicSpecies.h`): Manages chemical compositions, allowing representation in mass or number fractions. It interfaces with `atomicSpecies.h` which provides a database of atomic species properties (based on AME2020).
- Nuclear Reaction Networks (`network.h`, `approx8.h`): Defines a base `Network` class for nuclear reaction network calculations. `approx8.h` provides a specific implementation, `Approx8Network`, for an 8-isotope network (H1, He3, He4, C12, N14, O16, Ne20, Mg24) based on Frank Timmes' "aprox8". It includes functions for individual reaction rates and uses Boost.Numeric.Odeint for solving the ODE system.
- Physical Constants (`const.h`): A singleton class `Constants` that loads and provides access to a wide range of physical constants with their values, uncertainties, units, and references.
- Configuration Management (`config.h`): A singleton class `Config` for loading and accessing application settings from YAML configuration files.
- Probing and Logging (`probe.h`): The `Probe` namespace offers utility functions for debugging, such as GLVis visualization (`glVisView`), and a `LogManager` for handling application-wide logging using the Quill library.
- Mesh I/O (`meshIO.h`): The `MeshIO` class handles loading and basic manipulation (e.g., scaling) of computational meshes using MFEM's `mfem::Mesh`. It ensures that meshes are correctly loaded and accessible.
- Integrators (`integrators.h`): (Details inferred) Likely contains custom MFEM integrators or coefficients used in the finite element formulations.
- Custom Types (`4DSTARTypes.h`): Defines project-specific data type aliases within the `SSE` namespace, primarily for simplifying common `std::pair` combinations involving `mfem::Array<int>` and `mfem::Array<double>`. These include `SSE::MFEMArrayPair` and `SSE::MFEMArrayPairSet`, often used for managing collections of MFEM degree-of-freedom lists and their corresponding values, especially for boundary conditions.

# Future Development

Future work will focus on expanding the physics modules (e.g., equation of state, opacity), improving numerical solvers, and enhancing the parallelization capabilities for large-scale simulations.

# Contact and Contributions

For questions, bug reports, or contributions, please refer to the project's repository or contact the development team. (Details to be added)