@mainpage 4DSSE: A 4D Stellar Structure and Evolution Code

@section intro_sec 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.

@section build_sec 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.

@section usage_sec High-Level Usage Examples

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

@subsection usage_polysolver Solving for a Polytrope

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

```cpp
#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;
}
```

@subsection usage_composition Managing Chemical Compositions

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

```cpp
#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;
}
```

@subsection usage_network Nuclear Reaction Networks

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

```cpp
#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;
}
```

@subsection usage_constants Accessing Physical Constants

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

```cpp
#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;
}
```

@subsection usage_config Configuration Management

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

```cpp
#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;
}
```

@subsection usage_logging Logging

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

```cpp
#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;
}
```

@subsection usage_eos Equation of State (EOS)

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

```cpp
#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;
}
```

@subsection usage_meshio Mesh Handling

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

```cpp
#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;
}
```

@section modules_sec 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.

@section future_dev 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.

@section contact_sec 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)