warm lsearch

Author: sanathkeshav

include/reader.h

@@ -30,6 +30,9 @@ class Reader {
     json             materialProperties;
     double           TOL;
     json             errorParameters;
+    int              ls_max_iter{5};
+    double           ls_tol{1e-2};
+    bool             extrapolate_displacement{true};
     json             microstructure;
     int              n_it;
     vector<LoadCase> load_cases;

include/solver.h

@@ -702,7 +702,6 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, int load_idx, int time_
         }
         reader.writeData("homogenized_tangent", load_idx, time_idx, homogenized_tangent.data(), dims, 2);
     }
-    extrapolateDisplacement(); // prepare v_u for next time step
 }
 
 template <int howmany, int n_str>

include/solverCG.h

@@ -23,6 +23,7 @@ class SolverCG : public Solver<howmany, n_str> {
     RealArray s_real;
     RealArray d_real;
     RealArray rnew_real;
+    double    alpha_warm{0.1};
 
     void   internalSolve();
     void   LineSearchSecant();
@@ -64,6 +65,7 @@ void SolverCG<howmany, n_str>::internalSolve()
         printf("\n# Start FANS - Conjugate Gradient Solver \n");
 
     bool islinear = this->matmanager->all_linear;
+    alpha_warm    = 0.1;
 
     s_real.setZero();
     d_real.setZero();
@@ -118,11 +120,11 @@ template <int howmany, int n_str>
 void SolverCG<howmany, n_str>::LineSearchSecant()
 {
     double       err        = 10.0;
-    const int    MaxIter    = 5;
-    const double tol        = 1e-2;
+    const int    MaxIter    = this->reader.ls_max_iter;
+    const double tol        = this->reader.ls_tol;
     int          _iter      = 0;
     double       alpha_prev = 0.0;
-    double       alpha_curr = 0.1;
+    double       alpha_curr = alpha_warm;
 
     double rpd = dotProduct(v_r_real, d_real);
     v_u_real += d_real * alpha_curr;
@@ -151,7 +153,8 @@ void SolverCG<howmany, n_str>::LineSearchSecant()
         this->template compute_residual<0>(rnew_real, v_u_real);
         r1pd = dotProduct(rnew_real, d_real);
     }
-    v_r_real = rnew_real;
+    alpha_warm = (_iter == MaxIter && err > tol) ? 0.1 : alpha_curr;
+    v_r_real   = rnew_real;
     if (this->world_rank == 0)
         printf("line search iter %i, alpha %f - error %e - ", _iter, alpha_curr, err);
 }

src/main.cpp

@@ -34,6 +34,8 @@ void runSolver(Reader &reader)
             }
             solver->solve();
             solver->postprocess(reader, load_path_idx, time_step_idx);
+            if (reader.extrapolate_displacement)
+                solver->extrapolateDisplacement();
         }
         delete solver;
         delete matmanager;

src/reader.cpp

@@ -75,16 +75,11 @@ void Reader ::ReadInputFile(char input_fn[])
         microstructure = j["microstructure"];
         std::snprintf(ms_filename, sizeof(ms_filename), "%s", microstructure["filepath"].get<std::string>().c_str());
         // dataset name handling
-        std::string tmp_str = microstructure["datasetname"].get<std::string>();
-        if (tmp_str.empty()) {
-            throw std::invalid_argument(
-                "datasetname must not be empty and must refer to a valid HDF5 path");
-        }
-        // Ensure absolute HDF5 path, leeading slash
-        if (tmp_str.front() != '/') {
-            tmp_str.insert(tmp_str.begin(), '/');
-        }
-        std::snprintf(ms_datasetname, sizeof(ms_datasetname), "%s", tmp_str.c_str());
+        const auto tmp_str = microstructure["datasetname"].get<std::string>();
+        if (tmp_str.empty())
+            throw std::invalid_argument("datasetname must not be empty and must refer to a valid HDF5 path");
+        // Ensure absolute HDF5 path, leading slash
+        std::snprintf(ms_datasetname, sizeof(ms_datasetname), "%s%s", tmp_str.front() == '/' ? "" : "/", tmp_str.c_str());
         L = microstructure["L"].get<vector<double>>();
 
         if (j.contains("results_prefix")) {
@@ -100,6 +95,15 @@ void Reader ::ReadInputFile(char input_fn[])
         TOL             = errorParameters["tolerance"].get<double>();
         n_it            = j["n_it"].get<int>();
 
+        extrapolate_displacement = j.value("extrapolate_displacement", extrapolate_displacement);
+
+        if (j.contains("linesearch_parameters")) {
+            ls_max_iter = j["linesearch_parameters"].value("max_iter", ls_max_iter);
+            ls_tol      = j["linesearch_parameters"].value("tol", ls_tol);
+            if (ls_max_iter < 1 || ls_tol <= 0.0)
+                throw std::invalid_argument("linesearch_parameters: max_iter >= 1 and tol > 0 required");
+        }
+
         problemType = j["problem_type"].get<string>();
         method      = j["method"].get<string>();
 
@@ -173,7 +177,7 @@ void Reader ::ReadInputFile(char input_fn[])
         }
 
     } catch (const std::exception &e) {
-        fprintf(stderr, "ERROR trying to read input file '%s' for FANS\n", input_fn);
+        fprintf(stderr, "ERROR trying to read input file '%s' for FANS: %s\n", input_fn, e.what());
         exit(10);
     }
 }