diff --git a/src/bayesvalidrox/post_processing/post_processing.py b/src/bayesvalidrox/post_processing/post_processing.py
index 39ec81b75384ecc1a7038cef4c4679dc83172729..c32abd6267fd1038d64359ccf0bf1d439c11f39d 100644
--- a/src/bayesvalidrox/post_processing/post_processing.py
+++ b/src/bayesvalidrox/post_processing/post_processing.py
@@ -35,9 +35,6 @@ class PostProcessing:
 
     def __init__(self, engine, name='calib'):
         self.engine = engine
-        self.MetaModel = engine.MetaModel
-        self.ExpDesign = engine.ExpDesign
-        self.ModelObj = engine.Model
         self.name = name
 
     # -------------------------------------------------------------------------
@@ -63,20 +60,16 @@ class PostProcessing:
         """
 
         bar_plot = True if plot_type == 'bar' else False
-        meta_model_type = self.MetaModel.meta_model_type
-        Model = self.ModelObj
+        meta_model_type = self.engine.MetaModel.meta_model_type
 
         # Read Monte-Carlo reference
-        self.mc_reference = Model.read_observation('mc_ref')
+        self.mc_reference = self.engine.Model.read_observation('mc_ref')
 
         # Set the x values
         x_values_orig = self.engine.ExpDesign.x_values
 
         # Compute the moments with the PCEModel object
-        self.means, self.stds = self.MetaModel.calculate_moments()
-
-        # Get the variables
-        out_names = Model.Output.names
+        self.means, self.stds = self.engine.MetaModel.compute_moments()
 
         # Open a pdf for the plots
         newpath = (f'Outputs_PostProcessing_{self.name}/')
@@ -85,7 +78,7 @@ class PostProcessing:
 
         # Plot the best fit line, set the linewidth (lw), color and
         # transparency (alpha) of the line
-        for key in out_names:
+        for key in self.engine.out_names:
             fig, ax = plt.subplots(nrows=1, ncols=2)
 
             # Extract mean and std
@@ -173,9 +166,6 @@ class PostProcessing:
         None.
 
         """
-        MetaModel = self.MetaModel
-        Model = self.ModelObj
-
         if samples is None:
             self.n_samples = n_samples
             samples = self._get_sample()
@@ -189,12 +179,12 @@ class PostProcessing:
             self.model_out_dict = model_out_dict
         else:
             self.model_out_dict = self._eval_model(samples, key_str='valid')
-        self.pce_out_mean, self.pce_out_std = MetaModel.eval_metamodel(samples)
+        self.pce_out_mean, self.pce_out_std = self.engine.eval_metamodel(samples)
 
         try:
-            key = Model.Output.names[1]
+            key = self.engine.out_names[1]
         except IndexError:
-            key = Model.Output.names[0]
+            key = self.engine.out_names[0]
 
         n_obs = self.model_out_dict[key].shape[1]
 
@@ -232,9 +222,6 @@ class PostProcessing:
             Validation error for each output.
 
         """
-        MetaModel = self.MetaModel
-        Model = self.ModelObj
-
         # Set the number of samples
         if n_samples:
             self.n_samples = n_samples
@@ -252,12 +239,12 @@ class PostProcessing:
             outputs = self._eval_model(samples, key_str='validSet')
 
         # Run the PCE model with the generated samples
-        pce_outputs, _ = MetaModel.eval_metamodel(samples)
+        pce_outputs, _ = self.engine.eval_metamodel(samples)
 
         self.rmse = {}
         self.valid_error = {}
         # Loop over the keys and compute RMSE error.
-        for key in Model.Output.names:
+        for key in self.engine.out_names:
             # Root mena square
             self.rmse[key] = mean_squared_error(outputs[key], pce_outputs[key],
                                                 squared=False,
@@ -274,8 +261,8 @@ class PostProcessing:
                             in enumerate(zip(self.rmse[key],
                                              self.valid_error[key]))))
         # Save error dicts in PCEModel object
-        self.MetaModel.rmse = self.rmse
-        self.MetaModel.valid_error = self.valid_error
+        self.engine.MetaModel.rmse = self.rmse
+        self.engine.MetaModel.valid_error = self.valid_error
 
         return
 
@@ -296,7 +283,6 @@ class PostProcessing:
 
         """
         engine = self.engine
-        PCEModel = self.MetaModel
         n_init_samples = engine.ExpDesign.n_init_samples
         n_total_samples = engine.ExpDesign.X.shape[0]
 
@@ -564,13 +550,13 @@ class PostProcessing:
 
         """
         # This function currently only supports PCE/aPCE
-        if not hasattr(self.MetaModel, 'meta_model_type'):
+        PCEModel = self.engine.MetaModel
+        if not hasattr(PCEModel, 'meta_model_type'):
             raise AttributeError('Sobol indices currently only support PCE-type models!')
-        if self.MetaModel.meta_model_type.lower() not in ['pce', 'apce']:
+        if PCEModel.meta_model_type.lower() not in ['pce', 'apce']:
             raise AttributeError('Sobol indices currently only support PCE-type models!')
             
         # Extract the necessary variables
-        PCEModel = self.MetaModel
         basis_dict = PCEModel._basis_dict
         coeffs_dict = PCEModel._coeffs_dict
         n_params = PCEModel.ndim
@@ -583,7 +569,7 @@ class PostProcessing:
 
             sobol_cell_, total_sobol_ = {}, {}
 
-            for output in self.ModelObj.Output.names:
+            for output in self.engine.out_names:
 
                 n_meas_points = len(coeffs_dict[f'b_{b_i+1}'][output])
 
@@ -771,7 +757,7 @@ class PostProcessing:
                 total_sobol_all[k][i] = v
 
         self.total_sobol = {}
-        for output in self.ModelObj.Output.names:
+        for output in self.engine.out_names:
             self.total_sobol[output] = np.mean(total_sobol_all[output], axis=0)
 
         # ---------------- Plot -----------------------
@@ -784,7 +770,7 @@ class PostProcessing:
 
         fig = plt.figure()
 
-        for outIdx, output in enumerate(self.ModelObj.Output.names):
+        for outIdx, output in enumerate(self.engine.out_names):
 
             # Extract total Sobol indices
             total_sobol = self.total_sobol[output]
@@ -859,8 +845,6 @@ class PostProcessing:
         None.
 
         """
-        MetaModel = self.MetaModel
-
         if samples is None:
             self.n_samples = n_samples
             samples = self._get_sample()
@@ -869,7 +853,7 @@ class PostProcessing:
 
         # Evaluate the original and the surrogate model
         y_val = self._eval_model(samples, key_str='valid')
-        y_pce_val, _ = MetaModel.eval_metamodel(samples=samples)
+        y_pce_val, _ = self.engine.eval_metamodel(samples=samples)
 
         # Open a pdf for the plots
         newpath = f'Outputs_PostProcessing_{self.name}/'
@@ -968,8 +952,7 @@ class PostProcessing:
 
         self.n_samples = 1000
 
-        PCEModel = self.MetaModel
-        Model = self.ModelObj
+        PCEModel = self.engine.MetaModel
         n_samples = self.n_samples
 
         # Create 3D-Grid
@@ -1006,7 +989,7 @@ class PostProcessing:
                     pce_out_std[:, idx] = y_std
 
                     # Model evaluation
-                    model_out_dict, _ = Model.run_model_parallel(sample_mesh,
+                    model_out_dict, _ = self.engine.Model.run_model_parallel(sample_mesh,
                                                                  key_str='Valid3D')
                     model_outs[:, idx] = model_out_dict[Outkey].T
 
@@ -1066,7 +1049,7 @@ class PostProcessing:
         """
         if n_samples is None:
             n_samples = self.n_samples
-        self.samples = self.ExpDesign.generate_samples(
+        self.samples = self.engine.ExpDesign.generate_samples(
             n_samples,
             sampling_method='random')
         return self.samples
@@ -1090,15 +1073,13 @@ class PostProcessing:
             Dictionary of results.
 
         """
-        Model = self.ModelObj
-
         if samples is None:
             samples = self._get_sample()
             self.samples = samples
         else:
             self.n_samples = len(samples)
 
-        model_outs, _ = Model.run_model_parallel(samples, key_str=key_str)
+        model_outs, _ = self.engine.Model.run_model_parallel(samples, key_str=key_str)
 
         return model_outs
 
@@ -1189,8 +1170,6 @@ class PostProcessing:
         None.
 
         """
-        Model = self.ModelObj
-
         newpath = f'Outputs_PostProcessing_{self.name}/'
         if not os.path.exists(newpath):
             os.makedirs(newpath)
@@ -1201,7 +1180,7 @@ class PostProcessing:
 
         fig = plt.figure()
         # Plot the model vs PCE model
-        for keyIdx, key in enumerate(Model.Output.names):
+        for keyIdx, key in enumerate(self.engine.out_names):
 
             y_pce_val = self.pce_out_mean[key]
             y_pce_val_std = self.pce_out_std[key]
@@ -1247,4 +1226,4 @@ class PostProcessing:
             plt.close()
 
         # Zip the subdirectories
-        Model.zip_subdirs(f'{Model.name}valid', f'{Model.name}valid_')
+        self.engine.Model.zip_subdirs(f'{self.engine.Model.name}valid', f'{self.engine.Model.name}valid_')