diff --git a/kalman/graphics/figure2.py b/kalman/graphics/figure2.py index f9b03f1..8b71483 100644 --- a/kalman/graphics/figure2.py +++ b/kalman/graphics/figure2.py @@ -30,7 +30,7 @@ def load_data(file): return dat -def plot_parameters(dat, deparameterize=False, ref=None): +def plot_parameters(dat, input_file, deparameterize=False, ref=None): ''' Plot the parameters in separate subplots with uncertainties. Args: @@ -42,21 +42,29 @@ def plot_parameters(dat, deparameterize=False, ref=None): if is_ipython(): plt.ion() + idx_a = input_file.find('/') + idx_b = input_file[idx_a+1::].find('/') + name_file = input_file[idx_a+1:idx_b+idx_a+1] + inputfile_path = 'results/' + name_file + '/input.yaml' + - inputfile_path = 'results/aorta_C/input.yaml' with open(inputfile_path) as file: inputfile = yaml.full_load(file) #true_val = [10,250,250,250,30] - true_val = [94,250,683,615,30] + true_val = [10,60,220,160,100] + 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']) + for bnd_c in inputfile['estimation']['boundary_conditions']: + if 'windkessel' in bnd_c['type']: + for bnd_set in inputfile['boundary_conditions']: + if bnd_c['id'] == bnd_set['id']: + current_val.append(bnd_set['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)) @@ -72,6 +80,8 @@ def plot_parameters(dat, deparameterize=False, ref=None): col = cycle(['C0', 'C1', 'C2', 'C3','C4']) ls = cycle(['-', '-', '--', '--', ':', ':', '-.', '-.']) legends = cycle(['$R_3$','$R_4$','$R_5$','$R_6$','$U$']) + #legends = cycle(['$R_3$','$U$']) + col_ = next(col) ls_ = next(ls) @@ -100,7 +110,8 @@ def plot_parameters(dat, deparameterize=False, ref=None): # 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))) + print('Real values: \t {}'.format(true_val)) + print('Recon values: \t {}'.format(np.round(2**theta[-1, :]*current_val,2))) @@ -133,5 +144,5 @@ if __name__ == '__main__': args = get_parser().parse_args() dat = load_data(args.file) - - plot_parameters(dat, deparameterize=args.deparameterize, ref=args.ref) + + plot_parameters(dat, args.file,deparameterize=args.deparameterize, ref=args.ref) diff --git a/kalman/input_files/aorta.yaml b/kalman/input_files/aorta.yaml index 84ad1dc..8181195 100755 --- a/kalman/input_files/aorta.yaml +++ b/kalman/input_files/aorta.yaml @@ -6,7 +6,7 @@ fluid: dynamic_viscosity: 0.035 io: - write_path: 'results/aorta_C' + write_path: 'results/aorta_master' restart: path: '' # './projects/nse_coa3d/results/test_restart2/' time: 0 @@ -31,32 +31,32 @@ boundary_conditions: - id: 3 type: 'windkessel' - #value: [10,0,0] - #p0: [0,1333.223874] - value: [10,1000,0.01] - p0: [47,1333.223874] + value: [10,0,0] + p0: [0,1333.223874] + #value: [10,1000,0.01] + #p0: [47,1333.223874] - id: 4 type: 'windkessel' - #value: [250,0,0] - #p0: [0,1333.223874] - value: [250,8000,0.0001] - p0: [47,1333.223874] + value: [250,0,0] + p0: [0,1333.223874] + #value: [250,8000,0.0001] + #p0: [47,1333.223874] - id: 5 type: 'windkessel' - #value: [250,0,0] - #p0: [0,1333.223874] - value: [250,8000,0.0001] - p0: [47,1333.223874] + value: [250,0,0] + p0: [0,1333.223874] + #value: [250,8000,0.0001] + #p0: [47,1333.223874] - id: 6 type: 'windkessel' - #value: [250,0,0] - #p0: [0,1333.223874] - value: [250,8000,0.0001] - p0: [47,1333.223874] + value: [250,0,0] + p0: [0,1333.223874] + #value: [250,8000,0.0001] + #p0: [47,1333.223874] timemarching: velocity_pressure_coupling: 'fractionalstep' # monolithic, fractionalstep @@ -82,7 +82,7 @@ timemarching: flux_report_normalize_boundary: 1 T: 0.8 # end time - dt: 0.01 + dt: 0.001 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 @@ -120,22 +120,22 @@ linear_solver: estimation: boundary_conditions: - #- - # id: 3 - # type: 'windkessel' - # initial_stddev: 1 - #- - # id: 4 - # type: 'windkessel' - # initial_stddev: 1 - #- - # id: 5 - # type: 'windkessel' - # initial_stddev: 1 - #- - # id: 6 - # type: 'windkessel' - # initial_stddev: 1 + - + id: 3 + type: 'windkessel' + initial_stddev: 1 + - + id: 4 + type: 'windkessel' + initial_stddev: 1 + - + id: 5 + type: 'windkessel' + initial_stddev: 1 + - + id: 6 + type: 'windkessel' + initial_stddev: 1 - id: 2 type: 'dirichlet' @@ -147,13 +147,11 @@ estimation: - mesh: '/home/yeye/NuMRI/kalman/meshes/coaortaH3_leo2.0.h5' fe_degree: 1 - #xdmf_file: 'measurements/aorta_C/Perturbation/Ks12V120/u.xdmf' - #file_root: 'measurements/aorta_C/Perturbation/Ks12V120/u{i}.h5' - xdmf_file: 'measurements/aorta_C/u_all.xdmf' - file_root: 'measurements/aorta_C/u{i}.h5' + xdmf_file: 'measurements/aorta_master_s100/u_all.xdmf' + file_root: 'measurements/aorta_master_s100/u{i}.h5' indices: 0 # indices of checkpoints to be processed. 0 == all velocity_direction: ~ - noise_stddev: 15 # standard deviation of Gaussian noise + noise_stddev: 5 # standard deviation of Gaussian noise roukf: particles: 'simplex' # unique or simplex diff --git a/kalman/input_files/aorta_C.yaml b/kalman/input_files/aorta_C.yaml index 2e397f1..71b5913 100755 --- a/kalman/input_files/aorta_C.yaml +++ b/kalman/input_files/aorta_C.yaml @@ -6,14 +6,14 @@ fluid: dynamic_viscosity: 0.035 io: - write_path: 'results/aorta_C0_ten' + write_path: 'results/aorta_C' restart: path: '' # './projects/nse_coa3d/results/test_restart2/' time: 0 write_xdmf: True write_checkpoints: True write_hdf5_timeseries: False - write_velocity: 'tentative' # tentative or update + write_velocity: 'update' # tentative or update boundary_conditions: - @@ -24,41 +24,47 @@ boundary_conditions: id: 2 type: 'dirichlet' value: ['0','0','-U*sin(DOLFIN_PI*t/Th)*(t<=Th) + (t<=Tc)*(t>Th)*(-3.67949466208*U*sin(9*DOLFIN_PI*t/Th)*exp(-t*10)) - - U*sin(DOLFIN_PI*(t-Tc)/Th)*(t>Tc)*(t<= Tc + Th) + (t>Tc+Th)*(-3.67949466208*U*sin(9*DOLFIN_PI*(t-Tc)/Th)*exp(-(t-Tc)*10))'] + - U*sin(DOLFIN_PI*(t-Tc)/Th)*(t>Tc)*(t<= Tc + Th) + + (t<=2*Tc)*(t>Tc+Th)*(-3.67949466208*U*sin(9*DOLFIN_PI*(t-Tc)/Th)*exp(-(t-Tc)*10)) + - U*sin(DOLFIN_PI*(t-2*Tc)/Th)*(t>2*Tc)*(t<= 2*Tc + Th) + + (t<=3*Tc)*(t>2*Tc+Th)*(-3.67949466208*U*sin(9*DOLFIN_PI*(t-2*Tc)/Th)*exp(-(t-2*Tc)*10)) + - U*sin(DOLFIN_PI*(t-3*Tc)/Th)*(t>3*Tc)*(t<= 3*Tc + Th) + + (t<=4*Tc)*(t>3*Tc+Th)*(-3.67949466208*U*sin(9*DOLFIN_PI*(t-3*Tc)/Th)*exp(-(t-3*Tc)*10)) + - U*sin(DOLFIN_PI*(t-4*Tc)/Th)*(t>4*Tc)*(t<= 4*Tc + Th) + + (t<=5*Tc)*(t>4*Tc+Th)*(-3.67949466208*U*sin(9*DOLFIN_PI*(t-4*Tc)/Th)*exp(-(t-4*Tc)*10)) '] parameters: - U: 30 + U: 100 #100 safe Th: 0.35 Tc: 0.8 t: 0 - id: 3 type: 'windkessel' - #value: [10,1000,0.01] #paper corr - #value: [94,1794,0.0014] #paper windk - value: [94,0,0] - p0: [0,1333.223874] + value: [10,0.0008,2400] + #value: [10,0,2400] + p0: [80,1333.223874] + #p0: [0,1333.223874] - id: 4 type: 'windkessel' - #value: [250,8000,0.0001] #paper corr - #value: [250,10000,0.0004] #paper windk - value: [250,0,0] #paper - p0: [0,1333.223874] + value: [60,0.00034,4200] + #value: [60,0,4200] #safe + p0: [80,1333.223874] + #p0: [0,1333.223874] - id: 5 type: 'windkessel' - #value: [250,8000,0.0001] #paper corr - #value: [683,12960,0.0002] #paper windk - value: [683,0,0] - p0: [0,1333.223874] + value: [220,0.00034,11000] + #value: [220,0,11000] #safe + p0: [80,1333.223874] + #p0: [0,1333.223874] - id: 6 type: 'windkessel' - #value: [250,8000,0.0001] #paper corr - #value: [615,11664,0.0002] #paper windk - value: [615,0,0] - p0: [0,1333.223874] - + value: [160,0.00034,7800] + #value: [160,0,7800] #safe + p0: [80,1333.223874] + #p0: [0,1333.223874] timemarching: velocity_pressure_coupling: 'fractionalstep' # monolithic, fractionalstep @@ -83,8 +89,8 @@ timemarching: transpiration_bc_projection: 'robin' # robin, dirichlet flux_report_normalize_boundary: 1 - T: 1.6 # end time - dt: 0.001 + T: 0.8 # end time + dt: 0.002 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 @@ -129,7 +135,7 @@ estimation: #- # id: 4 # type: 'windkessel' - # initial_stddev: 1 + # initial_stddev: 1 #- # id: 5 # type: 'windkessel' @@ -148,8 +154,10 @@ estimation: - mesh: './meshes/coaortaH3_leo2.0.h5' fe_degree: 1 - xdmf_file: 'measurements/aorta_C/u_all.xdmf' - file_root: 'measurements/aorta_C/u{i}.h5' + #xdmf_file: 'measurements/aorta_C/Perturbation/Ks12V50/u_all.xdmf' + #file_root: 'measurements/aorta_C/Perturbation/Ks12V50/u{i}.h5' + xdmf_file: 'measurements/aorta_C_s100/u_all.xdmf' + file_root: 'measurements/aorta_C_s100/u{i}.h5' indices: 0 # indices of checkpoints to be processed. 0 == all velocity_direction: ~ noise_stddev: 15 # standard deviation of Gaussian noise @@ -160,4 +168,4 @@ estimation: reparameterize: True ODV_functional: enable: False - VENC: 102 # 102,120% 59,70% 42 50%, 21,25% + VENC: 172 # 241,172