CLARISSA Tutorial 05: Constraint Engine¶

Learning Objectives:

• Use Z3 solver for constraint satisfaction
• Implement physics constraints for reservoir simulation
• Validate and convert units automatically
• Apply constraint propagation to deck generation

Prerequisites: Notebooks 01-04

Estimated Time: 45 minutes

Why Constraints Matter¶

LLMs can generate plausible-looking simulation decks that are physically impossible:

Bad Generation	Why It's Wrong
Sw + So = 1.3	Saturations must sum to 1.0
P_res = -500 psi	Negative pressure impossible
k_v > k_h	Unusual (though not impossible)
BHP > P_res	Well can't flow

The Constraint Engine catches these errors before simulation.

In [ ]:

# Setup
try:
    from z3 import *
    Z3_AVAILABLE = True
    print('Z3 solver available')
except ImportError:
    Z3_AVAILABLE = False
    print('Z3 not installed - using fallback validation')
    print('Install with: pip install z3-solver')

from dataclasses import dataclass, field
from typing import List, Dict, Optional, Tuple, Any
from enum import Enum, auto
import math

# Setup try: from z3 import * Z3_AVAILABLE = True print('Z3 solver available') except ImportError: Z3_AVAILABLE = False print('Z3 not installed - using fallback validation') print('Install with: pip install z3-solver') from dataclasses import dataclass, field from typing import List, Dict, Optional, Tuple, Any from enum import Enum, auto import math

Section 1: Unit System and Conversion¶

ECLIPSE supports FIELD and METRIC units. We need automatic conversion and validation.

In [ ]:

class UnitSystem(Enum):
    FIELD = auto()   # ft, psi, stb, mscf, md, F
    METRIC = auto()  # m, bar, sm3, md, C

@dataclass
class UnitConversion:
    """Unit conversion factors (FIELD -> METRIC)."""
    length: float = 0.3048        # ft -> m
    pressure: float = 0.0689476   # psi -> bar
    volume_oil: float = 0.158987  # stb -> sm3
    volume_gas: float = 28.3168   # mscf -> sm3
    temperature: float = 5/9      # F offset -> C offset
    temperature_offset: float = -32 * 5/9 + 273.15  # F -> K
    permeability: float = 1.0     # md is same in both

CONVERSIONS = UnitConversion()

def convert_value(value: float, from_unit: str, to_system: UnitSystem) -> float:
    """Convert a value between unit systems."""
    if to_system == UnitSystem.FIELD:
        # METRIC -> FIELD (inverse)
        factors = {
            'length': 1/CONVERSIONS.length,
            'pressure': 1/CONVERSIONS.pressure,
            'volume_oil': 1/CONVERSIONS.volume_oil,
            'permeability': 1.0
        }
    else:
        # FIELD -> METRIC
        factors = {
            'length': CONVERSIONS.length,
            'pressure': CONVERSIONS.pressure,
            'volume_oil': CONVERSIONS.volume_oil,
            'permeability': 1.0
        }
    
    return value * factors.get(from_unit, 1.0)

# Demo conversion
depth_ft = 8500
depth_m = convert_value(depth_ft, 'length', UnitSystem.METRIC)
print(f"Depth: {depth_ft} ft = {depth_m:.1f} m")

pressure_psi = 3800
pressure_bar = convert_value(pressure_psi, 'pressure', UnitSystem.METRIC)
print(f"Pressure: {pressure_psi} psi = {pressure_bar:.1f} bar")

class UnitSystem(Enum): FIELD = auto() # ft, psi, stb, mscf, md, F METRIC = auto() # m, bar, sm3, md, C @dataclass class UnitConversion: """Unit conversion factors (FIELD -> METRIC).""" length: float = 0.3048 # ft -> m pressure: float = 0.0689476 # psi -> bar volume_oil: float = 0.158987 # stb -> sm3 volume_gas: float = 28.3168 # mscf -> sm3 temperature: float = 5/9 # F offset -> C offset temperature_offset: float = -32 * 5/9 + 273.15 # F -> K permeability: float = 1.0 # md is same in both CONVERSIONS = UnitConversion() def convert_value(value: float, from_unit: str, to_system: UnitSystem) -> float: """Convert a value between unit systems.""" if to_system == UnitSystem.FIELD: # METRIC -> FIELD (inverse) factors = { 'length': 1/CONVERSIONS.length, 'pressure': 1/CONVERSIONS.pressure, 'volume_oil': 1/CONVERSIONS.volume_oil, 'permeability': 1.0 } else: # FIELD -> METRIC factors = { 'length': CONVERSIONS.length, 'pressure': CONVERSIONS.pressure, 'volume_oil': CONVERSIONS.volume_oil, 'permeability': 1.0 } return value * factors.get(from_unit, 1.0) # Demo conversion depth_ft = 8500 depth_m = convert_value(depth_ft, 'length', UnitSystem.METRIC) print(f"Depth: {depth_ft} ft = {depth_m:.1f} m") pressure_psi = 3800 pressure_bar = convert_value(pressure_psi, 'pressure', UnitSystem.METRIC) print(f"Pressure: {pressure_psi} psi = {pressure_bar:.1f} bar")

Section 2: Physics Constraints with Z3¶

Z3 is a theorem prover that can verify constraint satisfaction.

In [ ]:

class ConstraintResult:
    """Result of constraint checking."""
    def __init__(self, satisfied: bool, message: str, fixes: List[str] = None):
        self.satisfied = satisfied
        self.message = message
        self.fixes = fixes or []
    
    def __repr__(self):
        status = 'PASS' if self.satisfied else 'FAIL'
        return f"[{status}] {self.message}"

class PhysicsConstraints:
    """Physics-based constraints for reservoir simulation."""
    
    @staticmethod
    def check_saturation(sw: float, so: float, sg: float = 0.0) -> ConstraintResult:
        """Saturations must sum to 1.0."""
        total = sw + so + sg
        if abs(total - 1.0) < 1e-6:
            return ConstraintResult(True, f"Saturations sum to {total:.4f}")
        else:
            # Suggest normalized values
            fixes = [
                f"Normalize: Sw={sw/total:.4f}, So={so/total:.4f}, Sg={sg/total:.4f}"
            ]
            return ConstraintResult(False, f"Saturations sum to {total:.4f}, not 1.0", fixes)
    
    @staticmethod
    def check_pressure_positive(pressure: float, name: str = "Pressure") -> ConstraintResult:
        """Pressure must be positive."""
        if pressure > 0:
            return ConstraintResult(True, f"{name} = {pressure:.1f} (positive)")
        else:
            return ConstraintResult(False, f"{name} = {pressure:.1f} is non-positive",
                                   [f"Set {name} to a positive value (e.g., 14.7 psi minimum)"])
    
    @staticmethod
    def check_pressure_gradient(depth: float, pressure: float, 
                                fluid: str = 'oil') -> ConstraintResult:
        """Check if pressure is reasonable for depth."""
        # Typical gradients (psi/ft)
        gradients = {'oil': 0.35, 'water': 0.45, 'gas': 0.1}
        expected_min = depth * gradients.get(fluid, 0.3)
        expected_max = depth * 0.6  # Overpressured
        
        if expected_min <= pressure <= expected_max:
            return ConstraintResult(True, 
                f"Pressure {pressure:.0f} psi reasonable for {depth:.0f} ft")
        else:
            return ConstraintResult(False,
                f"Pressure {pressure:.0f} psi unusual for {depth:.0f} ft depth",
                [f"Expected range: {expected_min:.0f} - {expected_max:.0f} psi"])
    
    @staticmethod
    def check_well_can_flow(p_reservoir: float, bhp: float, 
                            is_injector: bool = False) -> ConstraintResult:
        """Verify well can flow (producer: BHP < Pres, injector: BHP > Pres)."""
        if is_injector:
            if bhp > p_reservoir:
                return ConstraintResult(True, 
                    f"Injector BHP ({bhp:.0f}) > Pres ({p_reservoir:.0f}): Can inject")
            else:
                return ConstraintResult(False,
                    f"Injector BHP ({bhp:.0f}) <= Pres ({p_reservoir:.0f}): Cannot inject",
                    [f"Increase BHP above {p_reservoir + 100:.0f} psi"])
        else:
            if bhp < p_reservoir:
                return ConstraintResult(True,
                    f"Producer BHP ({bhp:.0f}) < Pres ({p_reservoir:.0f}): Can produce")
            else:
                return ConstraintResult(False,
                    f"Producer BHP ({bhp:.0f}) >= Pres ({p_reservoir:.0f}): Cannot flow",
                    [f"Decrease BHP below {p_reservoir - 100:.0f} psi"])
    
    @staticmethod
    def check_porosity(porosity: float) -> ConstraintResult:
        """Porosity must be between 0 and 1."""
        if 0 < porosity < 1:
            if porosity > 0.4:
                return ConstraintResult(True, 
                    f"Porosity {porosity:.2%} valid but unusually high",
                    ["Consider if this is realistic for your reservoir"])
            return ConstraintResult(True, f"Porosity {porosity:.2%} valid")
        else:
            return ConstraintResult(False, f"Porosity {porosity} out of range [0, 1]",
                                   ["Typical values: 0.05 - 0.35"])
    
    @staticmethod
    def check_permeability_anisotropy(kh: float, kv: float) -> ConstraintResult:
        """Check permeability anisotropy ratio."""
        if kv <= 0 or kh <= 0:
            return ConstraintResult(False, "Permeability must be positive")
        
        ratio = kv / kh
        if ratio > 1:
            return ConstraintResult(False, 
                f"kv/kh = {ratio:.2f} > 1 (kv > kh is unusual)",
                ["Typically kv < kh due to layering", "Verify this is intentional"])
        elif ratio < 0.01:
            return ConstraintResult(True,
                f"kv/kh = {ratio:.3f} (very low - tight barriers)",
                ["This will significantly limit vertical flow"])
        else:
            return ConstraintResult(True, f"kv/kh = {ratio:.2f} (reasonable)")

# Demo constraint checking
constraints = PhysicsConstraints()

print("Physics Constraint Checks:")
print("=" * 50)

# Good values
print(constraints.check_saturation(0.3, 0.7, 0.0))
print(constraints.check_pressure_gradient(8500, 3800))
print(constraints.check_porosity(0.22))

print("\nProblematic values:")
# Bad values
print(constraints.check_saturation(0.5, 0.6, 0.2))  # Sum > 1
print(constraints.check_well_can_flow(3800, 4000, is_injector=False))  # BHP > Pres
print(constraints.check_permeability_anisotropy(100, 150))  # kv > kh

class ConstraintResult: """Result of constraint checking.""" def __init__(self, satisfied: bool, message: str, fixes: List[str] = None): self.satisfied = satisfied self.message = message self.fixes = fixes or [] def __repr__(self): status = 'PASS' if self.satisfied else 'FAIL' return f"[{status}] {self.message}" class PhysicsConstraints: """Physics-based constraints for reservoir simulation.""" @staticmethod def check_saturation(sw: float, so: float, sg: float = 0.0) -> ConstraintResult: """Saturations must sum to 1.0.""" total = sw + so + sg if abs(total - 1.0) < 1e-6: return ConstraintResult(True, f"Saturations sum to {total:.4f}") else: # Suggest normalized values fixes = [ f"Normalize: Sw={sw/total:.4f}, So={so/total:.4f}, Sg={sg/total:.4f}" ] return ConstraintResult(False, f"Saturations sum to {total:.4f}, not 1.0", fixes) @staticmethod def check_pressure_positive(pressure: float, name: str = "Pressure") -> ConstraintResult: """Pressure must be positive.""" if pressure > 0: return ConstraintResult(True, f"{name} = {pressure:.1f} (positive)") else: return ConstraintResult(False, f"{name} = {pressure:.1f} is non-positive", [f"Set {name} to a positive value (e.g., 14.7 psi minimum)"]) @staticmethod def check_pressure_gradient(depth: float, pressure: float, fluid: str = 'oil') -> ConstraintResult: """Check if pressure is reasonable for depth.""" # Typical gradients (psi/ft) gradients = {'oil': 0.35, 'water': 0.45, 'gas': 0.1} expected_min = depth * gradients.get(fluid, 0.3) expected_max = depth * 0.6 # Overpressured if expected_min <= pressure <= expected_max: return ConstraintResult(True, f"Pressure {pressure:.0f} psi reasonable for {depth:.0f} ft") else: return ConstraintResult(False, f"Pressure {pressure:.0f} psi unusual for {depth:.0f} ft depth", [f"Expected range: {expected_min:.0f} - {expected_max:.0f} psi"]) @staticmethod def check_well_can_flow(p_reservoir: float, bhp: float, is_injector: bool = False) -> ConstraintResult: """Verify well can flow (producer: BHP < Pres, injector: BHP > Pres).""" if is_injector: if bhp > p_reservoir: return ConstraintResult(True, f"Injector BHP ({bhp:.0f}) > Pres ({p_reservoir:.0f}): Can inject") else: return ConstraintResult(False, f"Injector BHP ({bhp:.0f}) <= Pres ({p_reservoir:.0f}): Cannot inject", [f"Increase BHP above {p_reservoir + 100:.0f} psi"]) else: if bhp < p_reservoir: return ConstraintResult(True, f"Producer BHP ({bhp:.0f}) < Pres ({p_reservoir:.0f}): Can produce") else: return ConstraintResult(False, f"Producer BHP ({bhp:.0f}) >= Pres ({p_reservoir:.0f}): Cannot flow", [f"Decrease BHP below {p_reservoir - 100:.0f} psi"]) @staticmethod def check_porosity(porosity: float) -> ConstraintResult: """Porosity must be between 0 and 1.""" if 0 < porosity < 1: if porosity > 0.4: return ConstraintResult(True, f"Porosity {porosity:.2%} valid but unusually high", ["Consider if this is realistic for your reservoir"]) return ConstraintResult(True, f"Porosity {porosity:.2%} valid") else: return ConstraintResult(False, f"Porosity {porosity} out of range [0, 1]", ["Typical values: 0.05 - 0.35"]) @staticmethod def check_permeability_anisotropy(kh: float, kv: float) -> ConstraintResult: """Check permeability anisotropy ratio.""" if kv <= 0 or kh <= 0: return ConstraintResult(False, "Permeability must be positive") ratio = kv / kh if ratio > 1: return ConstraintResult(False, f"kv/kh = {ratio:.2f} > 1 (kv > kh is unusual)", ["Typically kv < kh due to layering", "Verify this is intentional"]) elif ratio < 0.01: return ConstraintResult(True, f"kv/kh = {ratio:.3f} (very low - tight barriers)", ["This will significantly limit vertical flow"]) else: return ConstraintResult(True, f"kv/kh = {ratio:.2f} (reasonable)") # Demo constraint checking constraints = PhysicsConstraints() print("Physics Constraint Checks:") print("=" * 50) # Good values print(constraints.check_saturation(0.3, 0.7, 0.0)) print(constraints.check_pressure_gradient(8500, 3800)) print(constraints.check_porosity(0.22)) print("\nProblematic values:") # Bad values print(constraints.check_saturation(0.5, 0.6, 0.2)) # Sum > 1 print(constraints.check_well_can_flow(3800, 4000, is_injector=False)) # BHP > Pres print(constraints.check_permeability_anisotropy(100, 150)) # kv > kh

Section 3: Z3 Solver Integration¶

Z3 can solve complex constraint systems and find valid parameter combinations.

In [ ]:

if Z3_AVAILABLE:
    def solve_saturations_z3(sw_min: float = 0.2, sw_max: float = 0.8,
                             so_target: float = None) -> Dict[str, float]:
        """Use Z3 to find valid saturation distribution."""
        solver = Solver()
        
        # Define real variables
        Sw = Real('Sw')
        So = Real('So')
        Sg = Real('Sg')
        
        # Constraints
        solver.add(Sw + So + Sg == 1)     # Must sum to 1
        solver.add(Sw >= sw_min)          # Above connate water
        solver.add(Sw <= sw_max)          # Below 1-Sor
        solver.add(So >= 0)               # Non-negative
        solver.add(So <= 1 - sw_min)      # Bounded
        solver.add(Sg >= 0)               # Non-negative
        solver.add(Sg <= 0.3)             # Limit gas (optional)
        
        if so_target:
            solver.add(So >= so_target - 0.05)
            solver.add(So <= so_target + 0.05)
        
        if solver.check() == sat:
            model = solver.model()
            # Convert Z3 rationals to floats
            def to_float(val):
                if hasattr(val, 'as_fraction'):
                    frac = val.as_fraction()
                    return float(frac.numerator) / float(frac.denominator)
                return float(val.as_decimal(10))
            
            return {
                'Sw': to_float(model[Sw]),
                'So': to_float(model[So]),
                'Sg': to_float(model[Sg])
            }
        else:
            return None
    
    # Demo Z3 saturation solver
    print("Z3 Saturation Solver:")
    result = solve_saturations_z3(sw_min=0.25, sw_max=0.75, so_target=0.5)
    if result:
        print(f"  Found valid saturations:")
        print(f"    Sw = {result['Sw']:.4f}")
        print(f"    So = {result['So']:.4f}")
        print(f"    Sg = {result['Sg']:.4f}")
        print(f"    Sum = {sum(result.values()):.4f}")
else:
    print("Z3 not available - skipping Z3 examples")
    print("The constraint validation above still works without Z3")

if Z3_AVAILABLE: def solve_saturations_z3(sw_min: float = 0.2, sw_max: float = 0.8, so_target: float = None) -> Dict[str, float]: """Use Z3 to find valid saturation distribution.""" solver = Solver() # Define real variables Sw = Real('Sw') So = Real('So') Sg = Real('Sg') # Constraints solver.add(Sw + So + Sg == 1) # Must sum to 1 solver.add(Sw >= sw_min) # Above connate water solver.add(Sw <= sw_max) # Below 1-Sor solver.add(So >= 0) # Non-negative solver.add(So <= 1 - sw_min) # Bounded solver.add(Sg >= 0) # Non-negative solver.add(Sg <= 0.3) # Limit gas (optional) if so_target: solver.add(So >= so_target - 0.05) solver.add(So <= so_target + 0.05) if solver.check() == sat: model = solver.model() # Convert Z3 rationals to floats def to_float(val): if hasattr(val, 'as_fraction'): frac = val.as_fraction() return float(frac.numerator) / float(frac.denominator) return float(val.as_decimal(10)) return { 'Sw': to_float(model[Sw]), 'So': to_float(model[So]), 'Sg': to_float(model[Sg]) } else: return None # Demo Z3 saturation solver print("Z3 Saturation Solver:") result = solve_saturations_z3(sw_min=0.25, sw_max=0.75, so_target=0.5) if result: print(f" Found valid saturations:") print(f" Sw = {result['Sw']:.4f}") print(f" So = {result['So']:.4f}") print(f" Sg = {result['Sg']:.4f}") print(f" Sum = {sum(result.values()):.4f}") else: print("Z3 not available - skipping Z3 examples") print("The constraint validation above still works without Z3")

In [ ]:

if Z3_AVAILABLE:
    def solve_well_constraints_z3(p_res: float, min_rate: float, max_rate: float,
                                  pi: float = 10.0) -> Dict[str, float]:
        """Find valid BHP and rate for a producer well.
        
        Uses simplified IPR: q = PI * (Pres - BHP)
        """
        solver = Solver()
        
        BHP = Real('BHP')
        Rate = Real('Rate')
        
        # Constraints
        solver.add(BHP > 0)              # Positive BHP
        solver.add(BHP < p_res)          # Must be below reservoir pressure
        solver.add(Rate == pi * (p_res - BHP))  # IPR equation
        solver.add(Rate >= min_rate)     # Meet minimum rate
        solver.add(Rate <= max_rate)     # Don't exceed max rate
        solver.add(BHP >= 500)           # Minimum flowing pressure
        
        if solver.check() == sat:
            model = solver.model()
            def to_float(val):
                if hasattr(val, 'as_fraction'):
                    frac = val.as_fraction()
                    return float(frac.numerator) / float(frac.denominator)
                return float(str(val))
            
            return {
                'BHP': to_float(model[BHP]),
                'Rate': to_float(model[Rate])
            }
        return None
    
    print("\nZ3 Well Constraint Solver:")
    result = solve_well_constraints_z3(p_res=3800, min_rate=500, max_rate=2000, pi=5.0)
    if result:
        print(f"  Valid operating point found:")
        print(f"    BHP = {result['BHP']:.0f} psi")
        print(f"    Rate = {result['Rate']:.0f} stb/d")
        print(f"    Drawdown = {3800 - result['BHP']:.0f} psi")

if Z3_AVAILABLE: def solve_well_constraints_z3(p_res: float, min_rate: float, max_rate: float, pi: float = 10.0) -> Dict[str, float]: """Find valid BHP and rate for a producer well. Uses simplified IPR: q = PI * (Pres - BHP) """ solver = Solver() BHP = Real('BHP') Rate = Real('Rate') # Constraints solver.add(BHP > 0) # Positive BHP solver.add(BHP < p_res) # Must be below reservoir pressure solver.add(Rate == pi * (p_res - BHP)) # IPR equation solver.add(Rate >= min_rate) # Meet minimum rate solver.add(Rate <= max_rate) # Don't exceed max rate solver.add(BHP >= 500) # Minimum flowing pressure if solver.check() == sat: model = solver.model() def to_float(val): if hasattr(val, 'as_fraction'): frac = val.as_fraction() return float(frac.numerator) / float(frac.denominator) return float(str(val)) return { 'BHP': to_float(model[BHP]), 'Rate': to_float(model[Rate]) } return None print("\nZ3 Well Constraint Solver:") result = solve_well_constraints_z3(p_res=3800, min_rate=500, max_rate=2000, pi=5.0) if result: print(f" Valid operating point found:") print(f" BHP = {result['BHP']:.0f} psi") print(f" Rate = {result['Rate']:.0f} stb/d") print(f" Drawdown = {3800 - result['BHP']:.0f} psi")

Section 4: Constraint Validator Class¶

A unified interface for validating deck parameters.

In [ ]:

@dataclass
class DeckParameters:
    """Parameters extracted from a deck for validation."""
    # Grid
    nx: int = 10
    ny: int = 10
    nz: int = 5
    dx: float = 100.0  # ft
    dy: float = 100.0
    dz: float = 20.0
    
    # Rock properties
    porosity: float = 0.2
    permx: float = 100.0  # md
    permy: float = 100.0
    permz: float = 10.0
    
    # Initial conditions
    pressure: float = 3800.0  # psi
    sw_init: float = 0.25
    depth: float = 8500.0  # ft
    
    # Wells
    wells: List[Dict] = field(default_factory=list)

class DeckValidator:
    """Validates deck parameters against physics constraints."""
    
    def __init__(self):
        self.physics = PhysicsConstraints()
        self.results: List[ConstraintResult] = []
    
    def validate(self, params: DeckParameters) -> Tuple[bool, List[ConstraintResult]]:
        """Run all validations on deck parameters."""
        self.results = []
        
        # Grid checks
        if params.nx * params.ny * params.nz > 1_000_000:
            self.results.append(ConstraintResult(False, 
                f"Grid too large: {params.nx}x{params.ny}x{params.nz} = {params.nx*params.ny*params.nz:,} cells",
                ["Consider coarsening for initial runs"]))
        else:
            self.results.append(ConstraintResult(True,
                f"Grid size OK: {params.nx*params.ny*params.nz:,} cells"))
        
        # Rock properties
        self.results.append(self.physics.check_porosity(params.porosity))
        self.results.append(self.physics.check_permeability_anisotropy(
            params.permx, params.permz))
        
        # Pressure checks
        self.results.append(self.physics.check_pressure_positive(
            params.pressure, "Initial pressure"))
        self.results.append(self.physics.check_pressure_gradient(
            params.depth, params.pressure))
        
        # Saturation
        so_init = 1 - params.sw_init  # Assuming no gas
        self.results.append(self.physics.check_saturation(
            params.sw_init, so_init, 0.0))
        
        # Well checks
        for well in params.wells:
            # Check well location within grid
            if well.get('i', 0) > params.nx or well.get('j', 0) > params.ny:
                self.results.append(ConstraintResult(False,
                    f"Well {well['name']} location ({well['i']},{well['j']}) outside grid"))
            
            # Check BHP
            if 'bhp' in well:
                self.results.append(self.physics.check_well_can_flow(
                    params.pressure, well['bhp'], 
                    well.get('type') == 'injector'))
        
        # Overall result
        all_passed = all(r.satisfied for r in self.results)
        return all_passed, self.results
    
    def print_report(self):
        """Print validation report."""
        passed = sum(1 for r in self.results if r.satisfied)
        total = len(self.results)
        
        print(f"\nValidation Report: {passed}/{total} checks passed")
        print("=" * 50)
        
        for r in self.results:
            print(r)
            if r.fixes:
                for fix in r.fixes:
                    print(f"    -> {fix}")

# Demo validation
params = DeckParameters(
    nx=20, ny=20, nz=5,
    porosity=0.22,
    permx=150, permy=150, permz=15,
    pressure=3800,
    sw_init=0.25,
    depth=8500,
    wells=[
        {'name': 'PROD1', 'i': 10, 'j': 10, 'bhp': 2000, 'type': 'producer'},
        {'name': 'INJ1', 'i': 1, 'j': 1, 'bhp': 4500, 'type': 'injector'},
    ]
)

validator = DeckValidator()
all_valid, results = validator.validate(params)
validator.print_report()

@dataclass class DeckParameters: """Parameters extracted from a deck for validation.""" # Grid nx: int = 10 ny: int = 10 nz: int = 5 dx: float = 100.0 # ft dy: float = 100.0 dz: float = 20.0 # Rock properties porosity: float = 0.2 permx: float = 100.0 # md permy: float = 100.0 permz: float = 10.0 # Initial conditions pressure: float = 3800.0 # psi sw_init: float = 0.25 depth: float = 8500.0 # ft # Wells wells: List[Dict] = field(default_factory=list) class DeckValidator: """Validates deck parameters against physics constraints.""" def __init__(self): self.physics = PhysicsConstraints() self.results: List[ConstraintResult] = [] def validate(self, params: DeckParameters) -> Tuple[bool, List[ConstraintResult]]: """Run all validations on deck parameters.""" self.results = [] # Grid checks if params.nx * params.ny * params.nz > 1_000_000: self.results.append(ConstraintResult(False, f"Grid too large: {params.nx}x{params.ny}x{params.nz} = {params.nx*params.ny*params.nz:,} cells", ["Consider coarsening for initial runs"])) else: self.results.append(ConstraintResult(True, f"Grid size OK: {params.nx*params.ny*params.nz:,} cells")) # Rock properties self.results.append(self.physics.check_porosity(params.porosity)) self.results.append(self.physics.check_permeability_anisotropy( params.permx, params.permz)) # Pressure checks self.results.append(self.physics.check_pressure_positive( params.pressure, "Initial pressure")) self.results.append(self.physics.check_pressure_gradient( params.depth, params.pressure)) # Saturation so_init = 1 - params.sw_init # Assuming no gas self.results.append(self.physics.check_saturation( params.sw_init, so_init, 0.0)) # Well checks for well in params.wells: # Check well location within grid if well.get('i', 0) > params.nx or well.get('j', 0) > params.ny: self.results.append(ConstraintResult(False, f"Well {well['name']} location ({well['i']},{well['j']}) outside grid")) # Check BHP if 'bhp' in well: self.results.append(self.physics.check_well_can_flow( params.pressure, well['bhp'], well.get('type') == 'injector')) # Overall result all_passed = all(r.satisfied for r in self.results) return all_passed, self.results def print_report(self): """Print validation report.""" passed = sum(1 for r in self.results if r.satisfied) total = len(self.results) print(f"\nValidation Report: {passed}/{total} checks passed") print("=" * 50) for r in self.results: print(r) if r.fixes: for fix in r.fixes: print(f" -> {fix}") # Demo validation params = DeckParameters( nx=20, ny=20, nz=5, porosity=0.22, permx=150, permy=150, permz=15, pressure=3800, sw_init=0.25, depth=8500, wells=[ {'name': 'PROD1', 'i': 10, 'j': 10, 'bhp': 2000, 'type': 'producer'}, {'name': 'INJ1', 'i': 1, 'j': 1, 'bhp': 4500, 'type': 'injector'}, ] ) validator = DeckValidator() all_valid, results = validator.validate(params) validator.print_report()

Section 5: Constraint Propagation¶

When one parameter changes, we may need to adjust others to maintain consistency.

In [ ]:

class ConstraintPropagator:
    """Propagates constraints to maintain consistency."""
    
    @staticmethod
    def adjust_saturations(sw: float, so: float, sg: float) -> Tuple[float, float, float]:
        """Normalize saturations to sum to 1.0."""
        total = sw + so + sg
        if total == 0:
            return 0.25, 0.75, 0.0  # Default
        return sw/total, so/total, sg/total
    
    @staticmethod
    def suggest_bhp_for_rate(p_res: float, target_rate: float, 
                             pi: float, is_injector: bool) -> float:
        """Calculate BHP needed to achieve target rate.
        
        Uses: q = PI * |Pres - BHP|
        """
        delta_p = target_rate / pi
        if is_injector:
            bhp = p_res + delta_p
        else:
            bhp = p_res - delta_p
        return max(bhp, 100)  # Minimum BHP
    
    @staticmethod
    def suggest_grid_refinement(area_acres: float, target_cells: int = 10000
                                ) -> Tuple[int, int, float]:
        """Suggest grid dimensions for target cell count.
        
        Returns (nx, ny, dx_ft).
        """
        area_ft2 = area_acres * 43560  # Convert acres to ft2
        
        # Assume square grid
        side = math.sqrt(area_ft2)
        n = int(math.sqrt(target_cells))
        dx = side / n
        
        return n, n, dx

# Demo propagation
prop = ConstraintPropagator()

print("Constraint Propagation Examples:")
print("=" * 50)

# Normalize bad saturations
sw, so, sg = prop.adjust_saturations(0.5, 0.6, 0.2)
print(f"\nNormalized saturations: Sw={sw:.3f}, So={so:.3f}, Sg={sg:.3f}")
print(f"Sum: {sw + so + sg:.4f}")

# Calculate BHP for target rate
bhp = prop.suggest_bhp_for_rate(p_res=3800, target_rate=1000, pi=5.0, is_injector=False)
print(f"\nBHP for 1000 stb/d at PI=5: {bhp:.0f} psi")

# Suggest grid dimensions
nx, ny, dx = prop.suggest_grid_refinement(area_acres=640, target_cells=10000)
print(f"\nGrid for 640 acres (1 section), ~10000 cells:")
print(f"  {nx} x {ny} cells, dx = {dx:.0f} ft")

class ConstraintPropagator: """Propagates constraints to maintain consistency.""" @staticmethod def adjust_saturations(sw: float, so: float, sg: float) -> Tuple[float, float, float]: """Normalize saturations to sum to 1.0.""" total = sw + so + sg if total == 0: return 0.25, 0.75, 0.0 # Default return sw/total, so/total, sg/total @staticmethod def suggest_bhp_for_rate(p_res: float, target_rate: float, pi: float, is_injector: bool) -> float: """Calculate BHP needed to achieve target rate. Uses: q = PI * |Pres - BHP| """ delta_p = target_rate / pi if is_injector: bhp = p_res + delta_p else: bhp = p_res - delta_p return max(bhp, 100) # Minimum BHP @staticmethod def suggest_grid_refinement(area_acres: float, target_cells: int = 10000 ) -> Tuple[int, int, float]: """Suggest grid dimensions for target cell count. Returns (nx, ny, dx_ft). """ area_ft2 = area_acres * 43560 # Convert acres to ft2 # Assume square grid side = math.sqrt(area_ft2) n = int(math.sqrt(target_cells)) dx = side / n return n, n, dx # Demo propagation prop = ConstraintPropagator() print("Constraint Propagation Examples:") print("=" * 50) # Normalize bad saturations sw, so, sg = prop.adjust_saturations(0.5, 0.6, 0.2) print(f"\nNormalized saturations: Sw={sw:.3f}, So={so:.3f}, Sg={sg:.3f}") print(f"Sum: {sw + so + sg:.4f}") # Calculate BHP for target rate bhp = prop.suggest_bhp_for_rate(p_res=3800, target_rate=1000, pi=5.0, is_injector=False) print(f"\nBHP for 1000 stb/d at PI=5: {bhp:.0f} psi") # Suggest grid dimensions nx, ny, dx = prop.suggest_grid_refinement(area_acres=640, target_cells=10000) print(f"\nGrid for 640 acres (1 section), ~10000 cells:") print(f" {nx} x {ny} cells, dx = {dx:.0f} ft")

Summary¶

In this tutorial, we learned:

1. Unit Systems: FIELD vs METRIC conversions for ECLIPSE
2. Physics Constraints: Saturation, pressure, permeability rules
3. Z3 Solver: Find valid parameter combinations automatically
4. Deck Validation: Comprehensive parameter checking
5. Constraint Propagation: Adjust related parameters together

Key Insight: Catching physics errors before simulation saves hours of debugging!

Next Tutorial: 06_Deck_Generator.ipynb - Template-based deck generation