modifications

This commit is contained in:
jeremias 2020-12-16 14:41:44 -03:00
parent ae04e3d0cc
commit 013d7e87ba
5 changed files with 414 additions and 16 deletions

View File

@ -6,14 +6,14 @@ fluid:
dynamic_viscosity: 0.035 dynamic_viscosity: 0.035
io: io:
write_path: 'results/aorta/' write_path: 'results/aorta'
restart: restart:
path: '' # './projects/nse_coa3d/results/test_restart2/' path: '' # './projects/nse_coa3d/results/test_restart2/'
time: 0 time: 0
write_xdmf: True write_xdmf: True
write_checkpoints: True write_checkpoints: True
write_hdf5_timeseries: False write_hdf5_timeseries: False
write_velocity: 'update' # tentative write_velocity: 'update' # update or tentative
boundary_conditions: boundary_conditions:
- -
@ -25,29 +25,29 @@ boundary_conditions:
type: 'dirichlet' type: 'dirichlet'
value: ['0','0','-U*sin(DOLFIN_PI*t/Th)*(t<=Th) + (Th<t)*(-3.67949466208*U*sin(9*DOLFIN_PI*t/Th)*exp(-t*10))'] value: ['0','0','-U*sin(DOLFIN_PI*t/Th)*(t<=Th) + (Th<t)*(-3.67949466208*U*sin(9*DOLFIN_PI*t/Th)*exp(-t*10))']
parameters: parameters:
U: 30 U: 60
Th: 0.35 Th: 0.35
t: 0 t: 0
- -
id: 3 id: 3
type: 'windkessel' type: 'windkessel'
value: [10,0.01,1000] value: [20,0,0]
p0: [47,1333.223874] p0: [0,1333.223874]
- -
id: 4 id: 4
type: 'windkessel' type: 'windkessel'
value: [250,0.0001,8000] value: [500,0,0]
p0: [47,1333.223874] p0: [0,1333.223874]
- -
id: 5 id: 5
type: 'windkessel' type: 'windkessel'
value: [250,0.0001,8000] value: [500,0,0]
p0: [47,1333.223874] p0: [0,1333.223874]
- -
id: 6 id: 6
type: 'windkessel' type: 'windkessel'
value: [250,0.0001,8000] value: [500,0,0]
p0: [47,1333.223874] p0: [0,1333.223874]
timemarching: timemarching:
velocity_pressure_coupling: 'fractionalstep' # monolithic, fractionalstep velocity_pressure_coupling: 'fractionalstep' # monolithic, fractionalstep
@ -111,12 +111,30 @@ linear_solver:
estimation: estimation:
boundary_conditions: boundary_conditions:
-
id: 3
type: 'windkessel'
initial_stddev: 1
- -
id: 4 id: 4
type: 'windkessel' type: 'windkessel'
initial_stddev: 1
-
id: 5
type: 'windkessel'
initial_stddev: 1
-
id: 6
type: 'windkessel'
initial_stddev: 1
-
id: 2
type: 'dirichlet'
parameters: 'U' parameters: 'U'
initial_stddev: 1 initial_stddev: 1
measurements: measurements:
- -
mesh: './meshes/coaortaH1.h5' mesh: './meshes/coaortaH1.h5'
@ -125,7 +143,7 @@ estimation:
file_root: 'results/aorta/measurements/u{i}.h5' file_root: 'results/aorta/measurements/u{i}.h5'
indices: 0 # indices of checkpoints to be processed. 0 == all indices: 0 # indices of checkpoints to be processed. 0 == all
velocity_direction: ~ velocity_direction: ~
noise_stddev: 1.5 # standard deviation of Gaussian noise noise_stddev: 25 # standard deviation of Gaussian noise
roukf: roukf:
particles: 'simplex' # unique or simplex particles: 'simplex' # unique or simplex

133
kalman/aorta_C.yaml Executable file
View File

@ -0,0 +1,133 @@
mesh: './meshes/coaortaH1.h5'
# Physical parameters of the fluid
fluid:
density: 1.2
dynamic_viscosity: 0.035
io:
write_path: 'results/aorta/'
restart:
path: '' # './projects/nse_coa3d/results/test_restart2/'
time: 0
write_xdmf: True
write_checkpoints: True
write_hdf5_timeseries: False
write_velocity: 'update' # tentative
boundary_conditions:
-
id: 1
type: 'dirichlet'
value: ['0','0','0']
-
id: 2
type: 'dirichlet'
value: ['0','0','-U*sin(DOLFIN_PI*t/Th)*(t<=Th) + (Th<t)*(-3.67949466208*U*sin(9*DOLFIN_PI*t/Th)*exp(-t*10))']
parameters:
U: 30
Th: 0.35
t: 0
-
id: 3
type: 'windkessel'
value: [10,1000,0.01]
p0: [47,1333.223874]
-
id: 4
type: 'windkessel'
value: [250,8000,0.0001]
p0: [47,1333.223874]
-
id: 5
type: 'windkessel'
value: [250,8000,0.0001]
p0: [47,1333.223874]
-
id: 6
type: 'windkessel'
value: [250,8000,0.0001]
p0: [47,1333.223874]
timemarching:
velocity_pressure_coupling: 'fractionalstep' # monolithic, fractionalstep
monolithic:
timescheme: 'gmp' # generalized midpoint, steady FIXME TODO
theta: 1 # 1: Euler, 0.5: implicit midpoint rule (one-legged)
nonlinear:
method: 'constant_extrapolation' # constant_extrapolation, linear_extrapolation, newton, picard, snes
maxit: 20
init_steps: 30
use_aitken: 1 # 0: False, 1: Picard only, 2: all
report: 1 # 0: None, 1: residuals, 2: residuals and energy (inflow/driving/forcing via ESSENTIAL Dbcs!)
atol: 1.e-6 # note: dot required!!
rtol: 1.e-16
stol: 0.0
fractionalstep:
scheme: 'CT' # CT, IPCS
coupled_velocity: False # False faster, True needed if robin_bc implicit
robin_bc_velocity_scheme: 'implicit' # explicit, semi-implicit, implicit
transpiration_bc_projection: 'robin' # robin, dirichlet
flux_report_normalize_boundary: 1
T: 0.8 # end time
dt: 0.01
write_dt: 0.04
checkpoint_dt: 0.04 # <= 0: only last; else value + last
report: 1 # 0: print nothing, 1: print time step and writeout, 2: 1 + flux
# solver setup
fem:
velocity_space: p1 # p1 p1b/p1+ p2
pressure_space: p1 # p1 p0/dg0 dg1
strain_symmetric: False
convection_skew_symmetric: True # aka Temam term
stabilization:
forced_normal:
enabled: True
boundaries: [6]
gamma: 10
backflow_boundaries: [3,4,5,6]
streamline_diffusion:
enabled: False
parameter: 'standard' # standard, shakib, codina, klr
length_scale: 'metric' # average, max, metric
parameter_element_constant: True
Cinv: ~
monolithic:
infsup: 'pspg' # pspg, pressure-stabilization
graddiv: False
consistent: False
pressure_stab_constant: 1.
fix_pressure: False
fix_pressure_point: [0., 0. , 0.]
linear_solver:
method: 'lu'
estimation:
boundary_conditions:
-
id: 2
type: 'dirichlet'
parameters: 'U'
initial_stddev: 1
measurements:
-
mesh: './meshes/coaortaH1.h5'
fe_degree: 1
xdmf_file: 'results/aorta/measurements/u_all.xdmf'
file_root: 'results/aorta/measurements/u{i}.h5'
indices: 0 # indices of checkpoints to be processed. 0 == all
velocity_direction: ~
noise_stddev: 12.5 # standard deviation of Gaussian noise
roukf:
particles: 'simplex' # unique or simplex
observation_operator: 'postprocessing' #state or postprocessing
reparameterize: True

View File

@ -5,12 +5,12 @@ fluid:
dynamic_viscosity: 0.035 dynamic_viscosity: 0.035
io: io:
write_path: 'results/channel3d/' write_path: 'results/channel3d/60/'
restart: restart:
path: '' # './projects/nse_coa3d/results/test_restart2/' path: '' # './projects/nse_coa3d/results/test_restart2/'
time: 0 time: 0
write_xdmf: True write_xdmf: True
write_checkpoints: True write_checkpoints: False
write_hdf5_timeseries: False write_hdf5_timeseries: False
write_velocity: 'update' # update or tentative write_velocity: 'update' # update or tentative
@ -26,7 +26,7 @@ boundary_conditions:
type: 'dirichlet' type: 'dirichlet'
value: ['0','0','U*(1-x[0]*x[0] - x[1]*x[1])*sin(DOLFIN_PI*t/T)'] value: ['0','0','U*(1-x[0]*x[0] - x[1]*x[1])*sin(DOLFIN_PI*t/T)']
parameters: parameters:
U: 30 U: 60
T: 0.9 T: 0.9
t: 0 t: 0
- -
@ -111,7 +111,7 @@ estimation:
file_root: 'results/channel3d/measurements/u{i}.h5' file_root: 'results/channel3d/measurements/u{i}.h5'
indices: 0 # indices of checkpoints to be processed. 0 == all indices: 0 # indices of checkpoints to be processed. 0 == all
velocity_direction: ~ velocity_direction: ~
noise_stddev: 1.0 # standard deviation of Gaussian noise noise_stddev: 10.0 # standard deviation of Gaussian noise
roukf: roukf:
particles: 'simplex' # unique or simplex particles: 'simplex' # unique or simplex

112
kalman/graphics/figure1.py Normal file
View File

@ -0,0 +1,112 @@
import matplotlib.pyplot as plt
import numpy as np
from itertools import cycle
import argparse
import pickle
def is_ipython():
''' Check if script is run in IPython.
Returns:
bool: True if IPython, else False '''
try:
get_ipython()
ipy = True
except NameError:
ipy = False
return ipy
def load_data(file):
''' Load numpy data from file.
Returns
dict: data dictionary
'''
dat = np.load(file)
return dat
def plot_parameters(dat):
''' Plot the parameters in separate subplots with uncertainties.
Args:
dat (dict): data dictionary
deparameterize (bool): flag indicating if parameters should be
deparameterized via 2**theta
ref: reference value to be plotted with parameters
'''
if is_ipython():
plt.ion()
fig1, axes = plt.subplots(1, 1,figsize=(8,6))
axes.set_ylabel(r'$\theta$',fontsize=18)
col = cycle(['C0', 'C1', 'C2', 'C3'])
for k in dat.keys():
t = dat[k]['times']
theta = dat[k]['theta']
P = dat[k]['P_theta']
theta = theta.reshape((-1, 1))
P = P.reshape((-1, 1, 1))
#theta = 2**theta*float(k)
sP = np.sqrt(P[:,0,0])
sP_up = 2**sP
sP_down = 2**(-sP)
col_ = next(col)
axes.plot(t, theta[:, 0], '-', c=col_)
#axes.fill_between(t, theta[:, 0]*sP_down[:], theta[:, 0]*sP_up[:], alpha=0.3,color=col_)
axes.fill_between(t, theta[:, 0] - np.sqrt(P[:, 0, 0]),theta[:, 0] + np.sqrt(P[:, 0, 0]), alpha=0.3,color=col_)
axes.set_xlabel(r'time',fontsize=18)
axes.plot(t, t*0 + np.log(30/float(k))/np.log(2), '-', c='black',ls='--')
axes.set_xlim([-0.01,0.2])
if not is_ipython():
plt.show()
def get_parser():
parser = argparse.ArgumentParser(
description='''
Plot the time evolution of the ROUKF estimated parameters.
To execute in IPython::
%run plot_roukf_parameters.py [-d] [-r N [N \
...]] file
''',
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('file', type=str, help='path to ROUKF stats file')
parser.add_argument('-d', '--deparameterize', action='store_true',
help='deparameterize the parameters by 2**theta')
parser.add_argument('-r', '--ref', metavar='N', nargs='+', default=None,
type=float, help='Reference values for parameters')
return parser
if __name__ == '__main__':
args = get_parser().parse_args()
files = ['10','30','60']
dat_array = {}
for ff in files:
path = args.file + ff + '/theta_stats.npz'
dat_array[ff] = load_data(path)
plot_parameters(dat_array)

135
kalman/graphics/figure2.py Normal file
View File

@ -0,0 +1,135 @@
import matplotlib.pyplot as plt
import numpy as np
from itertools import cycle
import argparse
import pickle
import yaml
def is_ipython():
''' Check if script is run in IPython.
Returns:
bool: True if IPython, else False '''
try:
get_ipython()
ipy = True
except NameError:
ipy = False
return ipy
def load_data(file):
''' Load numpy data from file.
Returns
dict: data dictionary
'''
dat = np.load(file)
return dat
def plot_parameters(dat, deparameterize=False, ref=None):
''' Plot the parameters in separate subplots with uncertainties.
Args:
dat (dict): data dictionary
deparameterize (bool): flag indicating if parameters should be
deparameterized via 2**theta
ref: reference value to be plotted with parameters
'''
if is_ipython():
plt.ion()
inputfile_path = 'results/aorta/input.yaml'
with open(inputfile_path) as file:
inputfile = yaml.full_load(file)
true_val = [10,250,250,250,30]
current_val = []
current_val.append(inputfile['boundary_conditions'][2]['value'][0])
current_val.append(inputfile['boundary_conditions'][3]['value'][0])
current_val.append(inputfile['boundary_conditions'][4]['value'][0])
current_val.append(inputfile['boundary_conditions'][5]['value'][0])
current_val.append(inputfile['boundary_conditions'][1]['parameters']['U'])
dim = dat['theta'].shape[-1]
fig1, axes = plt.subplots(1,1,figsize=(8,6))
axes.set_ylabel(r'$\theta$',fontsize=18)
t = dat['times']
theta = dat['theta']
P = dat['P_theta']
col = cycle(['C0', 'C1', 'C2', 'C3','C4'])
ls = cycle(['-', '-', '--', '--', ':', ':', '-.', '-.'])
legends = cycle(['$R_3$','$R_4$','$R_5$','$R_6$','$U$'])
col_ = next(col)
ls_ = next(ls)
legends_=next(legends)
if dim == 1:
theta = theta.reshape((-1, 1))
P = P.reshape((-1, 1, 1))
for i in range(dim):
axes.plot(t, theta[:, i] + 1.5*i, '-', color=col_,label=legends_)
axes.fill_between(t, theta[:, i] + 1.5*i - np.sqrt(P[:, i, i]),
theta[:, i] + 1.5*i + np.sqrt(P[:, i, i]), alpha=0.3,
color=col_)
true_level = np.log(true_val[i]/current_val[i])/np.log(2)
axes.plot(t,1.5*i + t*0 + true_level , color=col_,ls='--')
col_ = next(col)
legends_=next(legends)
axes.legend(fontsize=14,loc='lower right')
axes.set_xlim([-0.01,0.81])
axes.set_xlabel(r'time (s)',fontsize=18)
# print('theta_peak: \t {}'.format(theta[round(len(theta)/2), :]))
print('Final value theta: \t {}'.format(theta[-1, :]))
print('Deparameterized: 2^theta_end: \t {}'.format(2**theta[-1, :]))
print('Real values: \t {}'.format(np.round(2**theta[-1, :]*current_val,2)))
plt.savefig('windk_res')
if not is_ipython():
plt.show()
def get_parser():
parser = argparse.ArgumentParser(
description='''
Plot the time evolution of the ROUKF estimated parameters.
To execute in IPython::
%run plot_roukf_parameters.py [-d] [-r N [N \
...]] file
''',
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('file', type=str, help='path to ROUKF stats file')
parser.add_argument('-d', '--deparameterize', action='store_true',
help='deparameterize the parameters by 2**theta')
parser.add_argument('-r', '--ref', metavar='N', nargs='+', default=None,
type=float, help='Reference values for parameters')
return parser
if __name__ == '__main__':
args = get_parser().parse_args()
dat = load_data(args.file)
plot_parameters(dat, deparameterize=args.deparameterize, ref=args.ref)