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;
}

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;
}

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;
}

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;
}

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;
}

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;
}

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)