<jupyter_start><jupyter_text>Pressure deep learning surrogate model<jupyter_code>gpu_info = !nvidia-smi
gpu_info = '\n'.join(gpu_info)
if gpu_info.find('failed') >= 0:
  print('Select the Runtime > "Change runtime type" menu to enable a GPU accelerator, ')
  print('and then re-execute this cell.')
else:
  print(gpu_info)
from psutil import virtual_memory
ram_gb = virtual_memory().total / 1e9
print('Your runtime has {:.1f} gigabytes of available RAM\n'.format(ram_gb))

if ram_gb < 20:
  print('To enable a high-RAM runtime, select the Runtime > "Change runtime type"')
  print('menu, and then select High-RAM in the Runtime shape dropdown. Then, ')
  print('re-execute this cell.')
else:
  print('You are using a high-RAM runtime!')<jupyter_output>Your runtime has 27.3 gigabytes of available RAM

You are using a high-RAM runtime!<jupyter_text>1. Introduction In this notebook, we demonstrate the workflow to train a surrogate deep learning model which can predict pressure and saturation using the porosity and permeability map.Being accurate, fluid flow simulation through porous media may be really time consuming (several hours for a single run) and we have to perform hundreds of thousands of run to calibrate the models using observation data and to correctly quantifying uncertainty. / In oil industry a fast and accurate model are required for three-phase flow simulation on a 3D computational grid.Here we will be using the deep learning platform and library Tensorflow and Keras , due to its flexibility and applications to research.The Residual U-Net in combination with LSTM have been used as deep learning model. The modeule of this model network has been written as modeules. So, we need to upload the modules to use them.<jupyter_code>from google.colab import files
files.upload()
from tensorflow.keras.callbacks import EarlyStopping
import numpy as np
import unet_att as vae_util_att
import unet_uae as vae_util
import os
os.environ["CUDA_VISIBLE_DEVICES"]="0"
from tensorflow.keras import layers
from tensorflow.keras.models import load_model, Model
from tensorflow.python.keras import backend as K
from keras.optimizers import Adam
import tensorflow as tf 
config = tf.compat.v1.ConfigProto()
config.gpu_options.allow_growth = True
config.gpu_options.per_process_gpu_memory_fraction = 0.3
sess = tf.compat.v1.Session(config=config)
K.set_session(sess)
from google.colab import drive
drive.mount('/content/gdrive/')<jupyter_output>Mounted at /content/gdrive/<jupyter_text>2. The dataset First we start with loading in the dataset and evaluating its content.  2.1. Loading Datasetto do this, since the data are as .npy format we need to load them using load function. the following function was written to load data<jupyter_code>def load_data(path):
    """
    Function to load datasets in format .NPY
    
    Parameter:
    ----------
    path : string
        The absolute path of where data saved in local system
        
    Return:
    ----------
    loaded_data : ndarray
        The data which was loaded
    """
    loaded_data = np.load(path)
    return loaded_data
# define the absolute path of training datatsat
path_perm = '/content/gdrive/MyDrive/perm.npy'
path_press = '/content/gdrive/MyDrive/pressure.npy'
# use load_data function nd above path to loading data
X_data= load_data(path_perm)
target_data = load_data(path_press)<jupyter_output><empty_output><jupyter_text>2.2. Preprocessing data (Normalization and standardization)  Mean and standard deviation, Max and min of training and validation set were calcualted. Then pressure data standardized and normalized using these statistical values.<jupyter_code># Two common methods for feature scaling is : 1-Normalization & 2-Standardaisation

def normalize(data):
    """
    this function used for Max-Min Normalization (Min-Max scaling) by re-scaling
    features with a distribution value between 0 and 1. For every feature,the minimum
    value of that feature gets transformed into 0, and the maximum value 
    gets transformed into 1
    
    Parameter:
    ----------
    data : ndarray
        The numpy array which we want to normalize
        
    Return:
    ----------
    norm_data : ndarray
        The normalized data which transformed into 0 and 1
    """
    max_p = np.max(data[:, :, :, :])
    min_p = np.min(data[:, :, :, :])
    norm_data = (data - min_p)/(max_p - min_p)
    return norm_data

def standardize(data):
    """
    this function used for rescaling faetures to ensure the mean
    and the standard deviation to be 0 and 1, respectively.
    
    Parameter:
    ----------
    data : ndarray
        The numpy array which we want to normalize
        
    Return:
    ----------
    data : ndarray
        The standardized data which the mean
    and the standard deviation to be 0 and 1
    """
    data_mean = np.mean(data[:, :, :, :], axis = 0, keepdims = True)
    data_std = np.std(data[:, :, :, :], axis = 0, keepdims = True)
    std_data = (data - data_mean)/(data_std)
    return std_data
# Normalize data using abov normalize function
train_nr = 2250
test_nr = 750
p_t_mean = np.mean(target_data[:train_nr, ...], axis = 0, keepdims = True)
target_data = target_data - p_t_mean
print('max p is ', np.max(target_data[:train_nr, ...]), ', min p is ', np.min(target_data[:train_nr, ...]))
max_p = np.max(target_data[:train_nr, ...])
min_p = np.min(target_data[:train_nr, ...])
target_data = (target_data - min_p)/(max_p -min_p) - 0.5
print('max p is ', np.max(target_data), ', min p is ', np.min(target_data))
print('max p train is ', np.max(target_data[:train_nr, ...]), ', min p train is ', np.min(target_data[:train_nr, ...]))<jupyter_output>max p is  516.6753418240017 , min p is  -135.6267126736111
max p is  0.5039995270719013 , min p is  -0.500166926837037
max p train is  0.5 , min p train is  -0.5<jupyter_text>3. Training the model To train the model we defined and import the modules. Since we are dealing with encoder/decoder model we will be using loss and relative error metric as it can give better understanding of model by comparing the result from decoder with the target data.<jupyter_code>def vae_loss(x, t_decoded):
    '''Total loss for the plain UAE'''
    return K.mean(reconstruction_loss(x, t_decoded))


def reconstruction_loss(x, t_decoded):
    '''Reconstruction loss for the plain UAE'''

    return K.sum((K.batch_flatten(x) - K.batch_flatten(t_decoded)) ** 2, axis=-1)

def relative_error(x, t_decoded):
    return K.mean(K.abs(x - t_decoded) / x)<jupyter_output><empty_output><jupyter_text>summary of model layers and parameters.<jupyter_code>input_shape=(100, 100, 2)
depth = 10
vae_model,_ = vae_util.create_vae(input_shape, depth)
vae_model_att,_ = vae_util_att.create_vae(input_shape, depth)
#vae_model.summary()
depth = 10
nr = X_data.shape[0]
train_nr = 2250
test_nr = 750
train_x = np.concatenate([X_data[:train_nr,[0], ...],target_data[:train_nr,[0], ...]], axis = 1)
train_y = target_data[:train_nr, ...]

test_x = np.concatenate([X_data[nr-test_nr:,[0], ...], target_data[nr-test_nr:, [0], ...]], axis = 1)
test_y = target_data[nr-test_nr:,...]


train_x = train_x.transpose(0,2,3,1)
train_y = train_y[:,:,:,:,None]
test_x = test_x.transpose(0,2,3,1)
test_y = test_y[:,:,:,:,None]
#test_y = test_y.transpose(0,2,3,1)
print('train_x shape is ', train_x.shape)
print('train_y shape is ', train_y.shape)
print('test_x shape is ', test_x.shape)
print('test_y shape is ', test_y.shape)<jupyter_output>train_x shape is  (2250, 100, 100, 2)
train_y shape is  (2250, 10, 100, 100, 1)
test_x shape is  (750, 100, 100, 2)
test_y shape is  (750, 10, 100, 100, 1)<jupyter_text>5. Evaluation the model With the final model trained to a great accuracy and minimum loss we need to evaluate the performance of the model. Since we do not have access to the labels of the test data, we will use a version that is only trained on 80 % of the data.<jupyter_code>vae_model.load_weights('/content/gdrive/MyDrive/saved_models/pressure_le1-4.h5')
vae_model_att.load_weights('/content/gdrive/MyDrive/saved_models/pressure_10_points_attentionunet_1e-4batch8.h5')
p_pred_list = []
p_pred_att_list = []
test_batch_size = 50 
batch_num = int(test_nr / test_batch_size)
for i in range(batch_num):
    pred = vae_model.predict(test_x[i*test_batch_size: (i+1)*test_batch_size, ...])
    pred_att = vae_model_att.predict(test_x[i*test_batch_size: (i+1)*test_batch_size, ...])
    p_pred_list.append(pred[:, :, :, :, 0])
    p_pred_att_list.append(pred_att[:, :, :, :, 0])

p_pred = np.concatenate(p_pred_list, axis = 0)
p_pred_att = np.concatenate(p_pred_att_list, axis = 0)
print(p_pred.shape)
p_true = test_y[:,:,:,:,0]
print(p_true.shape)
# convert pressure back to correct range
p_true = ((max_p - min_p) * (p_true)+ min_p)  + p_t_mean[:,:, ...]
p_pred = ((max_p - min_p) * (p_pred) + min_p) + p_t_mean[:,:, ...]
p_pred_att = ((max_p - min_p) * (p_pred_att) + min_p) + p_t_mean[:,:, ...]
print(p_pred.shape)
print(p_pred_att.shape)
p_true = p_true.transpose(0,2,3,1)
p_pred = p_pred.transpose(0,2,3,1)
p_pred_att = p_pred_att.transpose(0,2,3,1)
print(p_true.shape)
print(p_pred.shape)
import matplotlib.pyplot as plt
%matplotlib inline
plt.set_cmap('jet')

ind_case = [100]
num_case = len(ind_case)
for k in ind_case:
    for i_tstep in range(1):
        plt.figure(figsize=(10,7))
        # True permeability
        plt.subplot(3, 1, 1)
        plt.imshow(test_x[k, ::-1, :, 0])
        cbar = plt.colorbar(fraction=0.046)
        cbar.ax.tick_params(labelsize=16)

        plt.tight_layout()
        # plt.subplot(3, 1, 2)
        # plt.imshow(test_x[k, ::-1, :, 0].reshape((100,100)), vmin=-0.09, vmax=1.07)
        # #plt.set_clim([-0.09, 1.07])
        # cbar = plt.colorbar(fraction=0.046)
        # cbar.ax.tick_params(labelsize=16)
        # plt.tight_layout()
        # # The difference between True permeability after prediction
        # plt.subplot(3, 1, 3)
        # plt.imshow(test_x[k, ::-1, :, 0].reshape((100,100)) - test_x[k, ::-1, :, i_tstep])
        # cbar = plt.colorbar(fraction=0.046)
        # cbar.ax.tick_params(labelsize=16)
        # #plt.clim([-0.2, 0.2])
        # plt.tight_layout()
        #plt.savefig('ind_case_%d_logk'%(k),bbox_inches='tight')
plt.show()
for k in ind_case:
    plt.figure(figsize=(20,15))
    tstep = [0,3,6,9]
    for i_tstep in range(len(tstep)):
        plt.subplot(5, 4, i_tstep+1)
        plt.imshow(p_true[k, i_tstep, :, :].reshape((100,100)), vmin = 30, vmax = 200)
        plt.colorbar(fraction=0.046)

        plt.subplot(5, 4, i_tstep+1+4)
        plt.imshow(p_pred[k, i_tstep, :, :],  vmin = 30, vmax = 200)
        plt.colorbar(fraction=0.046)

        plt.subplot(5, 4, i_tstep+1+8)
        plt.imshow(p_pred_att[k, i_tstep, :, :],  vmin = 30, vmax = 200)
        plt.colorbar(fraction=0.046)
        
        plt.subplot(5, 4, i_tstep+1+12)
        plt.imshow(p_true[k, i_tstep, :, :].reshape((100,100)) - p_pred[k, i_tstep, :, :], vmin = -10, vmax = 10)
        plt.colorbar(fraction=0.046)

        plt.subplot(5, 4, i_tstep+17)
        plt.imshow(p_true[k, i_tstep, :, :].reshape((100,100)) - p_pred_att[k, i_tstep, :, :], vmin = -10, vmax = 10)
        plt.colorbar(fraction=0.046)
        #plt.clim([-0.2, 0.2])
        #plt.tight_layout()

plt.savefig('/content/gdrive/MyDrive/Gridsearch_images/compare_pressure_map.png')
p_relative_error_time_att = np.mean(np.abs(p_true - p_pred_att) / (np.max(p_true, axis = (1, 2), keepdims =True) - np.min(p_true, axis = (1, 2), keepdims =True)), axis = (0, 1, 2))
p_relative_error_time = np.mean(np.abs(p_true - p_pred) / (np.max(p_true, axis = (1, 2), keepdims =True) - np.min(p_true, axis = (1, 2), keepdims =True)), axis = (0, 1, 2))
plt.figure(figsize=(8,6))
time = [3, 6, 12, 18, 24, 30, 36, 42, 51, 60]
plt.plot(time, p_relative_error_time, color='blue', linewidth=5, label='RR U-net')
plt.plot(time, p_relative_error_time_att, color='red', linewidth=5, label='A RR U-net')

plt.xlabel("Production time (months)")
plt.ylabel("Mean relative error")
plt.legend()
plt.savefig('/content/gdrive/MyDrive/Gridsearch_images/pressure_rel-error.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13352015165251793, 'fraction_numerical': 0.035522953927305694, 'mean_word_length': 3.8730923694779116, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '264020', 'license': 'Apache-2.0', 'n_tokens_mistral': 4661, 'n_tokens_neox': 4321, 'n_words': 1160, 'path': 'Notebooks/Compare_Pressure_model.ipynb', 'repo_name': 'LockhartTech/DeeplearningProxy', 'size': '12133'}