changes

kalman/compute_errors_kalman.py

@@ -0,0 +1,51 @@
+import numpy as np
+
+vref = [10, 250, 250, 250, 30]
+#case 2
+v_no = [9.37, 246.18, 251.68, 252.37, 29.79]
+v_70 = [18.21, 217.19, 198, 227.76, 29]
+vnew_70 = [9.42, 244.59, 247.93, 246.89, 29.79]
+v_50 = [64.05,  205.9, 195.55, 246.47, 22.58]
+vnew_50 =  [  9.46 ,244.47 ,248.09 ,244.69  ,29.76]
+
+#acse 1
+v_no_case1 = [  9.58 , 246.97,  252.41 , 252.47 , 29.81]
+v_70_case1 = [17.93, 215.25, 193.09, 226.87, 29.02]
+v_50_case1 = [64.59, 207.49, 193.08, 248.51, 22.63]
+vnew_70_case1 =  [  9.51 , 237.32, 238.26, 234.33,  29.79]
+vnew_50_case1 = [  9.48, 237.58, 239.46, 231.55,  29.76]
+
+
+def compute_relative_error(vcur,vref):
+    eps = 0
+    for i in range(len(vref)):
+        eps += ((vcur[i] - vref[i])/vref[i])**2
+    
+    return np.sqrt(eps)
+
+
+eps_ref = compute_relative_error(v_no,vref)
+eps_70 = compute_relative_error(v_70,vref)
+eps_new70 = compute_relative_error(vnew_70,vref)
+eps_50 = compute_relative_error(v_50,vref)
+eps_new50 = compute_relative_error(vnew_50,vref)
+
+eps_ref_case1 = compute_relative_error(v_no_case1,vref)
+eps_70_case1 = compute_relative_error(v_70_case1,vref)
+eps_new70_case1 = compute_relative_error(vnew_70_case1,vref)
+eps_50_case1 = compute_relative_error(v_50_case1,vref)
+eps_new50_case1 = compute_relative_error(vnew_50_case1,vref)
+
+print('reference error = {e}'.format(e=np.round(100*eps_ref,2)))
+print('venc 70%  error = {e}'.format(e=np.round(100*eps_70,2)))
+print('odv venc 70% error = {e}'.format(e=np.round(100*eps_new70,2)))
+print('venc 50%  error = {e}'.format(e=np.round(100*eps_50,2)))
+print('odv venc 50% error = {e}'.format(e=np.round(100*eps_new50,2)))
+
+print('-------------------')
+
+print('reference error = {e}'.format(e=np.round(100*eps_ref_case1,2)))
+print('venc 70%  error = {e}'.format(e=np.round(100*eps_70_case1,2)))
+print('odv venc 70% error = {e}'.format(e=np.round(100*eps_new70_case1,2)))
+print('venc 50%  error = {e}'.format(e=np.round(100*eps_50_case1,2)))
+print('odv venc 50% error = {e}'.format(e=np.round(100*eps_new50_case1,2)))

kalman/graphics/figure2.py

@@ -43,7 +43,7 @@ def plot_parameters(dat, deparameterize=False, ref=None):
         plt.ion()
 
 
-    inputfile_path = 'results/aorta_leo/input.yaml'
+    inputfile_path = 'results/aorta/input.yaml'
     with open(inputfile_path) as file:
         inputfile = yaml.full_load(file)
 

kalman/input_files/aorta.yaml

@@ -134,18 +134,20 @@ estimation:
             initial_stddev: 1
 
 
-
     measurements:
         -
-            mesh: './meshes/coaortaH1.h5'
+            mesh: './meshes/coaortaH3_leo2.0.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: 5     # standard deviation of Gaussian noise
+            noise_stddev: 15     # standard deviation of Gaussian noise
     
     roukf:
         particles: 'simplex' # unique or simplex
         observation_operator: 'postprocessing'         #state or postprocessing
         reparameterize: True
+        ODV_functional:
+            enable: False
+            VENC: 44    # cm/s  ref: 61 for 70% and 44 for 50%

kalman/input_files/aorta_leo.yaml

@@ -30,22 +30,22 @@ boundary_conditions:
     - 
         id: 3
         type: 'windkessel'
-        value: [1,0,0]
+        value: [10,0,0]
         p0: [0,1333.223874]
     - 
         id: 4
         type: 'windkessel'
-        value: [100,0,0]
+        value: [250,0,0]
         p0: [0,1333.223874]
     - 
         id: 5
         type: 'windkessel'
-        value: [100,0,0]
+        value: [250,0,0]
         p0: [0,1333.223874]
     - 
         id: 6
         type: 'windkessel'
-        value: [100,0,0]
+        value: [250,0,0]
         p0: [0,1333.223874]
 
 timemarching:
@@ -138,8 +138,8 @@ estimation:
         -
             mesh: './meshes/coaortaH3_leo2.0.h5'
             fe_degree: 1
-            xdmf_file: 'results/aorta/measurements/u_all.xdmf'
+            xdmf_file: 'results/aorta_leo/measurements/u_all.xdmf'
-            file_root: 'results/aorta/measurements/u{i}.h5'
+            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