From ec760808394694eb241039800c171ea9f50cdc72 Mon Sep 17 00:00:00 2001
From: jeremias <jeremias.garay.l@gmail.com>
Date: Fri, 28 May 2021 10:01:29 +0200
Subject: [PATCH] clean up

---
 kalman/graphics/figure2.py            |  51 ++++++---
 kalman/input_files/aorta.yaml         | 117 +++++++++++---------
 kalman/input_files/aorta_leo.yaml     | 150 --------------------------
 kalman/input_files/aorta_neumann.yaml | 129 ----------------------
 4 files changed, 105 insertions(+), 342 deletions(-)
 delete mode 100755 kalman/input_files/aorta_leo.yaml
 delete mode 100755 kalman/input_files/aorta_neumann.yaml

diff --git a/kalman/graphics/figure2.py b/kalman/graphics/figure2.py
index de42690..90e55d5 100644
--- a/kalman/graphics/figure2.py
+++ b/kalman/graphics/figure2.py
@@ -52,15 +52,41 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
         inputfile = yaml.full_load(file)
 
 
-    true_values = {
-            3: 2400,
-            4: 4200,
-            5: 11000,
-            6: 7800,
-            2: 100
-            }
+    #true_values = {
+    #        3: 3400,
+    #        4: 4200,
+    #        5: 11000,
+    #        6: 7800,
+    #        2: 100
+    #        }
   
-    
+    true_values = {
+            3: 4800,
+            4: 7020,
+            5: 11520,
+            6: 11520,
+            2: 75
+            }
+
+
+
+    true_values_c = {
+        3: 0.0008,
+        4: 0.00034,
+        5: 0.00034,
+        6: 0.00034,
+        2: 100
+        }
+
+    true_values_rp = {
+        3: 10,
+        4: 60,
+        5: 220,
+        6: 160,
+        2: 100
+        }
+
+
     current_val = []
     labels = []
     ids = []
@@ -75,10 +101,9 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
                     labels.append('$R_' + str(bnd_c['id']))
 
         elif 'dirichlet' in bnd_c['type']:
-            current_val.append(inputfile['boundary_conditions'][1]['parameters']['U'])
+            current_val.append(inputfile['boundary_conditions'][0]['parameters']['U'])
             ids.append(bnd_c['id'])
-            labels.append('U')
-
+            labels.append('$U')
 
     dim = dat['theta'].shape[-1]
     fig1, axes = plt.subplots(1,1,figsize=(8,6))
@@ -107,10 +132,10 @@ def plot_parameters(dat, input_file, deparameterize=False, ref=None):
     for i in range(dim):
         
         true_level =  np.log(true_values[ids[i]]/current_val[i])/np.log(2)
-        rec_value = np.round(2**theta[-1, i]*current_val[i],1)
+        rec_value = np.round(2**theta[-1, i]*current_val[i],2)
         cur_key = ids[i]
 
-        axes.plot(t, theta[:, i] + 1.5*i, '-', color=col_,label=legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key])  + '$')
+        axes.plot(t, theta[:, i] + 1.5*i, '-', color=col_,label= legends_ + '= ' + str(rec_value) + '/' + str(true_values[cur_key])  + '$')
         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_)
diff --git a/kalman/input_files/aorta.yaml b/kalman/input_files/aorta.yaml
index e007e9c..9807f85 100755
--- a/kalman/input_files/aorta.yaml
+++ b/kalman/input_files/aorta.yaml
@@ -4,10 +4,11 @@ mesh: './meshes/coaortaH1.h5'
 fluid:
     density: 1.2
     dynamic_viscosity: 0.035
-    stokes: True
+    stokes: False
+    implicit_windkessel: True
 
 io:
-    write_path: 'results/aorta_1ms'
+    write_path: 'results/aorta25_zPb_3R'
     restart:
         path: ''  # './projects/nse_coa3d/results/test_restart2/' 
         time: 0
@@ -17,53 +18,67 @@ io:
     write_velocity: 'update'     # update or tentative
 
 boundary_conditions:
+    - 
+        id: 2
+        type: 'dirichlet'
+        #value: ['0','0','-U*sin(DOLFIN_PI*t/Th)*(t<=Th) + (t<0.8)*(Th<t)*(U*DOLFIN_PI/Th*(t-Th)*exp(-(t-Th)*beta)) +
+        #                 -U*sin(DOLFIN_PI*(t-0.8)/Th)*(t<= 0.8+Th )*(t>0.8) + (t<1.6)*(0.8+Th<t)*(U*DOLFIN_PI/Th*(t-0.8-Th)*exp(-(t-0.8-Th)*beta)) +
+        #                 -U*sin(DOLFIN_PI*(t-1.6)/Th)*(t<= 1.6+Th )*(t>1.6) + (t<2.4)*(1.6+Th<t)*(U*DOLFIN_PI/Th*(t-1.6-Th)*exp(-(t-1.6-Th)*beta))' ]
+        value: ['0','0','-U*sin(DOLFIN_PI*t/Th)*(t<=Th) + (Th<t)*(U*DOLFIN_PI/Th*(t-Th)*exp(-(t-Th)*beta))']
+        parameters:
+            #U: 75
+            #U: 150 #Pa
+            U: 40 #Pb
+            Th: 0.36
+            beta: 70
+            t: 0
     - 
         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: 100
-            Th: 0.35
-            t: 0   
     - 
         id: 3
         type: 'windkessel'
         parameters:
-            R_p: 10
-            C: 0.0008
-            R_d: 2400
-            p0: 80
+            R_p: 200
+            C: 0.0005
+            #R_d: 4800
+            #R_d: 9000 #Pa
+            R_d: 2000 #Pb
+            p0: 85
             conv: 1333.223874   
     - 
         id: 4
         type: 'windkessel'
         parameters:
-            R_p: 60
-            C: 0.00034
-            R_d: 4200
-            p0: 80
+            R_p: 480
+            C: 0.00045
+            #R_d: 7200
+            #R_d: 12000 #Pa
+            R_d: 4000 #Pb
+            p0: 85
             conv: 1333.223874   
     - 
         id: 5
         type: 'windkessel'
         parameters:
-            R_p: 220
-            C: 0.00034
-            R_d: 11000
-            p0: 80
+            R_p: 520
+            C: 0.00045
+            #R_d: 11520
+            #R_d: 23000 #Pa
+            R_d: 6000 #Pb
+            p0: 85
             conv: 1333.223874   
     - 
         id: 6
         type: 'windkessel'
         parameters:
-            R_p: 160
-            C: 0.00034
-            R_d: 7800
-            p0: 80
+            R_p: 520
+            C: 0.00045
+            R_d: 11520
+            #R_d: 23000 #Pa
+            #R_d: 6000 #Pb
+            p0: 85
             conv: 1333.223874   
 
 timemarching:
@@ -91,8 +106,8 @@ timemarching:
 
     T: 0.8       # end time
     dt: 0.001
-    write_dt: 0.04
-    checkpoint_dt:  0.04 # <= 0: only last; else value + last
+    write_dt: 0.03
+    checkpoint_dt:  0.03 # <= 0: only last; else value + last
     report: 1   # 0: print nothing, 1: print time step and writeout, 2: 1 + flux
 
 # solver setup
@@ -104,9 +119,9 @@ fem:
     convection_skew_symmetric: True    # aka Temam term
     stabilization:
         forced_normal:
-            enabled: False
-            boundaries: [6]
-            gamma: 10
+            enabled: True
+            boundaries: [3,4,5,6]
+            gamma: 20
         backflow_boundaries: [3,4,5,6]
         streamline_diffusion:
             enabled: False
@@ -128,39 +143,41 @@ linear_solver:
 
 estimation:
     boundary_conditions:
-        #- 
-        #    id: 3
-        #    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: 5
+            type: 'windkessel'
+            initial_stddev: 1
         #- 
         #    id: 6
         #    type: 'windkessel'
         #    initial_stddev: 1
-        #- 
-        #    id: 2
-        #    type: 'dirichlet'
-        #    parameters: 'U'
-        #    initial_stddev: 1
+        - 
+            id: 2
+            type: 'dirichlet'
+            parameters: 'U'
+            initial_stddev: 1
 
 
     measurements:
         -
-            #mesh: '/home/yeye/NuMRI/kalman/meshes/coaortaH3_leo2.0.h5'
-            mesh: './meshes/coaortaH1.h5'
+            mesh: '/home/yeye/NuMRI/kalman/meshes/coaortaH3_leo2.0.h5'
+            #mesh: './meshes/coaortaH1.h5'
             fe_degree: 1
-            xdmf_file: 'measurements/aorta/u_all.xdmf'
-            file_root: 'measurements/aorta/u{i}.h5'
+            xdmf_file: 'measurements/aorta25_z/Perturbation/Mg15V120/u_all.xdmf'
+            file_root: 'measurements/aorta25_z/Perturbation/Mg15V120/u{i}.h5'
+            #xdmf_file: 'measurements/aorta25_z/u_all.xdmf'
+            #file_root: 'measurements/aorta25_z/u{i}.h5'
             indices: 0 # indices of checkpoints to be processed. 0 == all
-            velocity_direction: ~
-            noise_stddev: 0     # standard deviation of Gaussian noise
+            velocity_direction: [0,0,1]
+            noise_stddev: 55     # standard deviation of Gaussian noise
     
     roukf:
         particles: 'simplex' # unique or simplex
diff --git a/kalman/input_files/aorta_leo.yaml b/kalman/input_files/aorta_leo.yaml
deleted file mode 100755
index 53fe255..0000000
--- a/kalman/input_files/aorta_leo.yaml
+++ /dev/null
@@ -1,150 +0,0 @@
-mesh: './meshes/coaortaH1.h5'
-# Physical parameters of the fluid
-fluid:
-    density: 1.2
-    dynamic_viscosity: 0.035
-
-io:
-    write_path: 'results/aorta_leo'
-    restart:
-        path: ''  # './projects/nse_coa3d/results/test_restart2/' 
-        time: 0
-    write_xdmf: True
-    write_checkpoints: True
-    write_hdf5_timeseries: False
-    write_velocity: 'update'     # update or 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,0,0]
-        p0: [0,1333.223874]
-    - 
-        id: 4
-        type: 'windkessel'
-        value: [250,0,0]
-        p0: [0,1333.223874]
-    - 
-        id: 5
-        type: 'windkessel'
-        value: [250,0,0]
-        p0: [0,1333.223874]
-    - 
-        id: 6
-        type: 'windkessel'
-        value: [250,0,0]
-        p0: [0,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: 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'
-            parameters: 'U'
-            initial_stddev: 1
-
-
-
-    measurements:
-        -
-            mesh: './meshes/coaortaH3_leo2.0.h5'
-            fe_degree: 1
-            xdmf_file: 'results/aorta_leo/measurements/u_all.xdmf'
-            file_root: 'results/aorta_leo/measurements/u{i}.h5'
-            indices: 0 # indices of checkpoints to be processed. 0 == all
-            velocity_direction: ~
-            noise_stddev: 5     # standard deviation of Gaussian noise
-    
-    roukf:
-        particles: 'simplex' # unique or simplex
-        observation_operator: 'postprocessing'         #state or postprocessing
-        reparameterize: True
diff --git a/kalman/input_files/aorta_neumann.yaml b/kalman/input_files/aorta_neumann.yaml
deleted file mode 100755
index d38df81..0000000
--- a/kalman/input_files/aorta_neumann.yaml
+++ /dev/null
@@ -1,129 +0,0 @@
-
-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: 'neumann'
-        value: 0
-    - 
-        id: 4
-        type: 'neumann'
-        value: 0
-    - 
-        id: 5
-        type: 'neumann'
-        value: 0
-    - 
-        id: 6
-        type: 'neumann'
-        value: 0
-
-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: 1.5     # standard deviation of Gaussian noise
-    
-    roukf:
-        particles: 'simplex' # unique or simplex
-        observation_operator: 'postprocessing'         #state or postprocessing
-        reparameterize: True