Machine Determination of Better Representations of Multiphase Equations of State for Subsurface Flow Simulation¶

or¶

Solution of Differential Algebraic Equations with Local Constraints on the Latent Space of an Autoencoder¶

Alejandro Francisco Queiruga¶

Lawrence Berkeley National Lab¶

Outline¶

What problems am I trying to solve?
Phase transition difficulties
Rephrase our equations
The classical pendulum
Plugging autoencoders into equations
Unsupervised learning for phases
Batch training and testing of simulations

What kinds of problems?¶

Single phase Darcy's law is easy:

$\begin{equation} M \partial_t p = \nabla \cdot \frac{k}{\mu} \nabla p \end{equation}$

Multiphase Transport Equations¶

Mass balances and energy balance:

$\begin{align} \partial_t \rho_{H_2O} & = \nabla \cdot \frac{k_1(p,T)}{\mu} \nabla p \\ \partial_t \rho_{CO_2} & = \nabla \cdot \frac{k_2(p,T)}{\mu} \nabla p \\ \partial_t \rho_{CH_4} & = \nabla \cdot \frac{k_3(p,T)}{\mu} \nabla p \\ \partial_t \sum \rho u & = \nabla \cdot q \end{align}$

Complicated equations and empirical fits for everything
Mass fraction and multiple phases for each material:

$\begin{equation} \rho = X_{H_2O}\left(S_{gas} \rho_{gas} + S_{liq} \rho_{liq} + S_{ice} \rho_{ice} + S_{crit} \rho_{crit} \right) \end{equation}$

(Not shown: continuum mechanics and chemical reactions!)

How do we solve this? State Machines¶

Gridblocks have two unknowns X1 X2 and an additional phase tag.

Gas: $X = \{p,T\}$
$\rho_{gas}(p,T) =$ one fit
Liquid: $X = \{p,T\}$
$\rho_{liq}(p,T) =$ another fit
Liq-Gas: $X = \{S_{gas},T\}$
$\rho_{mix}(S,T) = S \rho_{gas}(p^*(T),T)+ (1-S) \rho_{liq}(p^*(T),T)$

Need potentially different variables for each state
Equilibria usually use a phase saturation $S$
Need logic to handle state of the material
Limits discretization choice

State machine for phase transitions:¶

Formulate balance equations as $\begin{equation} R(X ; phase) = 0 \end{equation}$ Solve the differentiable part: $\begin{equation} \frac{\partial R}{\partial X}\Delta X^{k+1/2} = -R(X^k ; phase^k) \end{equation}$ and iterate on the nondifferential part: $\begin{equation} phase^{k+1},X^{k+1} = statemachine(X^{k+1/2} ; phase^k) \end{equation}$

switch phase_old:
  case gas:
    if p,T crossing boundary:
      phase_new = liquid_gas
  case liquid:
    if p,T crossing boundary:
      phase_new = liquid_gas
  case liquid_gas:
    if S_gas >= 1:
      phase_new = gas
    if S_liquid >= 1:
      phase_new = gas

No quadratic convergence
Lots of bug ridden coding!
Easily unstable! Need a lot of hacks, like overshoots: $\begin{equation} p_{new} = (1+10^{-6})p_{old} \end{equation}$

How do we solve this better?¶

Take a step back:

Forget about phases and states.
Material is a manifold of possible intensive properties.
It's just empirical data.
We have observations and expect to stay on these observations.
We just need any representation of this constraint.

A Constrained Differential Algebraic Equation¶

Solve for $\rho(t), \, h(t), \, p(t),\, \text{and}\, T(t)$ satisfying: $\begin{align} \partial_t \rho & = \nabla \cdot \mathbf{k}\nabla p + r & \quad \}\text{mass balance}\\ \partial_t (\rho h-p) & = \nabla \cdot \mathbf{k'}\nabla T + s & \quad \}\text{energy balance}\\ \end{align}$ such that they lie on the material EOS, $\begin{equation} eos(\rho,h,p,T) = 0 \quad \text{or} \quad \{\rho,h,p,T\} \in \{ eos \} \end{equation}$ We just made the problem harder.

The Archetypical Constraint Problem: The Pendulum¶

Solve for $x(t),y(t)$ stationary on $\begin{equation} \mathcal{L}(x,y) = \frac{1}{2}m\left(\dot{x}^2 + \dot{y}^2\right) - m g y \end{equation}$ subject to $\begin{equation} x^2 + y^2 = R^2 \end{equation}$ Could do spring penalty, Lagrange multiplier, or...

The Lagrangian mechanics formulation lets us introduce a parameterization: $\begin{equation} x = R \cos \theta, \quad y = R \sin \theta \end{equation}$ to plug into the Lagrangian: $\begin{equation} \mathcal{L}(\theta) = \frac{R^2}{2} \dot{\theta}^2 + g R \sin \theta \end{equation}$ The equations of motion in terms of $\theta$ are then: $\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial L}{\partial \dot{\theta}} = \frac{\partial L}{\partial \theta} \end{equation}$

The Pendulum¶

The choice of $\theta$ is arbitrary. There are infinimately choices.

Suppose we only have data: (Imagine the EOS plot I showed earlier.)

In [9]:

import numpy as np
from matplotlib import pylab as plt
%matplotlib inline
theta = np.linspace(-np.pi*1.25,0.25*np.pi,50, dtype=np.float32)
data = np.vstack([np.cos(theta), np.sin(theta)]).T
plt.plot(data[:,0],data[:,1],'+')
plt.axis('square');

We're looking for a mapping from the set of points onto one variable.

Autoencoders¶

autoencoder

$\theta$ is one instance of an autoencoder. We can look for one without the geometric intuition.

Autoencoders¶

Solve the minimization problem on data set $x$ for parameters $a$ $\begin{equation} \min_a \sum_x \left( x - D(E(x;a);a) \right)^2 \end{equation}$ where $q = E(x)$ with len(q)<len(x).

$q$ represents a point on the Latent Space
Pick the dimension of $q$ based on prior intuition of the manifold (or tune it)
Pendulum is a 1D manifold, so look for a 1D $q$
It's very smooth and continuous so let's use polynomials:

$\begin{equation} \left\{\begin{array}{c} x\\y\end{array}\right\} \rightarrow W^1 x^n \rightarrow q \rightarrow W^2 x^n \rightarrow \left\{\begin{array}{c} x\\y\end{array}\right\} \end{equation}$

Caveats: Topology¶

We can't encode a full circle.
The inner space of this arbitrary architecture is a disk.
We would need multiple charts to span the circle (more dimensions to $q$ + chart selection), or pick an encoder that maps onto the right topology.

bad pendulum

This is why I chopped off the top.
Fortunately EOSes are disks so we don't need to solve this problem.

Solving the Lagrangian: Just do it.¶

Treat the autoencoder the same as $x(\theta)$ , plugging in $\begin{equation} L(x(q),v(q,\dot{q})) \end{equation}$ and build the DAE we know from physics: $\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial L}{\partial \dot{q}} = \frac{\partial L}{\partial q} \end{equation}$ where we get the velocity with the chain rule: $\begin{equation} v = \frac{\partial x}{\partial q} \dot{q} \end{equation}$

Implement it in the TensorFlow syntax (it's a general language!):

def L(i_q,i_qdot):                               # Input placeholders
    x = au.decode(i_q)                           # Trained decoder
    dxdq = atu.vector_gradient(x,i_q)            # chain rule
    v = tf.einsum("ikj,ij->ik",dxdq,i_qdot)      # velocity
    L = 0.5*tf.einsum("ij,ij->i",v,v) - g*x[:,1] # action
    return x,v, tf.expand_dims(L,-1)             # return TF graphs

Then plug it into a Runge Kutta (trapezoidal) to solve with Newton's method: $\begin{equation} \frac{\partial L}{\partial \dot{q}}_i - \Delta t \frac{\partial L}{\partial q}_i = \frac{\partial L}{\partial \dot{q}}_0 + \Delta t \frac{\partial L}{\partial q}_0 \end{equation}$

Solutions of the Pendulum¶

We solved for $q(t)$ and $\dot{q}(t)$
Decode $x,y$ and $v_x,v_y$ as a post-process

Back to Multiphase EOS¶

Solve for $q_1(t)$ and $q_2(t)$ such that: $\begin{align} \partial_t \rho(q_1,q_2) & = \nabla \cdot \mathbf{k} \nabla p(q_1,q_2) + r \\ \partial_t \rho h(q_1,q_2)-p & = \nabla \cdot \mathbf{k'}\nabla T(q_1,q_2) + s \end{align}$ where $\rho(q_1,q_2)$ etc. are the back of an autoencoder: $\begin{equation} \left\{ \begin{array}{c} T\\ p\\ \rho\\ h \end{array}\right\} \rightarrow E \rightarrow \left\{ \begin{array}{c} q_1\\q_2 \end{array} \right\}\rightarrow D \rightarrow \left\{ \begin{array}{c} T\\ p\\ \rho\\ h \end{array}\right\} \end{equation}$ Using the autoencoder to automatically remove the local constraint: $\begin{equation} eos(D(q_1,q_2)) = 0 \, \forall \,q_1,q_2 \end{equation}$

Only training the material representation, not the balance laws.

Solving the equations:¶

First, we had a constrained DAE: $\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} m(T,p,\rho,h) &= r(T,p,\rho,h)\\ eos(T,p,\rho,h) &= 0 \end{align}$ with no good representation of $eos$ everywhere.

The autoencoder parameterizes the constraint : $\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} m(D(q)) = r(D(q)) \end{equation}$ Apply a Runge-Kutta to get: (BW Euler) $\begin{equation} 0 = R(q^i) = m(q^{i}) - m(q^{0}) - \Delta t \, r(q^i) \end{equation}$ By design $R(q)$ is smooth enough to solve with Newton's method: $\begin{equation} \frac{\partial R}{\partial q} \Delta q^{k+1} = R(q^k) \end{equation}$ (Well, it's a hard nonlinear problem, but a smooth one!)

Methodology¶

Make a database of T,p,ρ,h,X1,...
- Piece together empirical fits for each phase from literature
- (Experimental data in the future)
Normalize (and log(p)) the database, shuffle it
- $\log(p)$ distributes the low- $p$ phases evenly w.r.t. high- $p$ phases
Train the autoencoder
- Don't use phase labels
- Batch test architectures
Load the models and generate physics code
Verify and grade architectures on tests
- Need more than autoencoder mean-squared-error
- Differentiability and numerical stability in simulation
- Evaluation speed (billions of times in a simulation!)
Pick best one to package into production code

Can be automated after step 1.

$T,p,rho,h$ are smooth in $q_1,q_2$ ¶

In [4]:

IFrame("slides/water_slgc_logp_64_qq_plots.html",1200,650)

Out[4]:

Testing the Method¶

Three different extents:

Linear Equation of State
- Test!
- Reduces to single phase Darcy's law problem (slide 1)
- $p = 10^5+[-10^3, 10^3] Pa,\quad T = [ 19, 21 ] ^o C$
Water Liquid-Gas Regime
- One phase boundary
- $p = [100,5\times 10^5] Pa, \quad T = [274,594] K$
Water Solid-Liquid-Gas-Supercritical Regimes
- No linear mapping to latent space
- Entire span
- $p = [6\times 10^{-6},3\times 10^8]Pa, \quad T = [150,1000] K$

.	name	train_goal	succ	min_err	max_err	run_time
x	Classifying_2,4,12,24,sigmoid_0	3.14e-04	3	2.61e-03	3.90e+04	138.3
x	Poly_3,7	5.00e-05	2	6.55e-05	1.45e-01	259.3
x	Classifying_2,6,18,36,sigmoid	1.39e-05	2	6.21e-10	1.67e-01	377.0
x	Deep_2,(12),3,(12,12,12)	2.54e-05	2	7.63e-06	2.09e-01	59.6
x	Deep_1,(12),1,(12,12,12)	2.92e-05	2	1.67e-08	3.38e+04	54.7
x	Classifying_2,1,12,24,sigmoid	8.82e-05	3	1.05e-05	4.05e+04	36.9
x	Classifying_2,4,12,24,sigmoid	-1.00e+00	2	1.22e-09	1.51e-01	158.3
x	Classifying_2,5,12,24,sigmoid	1.37e-05	3	2.44e-05	4.03e+04	227.5
x	Classifying_2,3,12,24,sigmoid	1.75e-05	2	9.05e-06	1.26e-01	88.4
x	Deep_1,(),1,(12,12,12)	7.31e-02	2	4.27e-06	4.08e+03	71.9
x	Poly_2,7	1.12e-04	3	8.47e-09	4.03e+04	253.6
x	Classifying_2,6,12,24,sigmoid	6.16e-05	2	4.56e-03	4.05e+04	314.2
x	Classifying_2,6,24,48,sigmoid	2.83e-06	2	6.34e-10	1.70e-01	368.4