philippslang

void demo(size_t nr, size_t nc, double *arr)
{
  for(int r(0); r < nr; ++r)
    for(int c(0); c < nc; ++c)
      arr[r*nc+c] = 1.; // fills arr[r][c]
}

philippslang

function [hmame, hmami] = Ellipsoid_discs_plot(vma, vme, vmi)
%Ellipsoid_discs_plot plots an ellipsoid as two
%elliptical discs.
%
%The center of the ellipsoid is assumed to be 0,0,0. The tool has been
%designed to visualize tensorial properties using the characteristic
%Eigenvectors. Using cross-section elliptical discs to visualize ellipsoids
%may provide improved visualization compared to volumetric ellipsoids.
%
% Parameters

philippslang

ic_load_tetin ./out.tin

set srfcs [ic_geo_get_objects surface]
ic_curve intersect STANDARD crv.00 $srfcs
set crvs [ic_geo_get_objects curve]
ic_geo_create_curve_ends $crvs
ic_save_tetin flow.tin

ic_geo_scale_meshing_params all 0.9
ic_set_meshing_params global 0 gref 0.9 gmax 9.0 gfast 1 gedgec 0.2 gnat 0 gcgap 1 gnatref 10

philippslang

import numpy as np
import keras.models as kem
import keras.layers as kel
import sklearn.preprocessing as preproc
np.random.seed(42)

FNAME_MODEL = 'tmp/decl_curve_rnn.h5'

philippslang

def calc_exponential_decline(p0, exp, time):
    return p0 * np.exp(-exp * time)

def calc_two_stage_decline(p0, exp_stage_zero, exp_stage_one, time_max, time_next_stage, time_min=0., 
  production_jumpfactor=4., num=50, noise=0.1, noise_mean=1.):
    time = np.linspace(time_min, time_max, num)
    stage_one = np.where(time > time_next_stage)    
    production = calc_exponential_decline(p0, exp_stage_zero, time) * np.random.normal(noise_mean, noise, time.shape) 
    time_stage_one_relative = np.linspace(time_min, time_max-time_next_stage, len(time[stage_one]))
    production[stage_one] = calc_exponential_decline(production[stage_one][0] + p0/production_jumpfactor, exp_stage_one, 

philippslang

def make_rnn(num_features, num_targets, num_timesteps, num_units):
    model = kem.Sequential()
    model.add(kel.LSTM(num_units, input_shape=(num_timesteps, num_features), return_sequences=True))
    model.add(kel.TimeDistributed(kel.Dense(num_targets, activation='sigmoid')))
    model.compile(loss='mse', optimizer='adam')
    return model

philippslang

num_features = 2 
ifeature_production = 0
ifeature_stage = 1

num_targets = 1 
itarget_production = 0 

# lstm units in the recurrent layer
num_units = 24 

philippslang

p0 = 50. # production rate at time zero
na = -p0 # numeric encoding of not available value
exp_stage_zero = 0.12 # exponent of production decline for stage zero
exp_stage_one = 0.1
time_max = 55.
bounds_stage_change_time = (20., 40.) # we'll generate training data with stage changes between these

num_timesteps = 50  # so many production values per sequence
num_discrete_stage_changes = 5 # we'll generate profiles with this many different stage change times
num_realizations_per_stage_change = 10 # and for each of these times, this many realizations (random noise differs)

philippslang

isample = 0
for time_stage_change in np.linspace(*bounds_stage_change_time, num_discrete_stage_changes):
    for _ in range(num_realizations_per_stage_change):
        _, production, stage = calc_two_stage_decline(p0, exp_stage_zero, exp_stage_one, time_max, t
          time_stage_change, num=num_timesteps) 
        for num_sample_points in range(1, num_timesteps):
            features[isample, :, ifeature_stage] = stage[:] 
            features[isample, :num_sample_points, ifeature_production] = production[:num_sample_points]
            targets[isample, 1:, itarget_production] = production[1:] 
            isample += 1

philippslang

Scaler = preproc.MinMaxScaler
normalizer_features = Scaler(copy=False)
normalizer_features.fit_transform(features.reshape(num_sequences*num_timesteps, num_features))
normalizer_targets = Scaler(feature_range=(0, 1), copy=False) # sigmoid
normalizer_targets.fit_transform(targets.reshape(num_sequences*num_timesteps, num_targets))