Skip to content

Repository files navigation

DiffusionEnv: High-Performance Finite Element Environment for Reinforcement Learning

A sophisticated Python environment that combines the numerical power of deal.II finite element computations with the flexibility of Python-based reinforcement learning frameworks. This package provides a high-performance solver for transient diffusion equations with multiple geometry support, wrapped in a clean Python interface optimized for RL research.

Table of Contents

Overview

DiffusionEnv solves the transient diffusion equation:

∂u/∂t = K∇²u

Where:

  • u(x,y,t) is the field variable (e.g., temperature, concentration)
  • K is the diffusion coefficient
  • The domain supports multiple geometry types with customizable dimensions and zero Dirichlet boundary conditions

New in v1.1: The environment now supports six different geometry types including squares, rectangles, circles, annuli, L-shaped domains, and quarter circles, each with customizable dimensions. This significantly expands the range of research scenarios and allows for testing algorithm robustness across different domain shapes.

This environment is specifically designed for reinforcement learning applications where agents need to:

  • Learn optimal control of PDE systems across different domain geometries
  • Understand spatio-temporal dynamics in various geometric contexts
  • Work with high-dimensional state spaces derived from physics simulations
  • Explore parameter optimization in computational physics with domain-specific challenges
  • Test algorithm robustness on non-rectangular and non-convex domains

Key Design Principles

  • Performance: C++ finite element core with minimal Python overhead
  • Geometric Flexibility: Support for multiple domain shapes with customizable dimensions
  • Compatibility: Works with standard RL frameworks (Stable-Baselines3, Ray RLlib, custom PyTorch)
  • Scientific Rigor: Built on deal.II's proven finite element implementation with proper curved boundary handling
  • Extensibility: Modular architecture for adding new physics, boundary conditions, or geometry types

Features

Core Functionality

  • High-Performance Finite Element Solver: Built on deal.II 9.5+ with optimized sparse linear algebra
  • Multiple Geometry Types: Support for squares, rectangles, circles, annuli, L-shaped domains, and quarter circles
  • Customizable Domain Dimensions: Full control over geometry parameters (sizes, radii, aspect ratios)
  • Flexible State Representations: Choose between compact physical quantities or full spatial solution fields
  • Configurable Physics Parameters: Adjustable diffusion coefficients, time steps, and simulation duration
  • Multiple Initial Conditions: Python function interface for arbitrary initial condition specification
  • Real-Time Data Access: Extract solution data, mesh geometry, and physical quantities at any time step
  • Geometry Introspection: Query domain bounds, geometry description, and coordinate systems
  • Visualization Support: Optional VTU output for detailed analysis in VisIt/ParaView

Geometry Support

  • Hyper Cube: Square domains with customizable side length
  • Hyper Rectangle: Rectangular domains with custom width and height
  • Hyper Ball: Circular domains with custom radius and center
  • Hyper Shell: Annular domains with custom inner and outer radii
  • L-Shaped: L-shaped domains with custom size (ideal for testing non-convex domains)
  • Quarter Hyper Ball: Quarter-circle domains with custom radius

RL Integration Features

  • Standard RL Interface: Reset/step paradigm compatible with OpenAI Gym conventions
  • Domain-Aware State Representation: Geometry information available for state augmentation
  • Batch Processing: Efficient repeated simulations for training across different geometries
  • Memory Efficient: Pre-compiled finite element structures reused across episodes
  • Configurable Rewards: Framework for physics-based and geometry-aware reward function design
  • Multi-Scale Analysis: Support for different mesh refinement levels across all geometries

Technical Features

  • CMake Integration: Leverages deal.II's robust build system
  • Backward Compatibility: Existing code continues to work unchanged
  • Python 3.7+ Support: Modern Python integration with type hints
  • NumPy Compatibility: All data returned as NumPy arrays for seamless ML integration
  • Error Handling: Comprehensive error reporting for debugging geometry configurations
  • Cross-Platform: Linux, macOS support (Windows with WSL)

Requirements

System Requirements

  • Operating System: Linux (recommended), macOS, Windows with WSL
  • Compiler: GCC 7+ or Clang 5+ with C++17 support
  • Memory: Minimum 4GB RAM (8GB+ recommended for larger problems)
  • CMake: Version 3.12 or higher

Software Dependencies

Required Dependencies

  1. deal.II 9.0+

    • High-performance finite element library
    • Must be compiled with shared library support
    • Installation guide: deal.II Documentation
  2. Python 3.7+

    • NumPy 1.15+
    • pybind11 2.6+

Optional Dependencies

  • Matplotlib: For visualization examples
  • SciPy: For advanced interpolation and analysis
  • Jupyter: For interactive examples
  • PyTorch/TensorFlow: For RL framework integration

Hardware Recommendations

Use Case CPU Memory Storage
Development 4+ cores 8GB 10GB
Research 8+ cores 16GB 50GB
Production 16+ cores 32GB 100GB

Installation

1. Install deal.II

Option A: Package Manager (Ubuntu/Debian)

sudo apt-get update
sudo apt-get install libdeal.ii-dev

Option B: From Source (Recommended for HPC)

# Download and compile deal.II
wget https://github.com/dealii/dealii/releases/download/v9.5.2/dealii-9.5.2.tar.gz
tar xzf dealii-9.5.2.tar.gz
cd dealii-9.5.2

mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=/path/to/dealii-install ..
make -j$(nproc)
make install

2. Set Environment Variables

# Add to your ~/.bashrc or ~/.zshrc
export DEAL_II_DIR=/path/to/dealii-install
export PATH=$DEAL_II_DIR/bin:$PATH

3. Create Python Environment

# Create virtual environment
python3 -m venv diffusion_env
source diffusion_env/bin/activate  # On Windows: diffusion_env\Scripts\activate

# Install Python dependencies
pip install --upgrade pip
pip install numpy matplotlib pybind11

4. Build DiffusionEnv

# Clone or download this repository
cd /path/to/diffusion_project

# Verify deal.II is found
echo $DEAL_II_DIR
deal.II-config --version  # Should work without errors

# Build the extension
python setup.py build_ext --inplace

# Test installation
python -c "import diffusion_env; print('Success!')"

Verification Script

# Run comprehensive verification including geometry features
python -c "
import diffusion_env
import numpy as np

# Test basic functionality
env = diffusion_env.DiffusionEnvironment(refinement_level=3)
print(f'✓ Environment created with {env.get_num_dofs()} DOF')

# Test geometry configuration
circle_config = diffusion_env.GeometryConfig.hyper_ball(radius=1.5)
circle_env = diffusion_env.DiffusionEnvironment(circle_config, refinement_level=3)
print(f'✓ Circular environment: {circle_env.get_geometry_description()}')

env.reset(lambda x, y: np.sin(np.pi * x) * np.sin(np.pi * y))
print('✓ Initial condition set')

env.step()
solution = env.get_solution_data()
print(f'✓ Solution extracted: {len(solution)} values')

# Test convenience functions
square_env = diffusion_env.create_square_env(size=2.0, refinement=3)
print(f'✓ Convenience function: {square_env.get_geometry_description()}')

print('Installation with geometry support verified successfully!')
"

Quick Start

Basic Usage (Backward Compatible)

import diffusion_env
import numpy as np

# Create environment with default unit square geometry
env = diffusion_env.DiffusionEnvironment(
    refinement_level=4,    # Mesh resolution (3-6 typical)
    diffusion_coeff=0.1,   # Physical parameter K
    dt=0.01,              # Time step size
    final_time=1.0        # Total simulation time
)

print(f"Environment has {env.get_num_dofs()} degrees of freedom")
print(f"Geometry: {env.get_geometry_description()}")

# Define initial condition
def gaussian_heat_source(x, y):
    """Gaussian heat source centered at (0.3, 0.7)"""
    return np.exp(-20 * ((x - 0.3)**2 + (y - 0.7)**2))

# Run simulation
env.reset(gaussian_heat_source)

while not env.is_done():
    env.step()
    
    # Get current state information
    solution = env.get_solution_data()      # Full spatial solution
    quantities = env.get_physical_quantities()  # Compact summary
    
    print(f"Time: {env.get_time():.3f}, "
          f"Energy: {quantities[0]:.6f}, "
          f"Max value: {quantities[1]:.6f}")

print("Simulation completed!")

New Geometry Features Usage

import diffusion_env
import numpy as np

# Example 1: Circular domain
circle_config = diffusion_env.GeometryConfig.hyper_ball(radius=2.0)
circle_env = diffusion_env.DiffusionEnvironment(circle_config, refinement_level=4)

print(f"Circle environment: {circle_env.get_geometry_description()}")
print(f"Domain bounds: {circle_env.get_domain_bounds()}")  # [x_min, x_max, y_min, y_max]

# Define initial condition appropriate for circular domain
def centered_gaussian(x, y):
    """Gaussian centered at origin for circular domain"""
    return np.exp(-(x**2 + y**2) / 0.5)

circle_env.reset(centered_gaussian)

# Example 2: Rectangular domain
rect_config = diffusion_env.GeometryConfig.hyper_rectangle(width=3.0, height=1.5)
rect_env = diffusion_env.DiffusionEnvironment(rect_config, refinement_level=4)

# Example 3: L-shaped domain (non-convex)
l_config = diffusion_env.GeometryConfig.l_shaped(domain_size=2.0)
l_env = diffusion_env.DiffusionEnvironment(l_config, refinement_level=4)

print(f"L-shaped environment: {l_env.get_geometry_description()}")

# Example 4: Annular domain (with hole in center)
annulus_config = diffusion_env.GeometryConfig.hyper_shell(inner_radius=0.5, outer_radius=1.5)
annulus_env = diffusion_env.DiffusionEnvironment(annulus_config, refinement_level=4)

# Example 5: Using convenience functions
square_env = diffusion_env.create_square_env(size=2.0, refinement=4)
circle_env_easy = diffusion_env.create_circle_env(radius=1.5, refinement=4)
rect_env_easy = diffusion_env.create_rectangle_env(width=3.0, height=2.0, refinement=4)

Geometry Types

Supported Geometries

1. Hyper Cube (Square)

config = diffusion_env.GeometryConfig.hyper_cube(side_length=2.0)
# Creates a square domain [0, 2] × [0, 2]

2. Hyper Rectangle

config = diffusion_env.GeometryConfig.hyper_rectangle(width=3.0, height=1.5)
# Creates a rectangular domain [0, 3] × [0, 1.5]

3. Hyper Ball (Circle)

config = diffusion_env.GeometryConfig.hyper_ball(radius=1.5, center=(0.0, 0.0))
# Creates a circular domain centered at origin with radius 1.5

4. Hyper Shell (Annulus)

config = diffusion_env.GeometryConfig.hyper_shell(inner_radius=0.5, outer_radius=1.2)
# Creates an annular domain with hole in center

5. L-Shaped Domain

config = diffusion_env.GeometryConfig.l_shaped(domain_size=1.5)
# Creates an L-shaped non-convex domain

6. Quarter Hyper Ball

config = diffusion_env.GeometryConfig.quarter_hyper_ball(radius=1.0)
# Creates a quarter-circle domain

Coordinate Systems

Different geometries use different coordinate systems:

Geometry Coordinate Range Center Location
Hyper Cube [0, size] × [0, size] (size/2, size/2)
Hyper Rectangle [0, width] × [0, height] (width/2, height/2)
Hyper Ball ≈ [-radius, radius] × [-radius, radius] (0, 0) by default
Hyper Shell Similar to hyper ball with hole (0, 0) by default
L-Shaped Complex, use get_domain_bounds() Variable
Quarter Ball First quadrant of circle Variable

Important: Always use env.get_domain_bounds() to get exact coordinate ranges for any geometry.

Choosing Geometry for Research

Geometry Best For Challenges
Square/Rectangle Algorithm development, baseline testing Simple, predictable
Circle Testing curved boundary handling No corners, symmetric
Annulus Obstacle avoidance, barrier problems Interior boundaries
L-Shaped Non-convex domains, corner effects Re-entrant corners, complex shape
Quarter Circle Reduced state space, rapid prototyping Boundary interactions

API Reference

GeometryConfig Class

class GeometryConfig:
    """Configuration for domain geometry and dimensions."""
    
    # Static factory methods (recommended approach)
    @staticmethod
    def hyper_cube(side_length: float = 1.0) -> GeometryConfig
        """Create square domain configuration."""
    
    @staticmethod
    def hyper_rectangle(width: float, height: float) -> GeometryConfig
        """Create rectangular domain configuration."""
    
    @staticmethod  
    def hyper_ball(radius: float, center: Tuple[float, float] = (0, 0)) -> GeometryConfig
        """Create circular domain configuration."""
    
    @staticmethod
    def hyper_shell(inner_radius: float, outer_radius: float, 
                   center: Tuple[float, float] = (0, 0)) -> GeometryConfig
        """Create annular domain configuration."""
    
    @staticmethod
    def l_shaped(domain_size: float = 1.0) -> GeometryConfig
        """Create L-shaped domain configuration."""
    
    @staticmethod
    def quarter_hyper_ball(radius: float, center: Tuple[float, float] = (0, 0)) -> GeometryConfig
        """Create quarter-circle domain configuration."""

DiffusionEnvironment Class

Constructor

DiffusionEnvironment(
    geometry_config: GeometryConfig = GeometryConfig.hyper_cube(),
    refinement_level: int = 4,
    diffusion_coeff: float = 0.1,
    dt: float = 0.01,
    final_time: float = 1.0
)

New Parameters:

  • geometry_config: Geometry configuration specifying domain type and dimensions

Backward Compatibility Constructor:

DiffusionEnvironment(
    refinement_level: int = 4,
    diffusion_coeff: float = 0.1,
    dt: float = 0.01,
    final_time: float = 1.0
)
# Creates unit hyper cube for backward compatibility

Mesh Size Reference:

Refinement Level Typical DOF (Square) Typical DOF (Circle) Typical Use
3 81 ~65 Development, debugging
4 289 ~250 Research, moderate resolution
5 1089 ~950 High-resolution studies
6 4225 ~3800 Production, detailed analysis

Core Methods

reset(initial_condition: Callable[[float, float], float]) -> None

Reset simulation with new initial condition.

Parameters:

  • initial_condition: Python function f(x, y) -> float that defines u(x,y,0)

Geometry Considerations:

  • Ensure initial condition is compatible with chosen geometry coordinate system
  • Use get_domain_bounds() to understand coordinate ranges
  • Consider symmetry and boundary conditions for optimal results

Example:

# Geometry-aware initial condition
bounds = env.get_domain_bounds()
center_x = (bounds[0] + bounds[1]) / 2
center_y = (bounds[2] + bounds[3]) / 2

def centered_gaussian(x, y):
    return np.exp(-10 * ((x - center_x)**2 + (y - center_y)**2))

env.reset(centered_gaussian)
step() -> None

Advance simulation by one time step. Solves the linear system for the next time level.

State Query Methods
get_time() -> float                    # Current simulation time
get_timestep() -> int                  # Current time step number
is_done() -> bool                      # True if t >= final_time
get_num_dofs() -> int                  # Number of degrees of freedom

Geometry Information Methods (New)

get_geometry_config() -> GeometryConfig

Returns the geometry configuration used to create the environment.

get_geometry_description() -> str

Returns a human-readable description of the current geometry.

Example Output:

  • "Hyper cube with side length 2.000000"
  • "Circle with radius 1.500000 centered at (0.000000, 0.000000)"
  • "Rectangle 3.000000 x 1.500000"
get_domain_bounds() -> List[float]

Returns the bounding box of the computational domain.

Returns: [x_min, x_max, y_min, y_max]

Use Cases:

  • Understanding coordinate system
  • Setting up appropriate initial conditions
  • Normalizing spatial coordinates for RL
bounds = env.get_domain_bounds()
x_range = bounds[1] - bounds[0]
y_range = bounds[3] - bounds[2]
domain_area_approx = x_range * y_range

Data Extraction Methods

get_solution_data() -> List[float]

Returns the complete finite element solution vector.

Returns: List of solution values, one per degree of freedom Use case: High-dimensional RL state representation, detailed analysis Memory: 8 bytes × number of DOF

solution = env.get_solution_data()
print(f"Solution dimension: {len(solution)}")
print(f"Value range: [{min(solution):.6f}, {max(solution):.6f}]")
get_mesh_points() -> List[List[float]]

Returns physical coordinates of all mesh points.

Returns: List of [x, y] coordinate pairs Use case: Spatial analysis, visualization, custom interpolation

Geometry Note: Coordinate ranges depend on chosen geometry. Use get_domain_bounds() for reference.

points = env.get_mesh_points()
x_coords = [p[0] for p in points]
y_coords = [p[1] for p in points]

# Verify coordinates are within expected bounds
bounds = env.get_domain_bounds()
assert all(bounds[0] <= x <= bounds[1] for x in x_coords)
assert all(bounds[2] <= y <= bounds[3] for y in y_coords)
get_physical_quantities() -> List[float]

Returns 5 physically meaningful summary quantities.

Returns: [total_energy, max_value, min_value, energy_center_x, energy_center_y]

Index Quantity Description Use Case
0 Total Energy ∫∫ u(x,y) dA Energy conservation analysis
1 Maximum Value max(u) Peak concentration tracking
2 Minimum Value min(u) Boundary condition verification
3 Energy Center X Weighted centroid x-coordinate Spatial dynamics
4 Energy Center Y Weighted centroid y-coordinate Spatial dynamics

Geometry Note: Energy center coordinates are in the coordinate system of the chosen geometry.

Use case: Compact RL state representation, reward function design

quantities = env.get_physical_quantities()
total_energy = quantities[0]
peak_value = quantities[1]
center_x, center_y = quantities[3], quantities[4]

# Check if energy center is within domain bounds
bounds = env.get_domain_bounds()
center_in_bounds = (bounds[0] <= center_x <= bounds[1] and 
                   bounds[2] <= center_y <= bounds[3])

Parameter Control Methods

set_diffusion_coefficient(K_new: float) -> None

Dynamically modify the diffusion coefficient during simulation.

Parameters:

  • K_new: New diffusion coefficient (must be > 0)

Note: This rebuilds internal matrices, so use sparingly within episodes.

# Start with slow diffusion
env = DiffusionEnvironment(
    diffusion_env.GeometryConfig.hyper_ball(radius=1.5),
    diffusion_coeff=0.05
)
env.reset(initial_condition)

# Speed up diffusion mid-simulation
for i in range(10):
    env.step()

env.set_diffusion_coefficient(0.2)  # Faster diffusion

for i in range(10):
    env.step()

Visualization Methods

write_vtk(filename: str) -> None

Write current solution to VTU file for visualization in VisIt or ParaView.

Parameters:

  • filename: Output filename (should end in .vtu)

Use case: Detailed visualization, publication figures, debugging

Geometry Note: VTU files correctly represent curved boundaries and complex geometries.

# Save solutions from different geometries
geometries = [
    ("square", diffusion_env.GeometryConfig.hyper_cube(1.0)),
    ("circle", diffusion_env.GeometryConfig.hyper_ball(0.56)),  # Same area
    ("l_shape", diffusion_env.GeometryConfig.l_shaped(1.0))
]

for name, config in geometries:
    env = diffusion_env.DiffusionEnvironment(config, refinement_level=4)
    env.reset(initial_condition)
    
    for i in range(20):
        env.step()
    
    env.write_vtk(f"solution_{name}_t{env.get_time():.3f}.vtu")

Convenience Functions

# Quick environment creation for common geometries
create_square_env(size: float = 1.0, refinement: int = 4) -> DiffusionEnvironment
create_circle_env(radius: float = 1.0, refinement: int = 4) -> DiffusionEnvironment  
create_rectangle_env(width: float = 1.0, height: float = 1.0, refinement: int = 4) -> DiffusionEnvironment

Example:

# These are equivalent
env1 = diffusion_env.DiffusionEnvironment(
    diffusion_env.GeometryConfig.hyper_cube(2.0), refinement_level=4
)

env2 = diffusion_env.create_square_env(size=2.0, refinement=4)

Usage Examples

Example 1: Basic Geometry Comparison

import diffusion_env
import numpy as np
import matplotlib.pyplot as plt

def compare_geometries():
    """Compare diffusion behavior in different geometries."""
    
    # Create environments with different geometries but similar "size"
    configs = [
        ("Square", diffusion_env.GeometryConfig.hyper_cube(1.0)),
        ("Circle", diffusion_env.GeometryConfig.hyper_ball(0.56)),  # Roughly same area
        ("Rectangle", diffusion_env.GeometryConfig.hyper_rectangle(1.2, 0.8)),  # Same area
    ]
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    for i, (name, config) in enumerate(configs):
        env = diffusion_env.DiffusionEnvironment(config, refinement_level=4)
        
        # Same initial condition for all: centered Gaussian
        bounds = env.get_domain_bounds()
        center_x = (bounds[0] + bounds[1]) / 2
        center_y = (bounds[2] + bounds[3]) / 2
        
        def centered_gaussian(x, y):
            return np.exp(-10 * ((x - center_x)**2 + (y - center_y)**2))
        
        env.reset(centered_gaussian)
        
        # Run for same number of steps
        for step in range(20):
            env.step()
        
        # Plot results
        points = np.array(env.get_mesh_points())
        solution = np.array(env.get_solution_data())
        
        scatter = axes[i].scatter(points[:, 0], points[:, 1], c=solution, 
                                cmap='hot', s=15, alpha=0.8)
        axes[i].set_aspect('equal')
        axes[i].set_title(f"{name}\nDoFs: {env.get_num_dofs()}")
        
        # Add geometry outline for circle
        if name == "Circle":
            circle = plt.Circle((center_x, center_y), 0.56, fill=False, color='black', linewidth=1)
            axes[i].add_patch(circle)
        
        print(f"{name}: {env.get_num_dofs()} DoFs, Max temp: {np.max(solution):.4f}")
    
    plt.tight_layout()
    plt.suptitle("Diffusion in Different Geometries (Same Area)", y=1.02)
    plt.show()

compare_geometries()

Example 2: Parameter Study Across Geometries

import diffusion_env
import numpy as np
import matplotlib.pyplot as plt

def geometry_parameter_study():
    """Study effect of diffusion coefficient across different geometries."""
    
    geometries = [
        ("Square", diffusion_env.GeometryConfig.hyper_cube(1.0)),
        ("Circle", diffusion_env.GeometryConfig.hyper_ball(0.56)),
        ("L-shaped", diffusion_env.GeometryConfig.l_shaped(1.0))
    ]
    
    diffusion_coeffs = [0.05, 0.1, 0.2]
    
    results = {}
    
    for geom_name, config in geometries:
        results[geom_name] = {}
        
        for K in diffusion_coeffs:
            env = diffusion_env.DiffusionEnvironment(
                config,
                diffusion_coeff=K, 
                dt=0.005, 
                final_time=0.5
            )
            
            # Geometry-aware initial condition
            bounds = env.get_domain_bounds()
            center_x = (bounds[0] + bounds[1]) / 2
            center_y = (bounds[2] + bounds[3]) / 2
            
            def standard_initial(x, y):
                return np.exp(-15 * ((x - center_x)**2 + (y - center_y)**2))
            
            env.reset(standard_initial)
            
            times, energies = [], []
            while not env.is_done():
                times.append(env.get_time())
                energies.append(env.get_physical_quantities()[0])
                env.step()
            
            results[geom_name][K] = {'times': times, 'energies': energies}
    
    # Plot results
    fig, axes = plt.subplots(1, len(geometries), figsize=(15, 5))
    
    for i, (geom_name, geom_results) in enumerate(results.items()):
        ax = axes[i]
        
        for K, data in geom_results.items():
            ax.plot(data['times'], data['energies'], 
                   label=f'K = {K}', linewidth=2)
        
        ax.set_xlabel('Time')
        ax.set_ylabel('Total Energy')
        ax.set_title(f'{geom_name} Geometry')
        ax.legend()
        ax.set_yscale('log')
        ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('geometry_parameter_study.png')
    plt.show()

geometry_parameter_study()

Example 3: Custom Initial Conditions for Different Geometries

import diffusion_env
import numpy as np

def create_geometry_specific_initial_conditions():
    """Create initial conditions tailored to each geometry type."""
    
    def square_initial(x, y):
        """Checkerboard pattern for square domain."""
        return 1.0 if (int(x*4) + int(y*4)) % 2 == 0 else 0.0
    
    def circle_initial(x, y):
        """Radial pattern for circular domain."""
        r = np.sqrt(x**2 + y**2)
        return np.cos(4 * np.pi * r) * np.exp(-2 * r**2)
    
    def rectangle_initial(x, y):
        """Wave pattern along length for rectangular domain."""
        return np.sin(2 * np.pi * x) * np.exp(-5 * (y - 0.75)**2)
    
    def l_shape_initial(x, y):
        """Multiple sources for L-shaped domain."""
        source1 = np.exp(-20 * ((x - 0.2)**2 + (y - 0.2)**2))
        source2 = np.exp(-20 * ((x - 0.8)**2 + (y - 0.5)**2))
        return source1 + 0.5 * source2
    
    def annulus_initial(x, y):
        """Ring pattern for annular domain."""
        r = np.sqrt(x**2 + y**2)
        target_r = 0.65  # Middle of annulus
        return np.exp(-50 * (r - target_r)**2)
    
    # Test each geometry with its custom initial condition
    test_cases = [
        ("Square", diffusion_env.GeometryConfig.hyper_cube(1.0), square_initial),
        ("Circle", diffusion_env.GeometryConfig.hyper_ball(0.8), circle_initial),
        ("Rectangle", diffusion_env.GeometryConfig.hyper_rectangle(2.0, 1.5), rectangle_initial),
        ("L-shaped", diffusion_env.GeometryConfig.l_shaped(1.0), l_shape_initial),
        ("Annulus", diffusion_env.GeometryConfig.hyper_shell(0.4, 0.9), annulus_initial)
    ]
    
    for name, config, initial_func in test_cases:
        env = diffusion_env.DiffusionEnvironment(config, refinement_level=4)
        env.reset(initial_func)
        
        print(f"{name} Environment:")
        print(f"  Description: {env.get_geometry_description()}")
        print(f"  DOFs: {env.get_num_dofs()}")
        print(f"  Domain bounds: {env.get_domain_bounds()}")
        
        # Run a few steps and check energy
        initial_energy = env.get_physical_quantities()[0]
        for _ in range(10):
            env.step()
        final_energy = env.get_physical_quantities()[0]
        
        print(f"  Initial energy: {initial_energy:.6f}")
        print(f"  Energy after 10 steps: {final_energy:.6f}")
        print(f"  Energy ratio: {final_energy/initial_energy:.6f}")
        print()
        
        # Optionally save for visualization
        env.write_vtk(f"initial_condition_{name.lower()}.vtu")

create_geometry_specific_initial_conditions()

Example 4: Adaptive Initial Conditions

import diffusion_env
import numpy as np

def create_adaptive_initial_condition(env):
    """Create initial condition that automatically adapts to any geometry."""
    
    bounds = env.get_domain_bounds()
    geometry_config = env.get_geometry_config()
    
    center_x = (bounds[0] + bounds[1]) / 2
    center_y = (bounds[2] + bounds[3]) / 2
    domain_scale = max(bounds[1] - bounds[0], bounds[3] - bounds[2])
    
    def adaptive_initial(x, y):
        # Adjust pattern based on geometry type
        if geometry_config.type == diffusion_env.GeometryType.HYPER_BALL:
            # For circles, use radial symmetry
            r = np.sqrt((x - center_x)**2 + (y - center_y)**2)
            max_radius = geometry_config.radius * 0.8
            return np.exp(-5 * (r / max_radius)**2) if r <= max_radius else 0.0
            
        elif geometry_config.type == diffusion_env.GeometryType.HYPER_SHELL:
            # For annulus, create ring pattern
            r = np.sqrt((x - center_x)**2 + (y - center_y)**2)
            target_radius = (geometry_config.inner_radius + geometry_config.outer_radius) / 2
            return np.exp(-20 * (r - target_radius)**2)
            
        elif geometry_config.type == diffusion_env.GeometryType.HYPER_RECTANGLE:
            # For rectangles, use aspect ratio aware pattern
            norm_x = (x - bounds[0]) / (bounds[1] - bounds[0])
            norm_y = (y - bounds[2]) / (bounds[3] - bounds[2])
            return np.sin(2 * np.pi * norm_x) * np.sin(np.pi * norm_y)
            
        elif geometry_config.type == diffusion_env.GeometryType.L_SHAPED:
            # For L-shape, multiple sources in different arms
            source1 = np.exp(-15/domain_scale * ((x - center_x*0.5)**2 + (y - center_y*0.5)**2))
            source2 = np.exp(-15/domain_scale * ((x - center_x*1.2)**2 + (y - center_y*1.2)**2))
            return source1 + 0.6 * source2
            
        else:
            # Default: centered Gaussian scaled to domain
            width = 10.0 / domain_scale
            return np.exp(-width * ((x - center_x)**2 + (y - center_y)**2))
    
    return adaptive_initial

# Test adaptive initial conditions
geometries = [
    diffusion_env.GeometryConfig.hyper_cube(1.5),
    diffusion_env.GeometryConfig.hyper_ball(1.0),
    diffusion_env.GeometryConfig.hyper_rectangle(2.0, 1.2),
    diffusion_env.GeometryConfig.hyper_shell(0.4, 1.0),
    diffusion_env.GeometryConfig.l_shaped(1.5)
]

print("Testing adaptive initial conditions:")
print("=" * 50)

for config in geometries:
    env = diffusion_env.DiffusionEnvironment(config, refinement_level=4)
    
    # Create and apply adaptive initial condition
    initial_condition = create_adaptive_initial_condition(env)
    env.reset(initial_condition)
    
    # Get initial state information
    quantities = env.get_physical_quantities()
    
    print(f"Geometry: {env.get_geometry_description()}")
    print(f"  Domain bounds: {env.get_domain_bounds()}")
    print(f"  Initial energy: {quantities[0]:.6f}")
    print(f"  Max value: {quantities[1]:.6f}")
    print(f"  Energy center: ({quantities[3]:.3f}, {quantities[4]:.3f})")
    
    # Run simulation and track energy evolution
    energy_history = [quantities[0]]
    for _ in range(15):
        env.step()
        energy_history.append(env.get_physical_quantities()[0])
    
    # Calculate energy decay rate
    initial_energy = energy_history[0]
    final_energy = energy_history[-1]
    decay_ratio = final_energy / initial_energy
    
    print(f"  Energy decay ratio: {decay_ratio:.6f}")
    print()

Performance Characteristics

Computational Complexity

Operation Time Complexity Memory Notes
Environment Creation O(N log N) O(N) One-time cost, N = DOF
Reset O(N) O(1) Per episode
Step O(N^1.5) O(1) Sparse solver
Data Extraction O(N) O(N) Copy overhead

Scaling Behavior by Geometry

# Typical performance on modern workstation (Intel i7, 16GB RAM)
# Times in milliseconds

Geometry Type       | Refinement | DOF  | Creation | Reset | Step | Extraction
Square             | 4          | 289  | 45 ms    | 2 ms  | 8 ms | 0.3 ms
Circle             | 4          | ~250 | 50 ms    | 2 ms  | 7 ms | 0.3 ms
Rectangle (2:1)    | 4          | 289  | 45 ms    | 2 ms  | 8 ms | 0.3 ms
L-shaped           | 4          | ~200 | 55 ms    | 2 ms  | 6 ms | 0.2 ms
Annulus            | 4          | ~220 | 60 ms    | 2 ms  | 7 ms | 0.3 ms
Quarter Circle     | 4          | ~120 | 40 ms    | 1 ms  | 4 ms | 0.1 ms

Memory Usage by Geometry

  • Core Environment: ~50 MB + 200 bytes × DOF
  • Per State: 24 bytes × DOF (solution + coordinates)
  • Temporary Storage: ~100 bytes × DOF during stepping
  • Geometry Overhead: <1 MB additional for curved boundary information

Optimization Tips

  1. Choose Appropriate Geometry and Mesh Size:

    • Start with quarter circles or small rectangles for rapid prototyping
    • Use refinement level 4 for development, 5+ for production
  2. Geometry-Specific Considerations:

    • Circular domains may be slightly more expensive due to curved boundaries
    • L-shaped domains have fewer DOFs but more complex assembly
    • Rectangular domains are most efficient for given refinement level
  3. Reuse Environments: Create once, reset many times for different episodes

  4. Batch Data Extraction: Extract data only when needed for RL updates

  5. Monitor Memory: Large replay buffers can consume significant memory with high-dimensional states

  6. Geometry Selection Strategy: Use simpler geometries during early training phases

Advanced Usage

Custom Geometry-Aware Initial Conditions

import diffusion_env
import numpy as np

def create_physics_informed_initial_condition(env, heat_sources=None):
    """
    Create physically meaningful initial conditions based on geometry.
    
    Args:
        env: DiffusionEnvironment instance
        heat_sources: List of (x, y, strength) tuples for heat sources
    """
    bounds = env.get_domain_bounds()
    config = env.get_geometry_config()
    
    if heat_sources is None:
        # Default: single centered source
        center_x = (bounds[0] + bounds[1]) / 2
        center_y = (bounds[2] + bounds[3]) / 2
        heat_sources = [(center_x, center_y, 1.0)]
    
    def physics_initial(x, y):
        total = 0.0
        
        for source_x, source_y, strength in heat_sources:
            # Distance from source
            distance = np.sqrt((x - source_x)**2 + (y - source_y)**2)
            
            # Geometry-aware width calculation
            if config.type == diffusion_env.GeometryType.HYPER_BALL:
                # Scale with radius
                width = 10.0 * config.radius
                
            elif config.type == diffusion_env.GeometryType.HYPER_SHELL:
                # Scale with shell thickness
                thickness = config.outer_radius - config.inner_radius
                width = 20.0 / thickness
                
            elif config.type in [diffusion_env.GeometryType.HYPER_CUBE, 
                               diffusion_env.GeometryType.HYPER_RECTANGLE]:
                # Scale with domain size
                domain_scale = max(bounds[1] - bounds[0], bounds[3] - bounds[2])
                width = 15.0 / domain_scale
                
            else:
                # Default scaling
                width = 15.0
            
            # Add source contribution
            contribution = strength * np.exp(-width * distance**2)
            total += contribution
        
        return total
    
    return physics_initial

# Example usage
geometries = [
    diffusion_env.GeometryConfig.hyper_cube(2.0),
    diffusion_env.GeometryConfig.hyper_ball(1.2),
    diffusion_env.GeometryConfig.hyper_shell(0.5, 1.5),
    diffusion_env.GeometryConfig.hyper_rectangle(3.0, 1.5)
]

for config in geometries:
    env = diffusion_env.DiffusionEnvironment(config, refinement_level=4)
    
    print(f"Geometry: {env.get_geometry_description()}")
    
    # Single source case
    initial_single = create_physics_informed_initial_condition(env)
    env.reset(initial_single)
    print(f"  Single source energy: {env.get_physical_quantities()[0]:.6f}")
    
    # Multiple sources case
    bounds = env.get_domain_bounds()
    multi_sources = [
        (bounds[0] + 0.3 * (bounds[1] - bounds[0]), 
         bounds[2] + 0.3 * (bounds[3] - bounds[2]), 1.0),
        (bounds[0] + 0.7 * (bounds[1] - bounds[0]), 
         bounds[2] + 0.7 * (bounds[3] - bounds[2]), 0.8)
    ]
    
    initial_multi = create_physics_informed_initial_condition(env, multi_sources)
    env.reset(initial_multi)
    print(f"  Multi source energy: {env.get_physical_quantities()[0]:.6f}")
    print()

Batch Processing with Different Geometries

import diffusion_env
import numpy as np
from concurrent.futures import ThreadPoolExecutor
from typing import List, Tuple, Dict

def run_geometry_simulation(params: Tuple) -> Dict:
    """Run simulation with specific geometry and parameters."""
    config, sim_params, initial_params = params
    
    env = diffusion_env.DiffusionEnvironment(config, **sim_params)
    
    # Geometry-adaptive initial condition
    bounds = env.get_domain_bounds()
    center_x, center_y, width_factor = initial_params
    
    # Normalize to domain
    actual_center_x = bounds[0] + center_x * (bounds[1] - bounds[0])
    actual_center_y = bounds[2] + center_y * (bounds[3] - bounds[2])
    domain_scale = max(bounds[1] - bounds[0], bounds[3] - bounds[2])
    actual_width = width_factor / domain_scale
    
    def initial_condition(x, y):
        return np.exp(-actual_width * ((x - actual_center_x)**2 + (y - actual_center_y)**2))
    
    env.reset(initial_condition)
    
    # Run simulation and collect data
    time_series = []
    while not env.is_done():
        quantities = env.get_physical_quantities()
        time_series.append({
            'time': env.get_time(),
            'energy': quantities[0],
            'max_value': quantities[1],
            'center_x': quantities[3],
            'center_y': quantities[4]
        })
        env.step()
    
    return {
        'geometry_type': config.type.name,
        'description': env.get_geometry_description(),
        'num_dofs': env.get_num_dofs(),
        'domain_bounds': env.get_domain_bounds(),
        'time_series': time_series,
        'final_energy': time_series[-1]['energy'] if time_series else 0.0
    }

def batch_geometry_comparison():
    """Run batch simulations across multiple geometries and parameters."""
    
    # Define geometries to test
    geometries = [
        diffusion_env.GeometryConfig.hyper_cube(1.0),
        diffusion_env.GeometryConfig.hyper_ball(0.56),  # Same area as unit square
        diffusion_env.GeometryConfig.hyper_rectangle(1.3, 0.77),  # Same area
        diffusion_env.GeometryConfig.l_shaped(1.0),
        diffusion_env.GeometryConfig.hyper_shell(0.3, 0.65)  # Same area
    ]
    
    # Simulation parameters
    sim_params = {
        'refinement_level': 4,
        'diffusion_coeff': 0.1,
        'dt': 0.005,
        'final_time': 0.4
    }
    
    # Initial condition parameters (normalized coordinates)
    initial_conditions = [
        [0.3, 0.3, 15.0],  # Lower-left
        [0.7, 0.7, 15.0],  # Upper-right
        [0.5, 0.5, 10.0],  # Centered, wide
        [0.5, 0.5, 25.0],  # Centered, narrow
    ]
    
    # Create parameter combinations
    params_list = []
    for geometry in geometries:
        for initial_params in initial_conditions:
            params_list.append((geometry, sim_params, initial_params))
    
    print(f"Running {len(params_list)} simulations...")
    
    # Run simulations in parallel
    with ThreadPoolExecutor(max_workers=4) as executor:
        results = list(executor.map(run_geometry_simulation, params_list))
    
    # Analyze results by geometry type
    geometry_stats = {}
    for result in results:
        geom_type = result['geometry_type']
        if geom_type not in geometry_stats:
            geometry_stats[geom_type] = {
                'count': 0,
                'dof_counts': [],
                'final_energies': [],
                'energy_decay_rates': []
            }
        
        stats = geometry_stats[geom_type]
        stats['count'] += 1
        stats['dof_counts'].append(result['num_dofs'])
        stats['final_energies'].append(result['final_energy'])
        
        # Calculate energy decay rate
        time_series = result['time_series']
        if len(time_series) > 1:
            initial_energy = time_series[0]['energy']
            final_energy = time_series[-1]['energy']
            total_time = time_series[-1]['time']
            
            if initial_energy > 0 and total_time > 0:
                decay_rate = -np.log(final_energy / initial_energy) / total_time
                stats['energy_decay_rates'].append(decay_rate)
    
    # Print summary
    print("\nGeometry Performance Summary:")
    print("=" * 70)
    print(f"{'Geometry':<15} {'Count':<6} {'Avg DOF':<8} {'Avg Final Energy':<16} {'Avg Decay Rate':<15}")
    print("-" * 70)
    
    for geom_type, stats in geometry_stats.items():
        avg_dof = np.mean(stats['dof_counts'])
        avg_energy = np.mean(stats['final_energies'])
        avg_decay = np.mean(stats['energy_decay_rates']) if stats['energy_decay_rates'] else 0.0
        
        print(f"{geom_type:<15} {stats['count']:<6} {avg_dof:<8.0f} "
              f"{avg_energy:<16.6f} {avg_decay:<15.6f}")
    
    return results, geometry_stats

# Run batch comparison
if __name__ == "__main__":
    results, stats = batch_geometry_comparison()
    
    # Additional analysis
    print("\nDetailed Analysis:")
    print("-" * 30)
    
    for geom_type, stat in stats.items():
        print(f"\n{geom_type}:")
        print(f"  DOF range: {min(stat['dof_counts'])} - {max(stat['dof_counts'])}")
        print(f"  Energy std: {np.std(stat['final_energies']):.6f}")
        if stat['energy_decay_rates']:
            print(f"  Decay rate std: {np.std(stat['energy_decay_rates']):.6f}")

Troubleshooting

Common Installation Issues

deal.II Not Found

Error: Could not find deal.II-config tool

Solutions:

  1. Verify installation: deal.II-config --version
  2. Set environment: export DEAL_II_DIR=/path/to/dealii
  3. Add to PATH: export PATH=$DEAL_II_DIR/bin:$PATH

Compilation Errors

CMake Error: target_link_libraries keyword signature conflict

Solution: Use the corrected CMakeLists.txt provided with this package that handles pybind11/deal.II compatibility.

Python Import Errors

ImportError: dynamic module does not define module export function

Solutions:

  1. Verify Python environment: which python
  2. Rebuild extension: python setup.py build_ext --inplace
  3. Check shared library dependencies: ldd diffusion_env*.so

Geometry-Related Issues

Invalid Geometry Configuration

std::invalid_argument: Shell radii must satisfy 0 < inner_radius < outer_radius

Solutions:

  1. Verify geometry parameters are physically valid
  2. Check that radii are positive for circular geometries
  3. Ensure outer radius > inner radius for shell geometries
  4. Verify rectangle dimensions are positive
# Correct geometry configurations
config1 = diffusion_env.GeometryConfig.hyper_shell(0.5, 1.2)  # Valid: inner < outer
config2 = diffusion_env.GeometryConfig.hyper_rectangle(2.0, 1.5)  # Valid: positive dimensions

# Invalid configurations (will raise errors)
# config_bad1 = diffusion_env.GeometryConfig.hyper_shell(1.2, 0.5)  # inner > outer
# config_bad2 = diffusion_env.GeometryConfig.hyper_ball(-1.0)  # negative radius

Initial Condition Out of Bounds

Warning: Initial condition evaluated outside domain bounds

Solutions:

  1. Use get_domain_bounds() to understand coordinate system
  2. Create geometry-aware initial conditions
  3. Test initial conditions with visualization
# Geometry-aware initial condition
def create_safe_initial_condition(env):
    bounds = env.get_domain_bounds()
    center_x = (bounds[0] + bounds[1]) / 2
    center_y = (bounds[2] + bounds[3]) / 2
    
    def safe_initial(x, y):
        # Ensure we're working within bounds
        if not (bounds[0] <= x <= bounds[1] and bounds[2] <= y <= bounds[3]):
            return 0.0
        return np.exp(-10 * ((x - center_x)**2 + (y - center_y)**2))
    
    return safe_initial

Runtime Issues

Memory Errors

std::bad_alloc: Cannot allocate memory

Solutions:

  1. Reduce refinement level
  2. Use simpler geometries during development
  3. Increase system memory
  4. Check for memory leaks in long-running scripts

Numerical Issues

Solution values become NaN or extremely large

Solutions:

  1. Reduce time step size
  2. Check initial condition validity for chosen geometry
  3. Verify diffusion coefficient is positive
  4. Ensure initial condition is compatible with boundary conditions
# Debugging numerical issues
env = diffusion_env.DiffusionEnvironment(config, dt=0.001)  # Smaller time step

def debug_initial(x, y):
    value = your_initial_condition(x, y)
    if not np.isfinite(value) or value < 0:
        print(f"Warning: Invalid initial value {value} at ({x}, {y})")
        return 0.0
    return value

env.reset(debug_initial)

Debugging Tips

Enable Geometry Debugging

import diffusion_env
import numpy as np

def debug_geometry_setup(config):
    """Debug geometry configuration and setup."""
    print(f"Geometry config: {config.type.name}")
    
    # Create environment and check setup
    env = diffusion_env.DiffusionEnvironment(config, refinement_level=3)
    
    print(f"Description: {env.get_geometry_description()}")
    print(f"DOFs: {env.get_num_dofs()}")
    print(f"Bounds: {env.get_domain_bounds()}")
    
    # Test with simple initial condition
    bounds = env.get_domain_bounds()
    center_x = (bounds[0] + bounds[1]) / 2
    center_y = (bounds[2] + bounds[3]) / 2
    
    def test_initial(x, y):
        return 1.0 if abs(x - center_x) < 0.1 and abs(y - center_y) < 0.1 else 0.0
    
    try:
        env.reset(test_initial)
        print("✓ Initial condition set successfully")
        
        env.step()
        quantities = env.get_physical_quantities()
        print(f"✓ First step completed, energy: {quantities[0]:.6f}")
        
    except Exception as e:
        print(f"✗ Error during simulation: {e}")
    
    return env

# Test all geometry types
geometry_configs = [
    diffusion_env.GeometryConfig.hyper_cube(1.0),
    diffusion_env.GeometryConfig.hyper_ball(0.8),
    diffusion_env.GeometryConfig.hyper_rectangle(1.5, 1.0),
    diffusion_env.GeometryConfig.hyper_shell(0.3, 0.8),
    diffusion_env.GeometryConfig.l_shaped(1.0),
    diffusion_env.GeometryConfig.quarter_hyper_ball(1.0)
]

for config in geometry_configs:
    print("=" * 50)
    debug_geometry_setup(config)

Contributing

We welcome contributions to DiffusionEnv! Here's how to get involved:

Development Setup

# Fork the repository and clone your fork
git clone https://github.com/yourusername/diffusion-env.git
cd diffusion-env

# Create development environment
python -m venv dev_env
source dev_env/bin/activate

# Install development dependencies
pip install -e .[dev]
pip install pytest black flake8 mypy

Code Style

We follow standard Python conventions:

  • Formatting: Use black for code formatting
  • Linting: Use flake8 for linting
  • Type Hints: Use mypy for type checking
  • Documentation: Use Google-style docstrings
# Run code quality checks
black diffusion_env/
flake8 diffusion_env/
mypy diffusion_env/

Testing

All contributions must include tests:

# Run test suite
pytest tests/ -v

# Run with coverage
pytest tests/ --cov=diffusion_env --cov-report=html

Submitting Changes

  1. Create a feature branch: git checkout -b feature-name
  2. Make your changes with tests
  3. Run quality checks: black, flake8, mypy, pytest
  4. Commit with clear message: git commit -m "Add feature X"
  5. Push branch: git push origin feature-name
  6. Open pull request with description of changes

Areas for Contribution

  • New Geometry Types: Additional shapes like ellipses, polygons, or custom domains
  • Enhanced Physics: Different PDEs, boundary conditions, or material properties
  • Performance: Optimization of critical loops, memory usage improvements
  • Visualization: Enhanced plotting utilities, interactive visualizations
  • Documentation: Examples, tutorials, API improvements
  • Testing: Extended test coverage, performance benchmarks

Citation

If you use DiffusionEnv in your research, please cite:

@software{diffusion_env,
  title={DiffusionEnv: High-Performance Finite Element Environment for Reinforcement Learning},
  author={Your Name},
  year={2024},
  url={https://github.com/yourusername/diffusion-env},
  note={Python package for physics-informed reinforcement learning with multiple geometry support}
}

License

This project is licensed under the MIT License - see the LICENSE file for details.

Acknowledgments

  • deal.II development team for the excellent finite element library
  • pybind11 developers for seamless Python-C++ integration
  • The computational science community for inspiration and feedback
  • Contributors and users who help improve this package

Related Work

  • deal.II: The finite element library powering this environment
  • FEniCS: Alternative FE library with Python bindings
  • OpenAI Gym: Standard RL environment interface
  • Stable-Baselines3: Popular RL algorithms implementation
  • Physics-Informed Neural Networks: Related ML approach for PDEs

For questions, bug reports, or feature requests, please open an issue on GitHub or contact the maintainers.

About

Auto regressive model on Diffusion

Resources

Stars

0 stars

Watchers

0 watching

Forks

Releases

Packages

Used by

Contributors

Languages