[Data Assimilation] Getting ready for DA cycles

.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H

@@ -129,20 +129,22 @@ read_plot_file(amrex::PlotFileData& pf,
                amrex::MultiFab& mf);
 
 amrex::MultiFab
-read_member_multifab(int n,
+read_member_multifab(const int da_iter,
+                     const int n,
                      const std::string& pf_name,
                      const amrex::Vector<std::string>& varnames);
 
 amrex::MultiFab
-compute_ensemble_mean(int Nens,
+compute_ensemble_mean(const int da_iter,
+                      const int Nens,
                       const std::string& pf_name,
                       const amrex::Vector<std::string>& varnames);
 
 void Apply_H(const amrex::MultiFab& x_mf, amrex::MultiFab& y_mf);
 
 void read_in_observations(const int& da_iter, amrex::MultiFab& y_obs);
 
-std::vector<std::string> get_plotfile_list();
+std::vector<std::string> get_plotfile_list(const int da_iter);
 void read_in_observations(const int& da_iter,
                           const amrex::Vector<std::string>& varnames,
                           amrex::MultiFab& y_obs);
@@ -156,12 +158,14 @@ void compute_d_prime_vec(amrex::MultiFab& d_prime_vec,
                     const amrex::MultiFab& d_vec,
                     const amrex::Vector<amrex::Real>& R_diag);
 
-void compute_mean_H_xf(amrex::MultiFab& mean_H_xf,
+void compute_mean_H_xf(const int da_iter,
+                       amrex::MultiFab& mean_H_xf,
                        const int Nens,
                        const std::string& last_pf_name,
                        const amrex::Vector<std::string>& varnames);
 
-void compute_r_vec (int Nens,
+void compute_r_vec (const int da_iter,
+                    int Nens,
                     const std::string& last_pf_name,
                     const amrex::Vector<std::string>& varnames,
                     const amrex::MultiFab& mean_H_xf,
@@ -173,14 +177,16 @@ void compute_alpha_vec (const int& Nens,
                         const amrex::Vector<amrex::Real>& r_vec,
                         amrex::Vector<amrex::Real>& alpha_vec);
 
-void compute_Xf_prime_times_vector (const int Nens,
-                            const std::string& last_pf_name,
-                            const amrex::Vector<std::string>& varnames,
-                            const amrex::MultiFab& xf_bar,
-                            const amrex::Vector<amrex::Real>& vec_in,
-                            amrex::MultiFab& result);
+void compute_Xf_prime_times_vector (const int da_iter,
+                                    const int Nens,
+                                    const std::string& last_pf_name,
+                                    const amrex::Vector<std::string>& varnames,
+                                    const amrex::MultiFab& xf_bar,
+                                    const amrex::Vector<amrex::Real>& vec_in,
+                                    amrex::MultiFab& result);
 
-void update_ensemble (const int Nens,
+void update_ensemble (const int da_iter,
+                      const int Nens,
                       const std::string& last_pf_name,
                       const amrex::Vector<std::string>& varnames,
                       const amrex::MultiFab& xf_bar,
@@ -200,12 +206,16 @@ Matrix matrix_multiply(const Matrix& A, const Matrix& B);
 void matrix_print(const Matrix& A);
 void lapack_testing();
 
-void compute_S_matrix(Matrix& S,
+void compute_S_matrix(const int da_iter,
+                      Matrix& S,
                       const int& Nens,
                       const amrex::MultiFab& mean_H_xf,
                       const amrex::Vector<amrex::Real>& R_diag,
                       const std::string& last_pf_name,
                       const amrex::Vector<std::string>& varnames);
 std::string
 MakeEnsembleCheckpointName (int da_iter, int ens_no);
+
+std::string
+MakeFullPath (const std::string& filename, const int da_iter);
 #endif

.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp

@@ -140,10 +140,10 @@ WriteUpdatedEnsembleToERFClassData (const MultiFab& mf_cc_fine,
 }
 
 void
-ERF::ComputeAndWriteEnsemblePerturbations()
+ERF::ComputeAndWriteEnsemblePerturbations(const int da_iter)
 {
 
-    auto pltfiles = get_plotfile_list();
+    auto pltfiles = get_plotfile_list(da_iter);
 
     // Step 2: loop over all plotfiles (timestamps at which the plotfiles are written)
     // ie.find the ensemble mean at iteration 100, loop over the plt00100 file in each of the
@@ -153,13 +153,13 @@ ERF::ComputeAndWriteEnsemblePerturbations()
     const std::string member_prefix = "member_";
     for (const auto& pf_name : pltfiles)
     {
-        MultiFab mf_mean = compute_ensemble_mean(Nens, pf_name, varnames);
+        MultiFab mf_mean = compute_ensemble_mean(da_iter, Nens, pf_name, varnames);
         // -------------------------------
         // Step 3 & 4: Compute perturbations and write plotfiles
         // -------------------------------
         for (int n = 0; n < Nens; ++n)
         {
-            MultiFab mf_pert = read_member_multifab(n, pf_name, varnames);
+            MultiFab mf_pert = read_member_multifab(da_iter, n, pf_name, varnames);
             MultiFab::Subtract(mf_pert, mf_mean, 0, 0, mf_mean.nComp(), mf_mean.nGrow());
 
             // Create output directory
@@ -187,7 +187,7 @@ ERF::PerformDataAssimilation(int da_iter)
 {
     //lapack_testing();
 
-    auto pltfiles = get_plotfile_list();
+    auto pltfiles = get_plotfile_list(da_iter);
     std::string last_pf_name;
     if (!pltfiles.empty()) {
         last_pf_name = pltfiles.back();
@@ -205,12 +205,12 @@ ERF::PerformDataAssimilation(int da_iter)
     Vector<std::string> varnames = {"density","theta", "x_velocity","y_velocity","z_velocity", "qv", "qc", "qrain"};
 
     // Compute the ensemble mean
-    MultiFab xf_bar = compute_ensemble_mean(Nens, last_pf_name, varnames);
+    MultiFab xf_bar = compute_ensemble_mean(da_iter, Nens, last_pf_name, varnames);
 
     Print() << "Computing ensemble mean complete" << std::endl;
 
     // Construct perturbation plotfile name
-    std::string pltname = "plt_ens_mean";
+    std::string pltname = MakeFullPath("plt_ens_mean", da_iter);
     WriteSingleLevelPlotfile(pltname,
                              xf_bar,
                              varnames,
@@ -220,13 +220,13 @@ ERF::PerformDataAssimilation(int da_iter)
 
     // Compute the mean of forecast observations yf_bar = Hx_f
     MultiFab mean_H_xf;
-    compute_mean_H_xf(mean_H_xf, Nens, last_pf_name, varnames);
+    compute_mean_H_xf(da_iter, mean_H_xf, Nens, last_pf_name, varnames);
 
     Print() << "compute_mean_H_xf complete" << std::endl;
 
     Vector<std::string> varnames1 = {"density","theta", "x_velocity","y_velocity","z_velocity", "qv", "qc", "qrain"};
     // Construct perturbation plotfile name
-    std::string pltname1 = "plt_mean_H_xf";
+    std::string pltname1 = MakeFullPath("plt_mean_H_xf", da_iter);
     WriteSingleLevelPlotfile(pltname1,
                              mean_H_xf,
                              varnames1,
@@ -240,7 +240,7 @@ ERF::PerformDataAssimilation(int da_iter)
 
     Print() << "Observation reading complete" << std::endl;
 
-    /*std::string pltname2 = "plt_y_obs";
+    /*std::string pltname2 = WriteSingleLevelPlotfile("plt_y_obs", da_iter);
     Vector<std::string> varnames2 = {"x_velocity", "y_velocity"};
     WriteSingleLevelPlotfile(pltname2,
                             y_obs,
@@ -279,13 +279,13 @@ ERF::PerformDataAssimilation(int da_iter)
 
     // Compute r = Y'^Td'
     Vector<Real> r_vec;
-    compute_r_vec(Nens, last_pf_name, varnames, mean_H_xf, d_prime_vec, r_vec);
+    compute_r_vec(da_iter, Nens, last_pf_name, varnames, mean_H_xf, d_prime_vec, r_vec);
 
     Print() << "compute_r_vec complete" << std::endl;
 
     // Compute the S matrix
     Matrix S_mat(Nens);
-    compute_S_matrix(S_mat, Nens, mean_H_xf, R_diag, last_pf_name, varnames);
+    compute_S_matrix(da_iter, S_mat, Nens, mean_H_xf, R_diag, last_pf_name, varnames);
 
     Print() << "compute_S_matrix complete" << std::endl;
 
@@ -295,11 +295,11 @@ ERF::PerformDataAssimilation(int da_iter)
     Print() << "compute_alpha_vec complete" << std::endl;
 
     MultiFab Xf_prime_alpha;
-    compute_Xf_prime_times_vector(Nens, last_pf_name, varnames, xf_bar, alpha_vec, Xf_prime_alpha);
+    compute_Xf_prime_times_vector(da_iter, Nens, last_pf_name, varnames, xf_bar, alpha_vec, Xf_prime_alpha);
 
     Print() << "compute_Xf_prime_times_vector complete" << std::endl;
 
-    /*std::string pltname3 = "plt_Xf_prime_alpha";
+    /*std::string pltname3 = WriteSingleLevelPlotfile("plt_Xf_prime_alpha", da_iter);
     Vector<std::string> varnames3 = {"density", "theta", "x_velocity", "y_velocity", "z_velocity"};
     WriteSingleLevelPlotfile(pltname3,
                             Xf_prime_alpha,
@@ -314,11 +314,10 @@ ERF::PerformDataAssimilation(int da_iter)
 
     Print() << "add_multifabs(xf_bar, Xf_prime_alpha, xf_bar_updated) complete" << std::endl;
 
-    /*std::string pltname4 = "plt_xf_bar_updated";
-    Vector<std::string> varnames4 = {"density", "theta", "x_velocity", "y_velocity", "z_velocity"};
+    std::string pltname4 = MakeFullPath("plt_xf_bar_updated", da_iter);
     WriteSingleLevelPlotfile(pltname4,
                             xf_bar_updated,
-                            varnames4,
+                            varnames,
                             geom[0],
                             0.0,   // time
                             0);    // level*/
@@ -332,7 +331,7 @@ ERF::PerformDataAssimilation(int da_iter)
     for(int n=0; n< Nens; n++) {
         Print() << "Updating for ensemble " << n << std::endl;
         MultiFab mf_ens_pert;
-        update_ensemble(Nens, last_pf_name, varnames, xf_bar, T_mat, n, mf_ens_pert);
+        update_ensemble(da_iter, Nens, last_pf_name, varnames, xf_bar, T_mat, n, mf_ens_pert);
 
         // Update the ensemble
         MultiFab mf_ens_updated;
@@ -348,6 +347,16 @@ ERF::PerformDataAssimilation(int da_iter)
 
         m_plot3d_int_1 = -1;
         m_check_int = -1;
+        // InitData() handles all the initialization including when doing a
+        // restart. But when doing data assimilation, in the inputs file, we do not specify a
+        // erf.restart. So, the ERF class does not know it is a restart.
+        // So, whenever InitData() is called without assigning a restart_chkfile,
+        // it just does the from scratch initialization. This function is called by the
+        // temporararily created erf class - tmp_erf, from main.cpp. We simply want to fill all the data
+        // structures here by calling InitData(), and then the WriteUpdatedEnsembleToERFClassData will write
+        // all the updated ensemble data into the lev_new data structure, which will then be
+        // written into the checkpoint files by MakeEnsembleCheckpointName.
+
         InitData();
         auto& lev_new = vars_new[0];
 
@@ -363,31 +372,3 @@ ERF::PerformDataAssimilation(int da_iter)
         WriteCheckpointFile();
     }
 }
-
-
-void
-ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
-{
-    // --------------------------------------------------------
-    // Set up member-specific output directories
-    // --------------------------------------------------------
-    std::string member_dir;
-
-    std::stringstream ss;
-    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
-    member_dir = ss.str();
-
-    if (ParallelDescriptor::IOProcessor())
-    {
-        fs::create_directories(member_dir + "/plotfiles");
-        fs::create_directories(member_dir + "/chkfiles");
-        fs::create_directories(member_dir + "/pertfiles");
-    }
-
-    ParallelDescriptor::Barrier();
-
-    // Only change the basename/prefix.
-    // ERF will handle the step numbering as before.
-    check_file = member_dir + "/chkfiles/chk";
-    plot3d_file_1  = member_dir + "/plotfiles/plt";
-}

.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_SquallLine2D.cpp

@@ -96,6 +96,12 @@ void compute_R_diag_vals(Vector<Real>& R_diag)
     // Set the size to 2
     R_diag.resize(8,0.01);
 
-    R_diag[0] = 0.01;
-    R_diag[1] = 0.01;
+    R_diag[0] = 4e-4; //density
+    R_diag[1] = 1.0; // theta
+    R_diag[2] = 1.0; // x_velocity
+    R_diag[3] = 1.0; // y_velocity
+    R_diag[4] = 0.25; // z_velocity
+    R_diag[5] = 2.5e-7; // qv
+    R_diag[6] = 1e-8; // qc
+    R_diag[7] = 1e-8; // qrain
 }