CoolSolve is a parser, structural analyzer, and equation evaluator for the EES (Engineering Equation Solver) language, designed to be the foundation of an open-source equation solver that integrates with CoolProp for thermodynamic property calculations.

• EES Parser: Parses EES source code (.eescode files) into an Abstract Syntax Tree

  • Variables (scalars and arrays like P[1])
  • Equations with operators (+, -, *, /, ^)
  • Procedural Statements: Support for = and := in procedural blocks
  • Functions and Procedures: FUNCTION and PROCEDURE blocks with local scoping
  • Procedure Calls: CALL name(inputs : outputs) syntax
  • Conditional Logic: IF-THEN-ELSE support within procedural blocks
  • Comments ("...", {...}, //...)
  • Directives ($ifnot, $endif, etc.)
  • Units annotations ("[Pa]", "[kJ/kg]")
  • Function calls with named arguments (enthalpy(R134a, T=300, P=1E5))

• Structural Analysis: Decomposes the equation system into solvable blocks

  • Hopcroft-Karp algorithm for variable-equation matching
  • Tarjan's algorithm for Strongly Connected Components (SCCs)
  • Block decomposition (identifies algebraic loops)

• Automatic Differentiation: Forward-mode AD for efficient Jacobian computation

  • Full support for arithmetic operators (+, -, *, /, ^)
  • Transcendental functions (sin, cos, exp, log, sqrt, etc.)
  • Exact analytical derivatives (no finite differences)
  • Prepares for Newton-based numerical solvers

• Equation Evaluator: Evaluates residuals and Jacobians for each block

  • Block-level evaluation with external variable support
  • System-level orchestration for sequential block solving
  • CoolProp integration for thermodynamic property calculations
  • Automatic temperature conversion (Celsius ↔ Kelvin)
  • Configurable Solver Pipeline with multiple algorithms and execution modes
  • Explicit solve for size-1 blocks: Bypasses Newton entirely for structurally explicit assignments (one residual evaluation, no Jacobian)
  • Newton1D root-finder for size-1 implicit blocks with multi-probe exploration and bisection fallback
  • Newton + Line Search for fast convergence on well-conditioned blocks
  • Trust-Region Dogleg for robust convergence on stiff nonlinear blocks
  • Levenberg-Marquardt for improved convergence when initial guesses are poor
  • Partitioned Block Updates as a fallback for ill-conditioned algebraic loops
  • Parallel solver execution (multithreaded, first-to-converge wins)

• Output Formats:

  • JSON (for automated testing and integration)
  • EES-compatible residuals format
  • LaTeX equations
  • Evaluator report (block summary and state)
  • Solution/Initials files: Support for loading initial values and comparing solutions

• Debug Mode: Creates a comprehensive output folder with all analysis information

For details on current performance optimizations and the future roadmap for solver integration, see the Performance Improvement and Solver Integration Plan.

CoolSolve uses several file formats for input and verification:

• .eescode: The EES source code to be parsed and solved.
• .initials: Initial values for variables, used to seed the solver or evaluator. Format: variable=value (one per line).
• coolsolve.conf: Optional solver configuration. Place in the same folder as your .eescode file (not in subfolders). Format: key = value per line; lines starting with # are comments. Only the options you set override the defaults (see include/coolsolve/solver.h for SolverOptions). An example with all keys and comments is in examples/coolsolve.conf. In debug mode (-d), this file is copied into the debug folder. Key pipeline options:
  • solverPipeline: Comma-separated list of solvers to try (e.g. Newton, LM, TrustRegion, Partitioned)
  • pipelineMode: sequential (default) or parallel (first-to-converge wins)
  • enableTearing: When true, use structural tearing for blocks of size ≥ tearingMinBlockSize: a greedy feedback-vertex-set breaks the block into tear variables (solved with Newton) and an acyclic part (solved sequentially). Can help on stiff or ill-conditioned algebraic loops. See also tearingMaxIterations, tearingMinBlockSize, tearingInnerIterations.

• C++17 compatible compiler (GCC 7+, Clang 5+)
• CMake 3.14 or later
• Git (for fetching dependencies)

CoolSolve is a standalone library that automatically downloads all its dependencies using CMake's FetchContent mechanism:

Dependencies are cached in .fetchcontent_cache/ at the project root. This cache persists across build folder deletions, so you won't need to re-download dependencies if you clean and rebuild. The first cmake run will take several minutes to fetch CoolProp and its submodules; subsequent runs take only a few seconds.

mkdir build && cd build
cmake ..
make -j$(nproc)

This will build:

• coolsolve - The main executable
• coolsolve_tests - The test suite

Important: CoolSolve defaults to Release mode for optimal performance. CoolProp property calculations are computationally intensive, and Debug mode can be 10-50x slower than Release.

# Default: Release mode (recommended)
cmake ..

# Explicitly set Release mode
cmake -DCMAKE_BUILD_TYPE=Release ..

# Debug mode (for development only - very slow!)
cmake -DCMAKE_BUILD_TYPE=Debug ..

If you experience slow solve times, verify you're building in Release mode:

# Check current build type
grep CMAKE_BUILD_TYPE CMakeCache.txt

By default, CoolSolve fetches the latest master branch of CoolProp. To use a specific version:

# Use a specific tag or commit
cmake -DCOOLSOLVE_COOLPROP_TAG=v6.6.0 ..

# Or a specific commit hash
cmake -DCOOLSOLVE_COOLPROP_TAG=abc123def ..

# Clean build artifacts only (keeps dependency cache)
rm -rf build && mkdir build && cd build && cmake .. && make -j$(nproc)

# Full clean including dependency cache (will re-download everything)
rm -rf build .fetchcontent_cache

# Parse an EES file and output JSON analysis
./coolsolve input.eescode

# Output in EES residuals format
./coolsolve -f residuals input.eescode

# Output LaTeX equations
./coolsolve -f latex input.eescode

# Save output to a file
./coolsolve -o output.json input.eescode

Debug mode creates a folder containing all analysis information, useful for understanding the equation structure and debugging. For a complete guide on diagnosing and fixing solver failures (including the simplified-model workflow), see Debugging Models.

# Create debug output in <input>_coolsolve/ folder
./coolsolve -d input.eescode

# Specify custom debug output directory
./coolsolve -d my_debug_folder input.eescode

The debug folder contains:

To validate the structural analysis against EES output:

./coolsolve -c reference.residuals input.eescode

CoolProp's superancillary functions provide high-accuracy VLE (vapor-liquid equilibrium) calculations, especially near the critical point. However, they add computational overhead during solving.

# Default: superancillaries enabled 
./coolsolve model.eescode

# superancillaries disabled
./coolsolve --no-superancillary model.eescode

# Alternative: use environment variable
COOLPROP_ENABLE_SUPERANCILLARIES=false ./coolsolve model.eescode

cd build

# Run all unit tests
./coolsolve_tests

# Run comprehensive example file tests
./coolsolve_tests "[examples-comprehensive]"

Or using CTest:

cd build
ctest --output-on-failure

The comprehensive test runs all .eescode files in the examples/ folder and validates solutions against known expected values (with 1% tolerance). A detailed report is written to examples/test_examples.md.

CoolSolve/
├── CMakeLists.txt              # Build configuration (handles all dependencies)
├── README.md                   # This file
├── main.cpp                    # CLI entry point
├── .fetchcontent_cache/        # Cached dependencies (auto-generated, git-ignored)
├── include/coolsolve/
│   ├── ast.h                   # Abstract Syntax Tree definitions
│   ├── parser.h                # Parser interface
│   ├── ir.h                    # Intermediate Representation
│   ├── structural_analysis.h   # Analysis algorithms
│   ├── autodiff_node.h         # Forward-mode AD types and operations
│   ├── evaluator.h             # Block and system evaluators
│   └── solver.h                # Solver pipeline, Newton, TR, LM, Partitioned
├── src/
│   ├── parser.cpp              # EES parser implementation
│   ├── ir.cpp                  # IR building and LaTeX generation
│   ├── structural_analysis.cpp # Matching and SCC algorithms
│   ├── autodiff_node.cpp       # AD function implementations
│   ├── evaluator.cpp           # Evaluator implementations
│   └── solver.cpp              # All solver implementations + pipeline orchestrator
├── tests/
│   ├── test_parser.cpp         # Parser/IR unit tests (Catch2)
│   ├── test_evaluator.cpp      # AD/Evaluator unit tests (Catch2)
│   ├── test_newton.cpp         # Newton solver unit tests
│   ├── test_solver_pipeline.cpp # Pipeline, LM, config tests
│   ├── test_solver_integration.cpp # Full-system solver integration tests
│   └── test_examples.cpp       # Integration tests with example files
└── examples/                   # Example .eescode files for testing

The parser reads EES source code and builds an Abstract Syntax Tree (AST). It handles EES-specific syntax including:

• Case-insensitive keywords
• Multiple comment styles
• Thermodynamic function calls with named parameters
• Array subscript notation

The AST is transformed into an Intermediate Representation (IR) that:

• Extracts all variables from equations
• Builds the incidence matrix (which variables appear in which equations)
• Identifies thermodynamic function calls (first argument is fluid name, not a variable)

The equation system is analyzed using graph algorithms:

1. Matching: The Hopcroft-Karp algorithm finds a maximum bipartite matching, assigning each equation to one "output" variable it will determine.

2. SCC Detection: Tarjan's algorithm finds Strongly Connected Components in the dependency graph. Each SCC becomes a "block" that must be solved simultaneously.

3. Block Ordering: Blocks are topologically sorted so they can be solved in sequence.

Forward-mode automatic differentiation computes exact derivatives alongside values:

1. ADValue: A dual number storing both value and gradient (partial derivatives w.r.t. all block variables).

2. Propagation Rules: Each operation propagates gradients using the chain rule:

   • Addition: (x + y).grad = x.grad + y.grad
   • Multiplication: (x * y).grad = x.grad * y + y.grad * x
   • Power: (x^n).grad = n * x^(n-1) * x.grad
   • Functions: sin(x).grad = cos(x) * x.grad

3. Jacobian Construction: For each equation residual F_i, the gradient gives row i of the Jacobian matrix.

The evaluator system provides numerical computation:

1. ExpressionEvaluator: Traverses the AST and computes ADValues for any expression.

2. BlockEvaluator: Evaluates a block of equations, returning residuals F(x) and Jacobian J.

3. SystemEvaluator: Orchestrates evaluation across all blocks with proper variable handling.

CoolSolve solves each block using a strategy adapted to the block's size and structure, with several layers of fallback:

Block to solve
  ├── Size 1, explicit?  →  Direct evaluation (no iteration)
  ├── Size 1, implicit?  →  Newton1D solver (see below)
  ├── Size ≥ tearingMinBlockSize & tearing enabled?
  │     →  Structural tearing (FVS + acyclic solve + Newton on tears)
  │     └── On failure: restore initial guess, fall through ↓
  └── Solver pipeline (configurable fallback chain)
        →  Newton → TrustRegion → LM → Partitioned  (default, sequential)

The pipeline is configured via coolsolve.conf (see below) or programmatically through SolverOptions::solverPipeline and SolverOptions::pipelineMode.

0. Newton1D (automatic for size-1 implicit blocks)

   • Specialized root-finding solver for single-equation blocks where the unknown cannot be isolated symbolically.
   • Phase 1 — Trust-region Newton: Standard Newton steps with adaptive radius limiting and sign-change detection for bisection fallback.
   • Phase 2 — Multi-probe exploration: If Phase 1 stalls (e.g. initial guess far from root), evaluates the residual at ~900 probe points: log-spaced values across ±1e8, plus values near every external variable (×0.5, ×0.9, …, ×2.0). This finds narrow sign-change regions even when they occur near poles in the residual. All sign changes are scored by midpoint residual to prefer true roots over poles.
   • Phase 3 — Bisection + Newton hybrid: Once a bracket is found, alternates bisection and Newton steps within the bracket for fast, guaranteed convergence.
   • Phase 4 — Final Newton polish: A short Newton loop from the best point found, with relaxed tolerance acceptance.
   • Falls through to the standard pipeline if all phases fail.

1. Newton + Line Search (Newton)

   • Solves J(x) * dx = -F(x) and applies backtracking to ensure descent.
   • Efficient when the Jacobian is well-conditioned and the model is smooth.

2. Trust-Region Dogleg (TrustRegion)

   • Uses a dogleg step that blends steepest descent with the Newton step to keep updates inside a safe radius.
   • Helps avoid oversized steps that drive thermodynamic calls into invalid regions (e.g., non-physical pressure/temperature).

3. Levenberg-Marquardt (LM or LevenbergMarquardt)

   • Solves (J^T J + λ D) dx = -J^T F with adaptive damping parameter λ.
   • When λ is large → gradient descent (safe, slow); when λ is small → Gauss-Newton (fast, quadratic convergence near solution).
   • Particularly effective when the initial guess is far from the solution, because the damping prevents oversized steps that would cause divergence.

4. Partitioned Block Updates (Partitioned)

   • Uses the equation-to-output-variable mapping from structural matching to apply per-variable diagonal updates inside a block: x_i <- x_i - w * F_i / (dF_i/dx_i).
   • This mimics a DAE-style "tear" without changing the block structure: each equation directly updates its matched variable, reducing coupling and improving stability in stiff or highly nonlinear loops.
   • Designed as a last-resort stabilizer when full Newton steps are unreliable.

5. Structural Tearing (option enableTearing)

   • When enabled, blocks of size ≥ tearingMinBlockSize are first solved via equation tearing: a greedy feedback vertex set (FVS) is computed so that removing the corresponding equations (and their output variables) makes the block acyclic. The acyclic part is then solved sequentially (one equation, one unknown per step), and the tear variables are updated with Newton on the tear residuals. This reduces the simultaneous system to the tear set only and can improve robustness on stiff or ill-conditioned loops.
   • Config: enableTearing = true, tearingMaxIterations, tearingMinBlockSize, tearingInnerIterations. In debug mode (-d), a tearing.md file lists tear sets and acyclic order per block.

Models operating near phase boundaries or critical points (e.g. supercritical CO2) can produce unphysical trial points during Newton iteration. Two layers of protection are applied:

• Input clamping: Pressure, temperature, and density are clamped to physically valid floors ($P \geq 1000$ Pa, $T \geq 50$ K, $\rho \geq 10^{-4}$) before CoolProp calls, preventing the most common crash-inducing inputs.
• Penalty-based error handling: If CoolProp still throws, a finite penalty value (10⁴) is returned instead of propagating the exception, keeping the residual landscape smooth for TR and LM solvers.

• Sequential (default): Solvers are tried one after another. Each subsequent solver warm-starts from the best solution found so far (rather than resetting to the original initial guess). If a full pipeline pass reduces the residual by ≥5% without converging, the pipeline restarts from the best point for up to 10 rounds. The Partitioned solver is automatically skipped in later rounds if it worsens the solution.
• Parallel: All solvers are launched concurrently in separate threads. The first solver to converge wins and its solution is used. This can save time when it is unclear which algorithm will work best for a given block.

On failure, the error message reports both the initial and best-achieved residual norms with the percentage of reduction, helping diagnose whether the problem is a poor initial guess or a genuinely unsolvable system.

solverPipeline = Newton, TrustRegion, LevenbergMarquardt, Partitioned
pipelineMode = sequential

To use only Levenberg-Marquardt (no fallback):

solverPipeline = LM

solverPipeline = Newton, TrustRegion, LM
pipelineMode = parallel

The analysis results can be exported in various formats for:

• Integration with numerical solvers (JSON)
• Documentation (LaTeX)

Given this EES code:

T_in = 25            // Temperature in Celsius
P = 101325           // Pressure in Pa
h = enthalpy(Water, T=T_in, P=P)
s = entropy(Water, T=T_in, P=P)

CoolSolve will:

1. Parse 4 equations
2. Identify 4 variables: T_in, P, h, s
3. Create 4 blocks (all explicit, no algebraic loops)
4. Determine solution order: T_in → P → h, s
5. Evaluate thermodynamic properties using CoolProp (with automatic °C→K conversion)

CoolSolve integrates with CoolProp for thermodynamic property calculations:

• Supported functions: enthalpy, entropy, density, pressure, temperature, quality, cp, cv
• Supported input pairs: T,P, T,x, P,h, P,s, H,P, and others
• Temperature units: Automatically converts Celsius (EES convention) to Kelvin (CoolProp)
• Derivatives: Computed using finite differences from CoolProp's PropsSI function
• Fluids: All CoolProp fluids are supported (Water, R134a, Nitrogen, etc.)

CoolSolve assumes the following default units:

• Temperature: °C (automatically converted to K for CoolProp calls)
• Pressure: Pa
• Energy: J
• Mass: kg
• Force: N
• Length: m

Built-in constants like g# use these SI units (e.g., g# = 9.80665 m/s^2).

Example usage in EES code:

T_ev = -10           // Celsius
P_ev = pressure(R134a, T=T_ev, x=1)   // Saturation pressure
h_1 = enthalpy(R134a, P=P_ev, T=T_1)  // Enthalpy at state 1

CoolSolve supports user-defined functions and procedures for modular code:

PROCEDURE single_phase_HX(cf$, hf$, t_cf_su : Q_dot)
    "Procedural block with local scope"
    cp_cf = specheat(cf$, T=t_cf_su, P=101325)
    Q_dot := alpha * cp_cf * (t_hf_su - t_cf_su)
END

CALL single_phase_HX('Air_ha', 'Water', 20 : Q_total)

• Local Scope: Variables inside functions and procedures are local and do not interfere with the main program.
• Assignment: Use := for assignment inside procedural blocks (though = is also supported for compatibility).
• Control Flow: IF-THEN-ELSE statements are supported within these blocks.
• Automatic Differentiation: CoolSolve automatically propagates derivatives through procedural calls, ensuring accurate Jacobians.

Important: EES and CoolProp use different reference states for thermodynamic properties:

This means:

• Enthalpy values from EES may differ from CoolProp by ~100-200 kJ/kg
• Entropy values from EES may differ from CoolProp by ~0.8-1 kJ/kg/K
• Differences (e.g., h2 - h1) are the same in both reference states
• When loading initial values from .initials files, some CoolProp function calls may fail if the enthalpy/entropy is outside the expected range for the given pressure

For new equation systems, use CoolProp-computed values for initial guesses to ensure consistency.

The next steps in the implementation plan include:

• KINSOL (SUNDIALS) integration: For large-scale nonlinear systems requiring robust preconditioning
• alternative solvers: https://www.osti.gov/servlets/purl/876345
• Homotopy / Continuation methods: Gradually transform from an easy problem to the target for highly nonlinear systems
• Analytical Derivatives: Use CoolProp's AbstractState::first_partial_deriv for exact derivatives
• Profiling and optimization: Improve the solving time