Bug fix: set_bc_dir() needed the current solution arrays (solutions.current.A/Y/D) to be allocated and initialized before it was called

Author: shiyi

Code/Source/solver/Integrator.cpp

@@ -50,17 +50,12 @@ void Integrator::initialize_arrays() {
   Ag_.resize(tDof, tnNo);
   Yg_.resize(tDof, tnNo);
   Dg_.resize(tDof, tnNo);
-  solutions_.current.A.resize(tDof, tnNo);
-  solutions_.current.D.resize(tDof, tnNo);
-  solutions_.current.Y.resize(tDof, tnNo);
   res_.resize(nFacesLS);
   incL_.resize(nFacesLS);
 
-  // old solution already initialized via move in constructor
-  // Initialize current solution from old solution
-  solutions_.current.A = solutions_.old.A;
-  solutions_.current.D = solutions_.old.D;
-  solutions_.current.Y = solutions_.old.Y;
+  // Solution arrays (old and current) are already initialized in initialize()
+  // and boundary conditions have been applied via set_bc_dir()
+  // Do NOT reinitialize them here to preserve those boundary values
 }
 
 //------------------------

Code/Source/solver/initialize.cpp

@@ -643,6 +643,12 @@ void initialize(Simulation* simulation, Vector<double>& timeP)
   solutions.current.Y.resize(tDof, tnNo);
   solutions.current.D.resize(tDof, tnNo);
 
+  // Initialize current solution from old solution
+  // (will be modified by set_bc_dir later, before Integrator creation)
+  solutions.current.A = solutions.old.A;
+  solutions.current.Y = solutions.old.Y;
+  solutions.current.D = solutions.old.D;
+
   com_mod.Bf.resize(nsd,tnNo);
 
   // [TODO] DaveP not implemented?
@@ -770,8 +776,9 @@ void initialize(Simulation* simulation, Vector<double>& timeP)
     rmsh.D0 = solutions.old.D;
   } // resetSim
 
-  // Initialize new variables
-  // An, Yn, Dn moved to Integrator class (initialized there from Ao, Yo, Do)
+  // An, Yn, Dn are now part of SolutionStates struct (solutions.current)
+  // They are initialized from Ao, Yo, Do earlier in this function
+  // and will be modified by set_bc_dir() before being passed to Integrator
 
   for (int iM = 0; iM < nMsh; iM++) { 
     if (cm.mas(cm_mod)) {
@@ -840,13 +847,13 @@ void initialize(Simulation* simulation, Vector<double>& timeP)
     post::ppbin2vtk(simulation);
   }
 
-  // Making sure the old solution satisfies BCs
+  // Making sure the solution satisfies BCs
   //
-  // Modifies Ao, Yo, Do (via solutions struct)
+  // Modifies both old (Ao, Yo, Do) and current (An, Yn, Dn) solution arrays
   //
   set_bc::set_bc_dir(com_mod, solutions);
 
-  // Preparing TXT files (pass solutions since An, Dn, and Yn haven't been created in Integrator yet)
+  // Preparing TXT files
   txt_ns::txt(simulation, true, solutions);
 
   // Printing the first line and initializing timeP