DataAnalyticsEngineering · sanathkeshav · Oct 30, 2025 · Oct 29, 2025 · Oct 29, 2025 · Oct 29, 2025
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # FANS Changelog
 
+## latest
+
+- Bugfix: Fix bad initial guess for time steps > 0 [#109](https://github.com/DataAnalyticsEngineering/FANS/pull/109)
+
 ## v0.5.1
 
 - Optimize HDF5 I/O performance by using collective I/O operations [#105](https://github.com/DataAnalyticsEngineering/FANS/pull/105)

diff --git a/include/material_models/CompressibleNeoHookean.h b/include/material_models/CompressibleNeoHookean.h
@@ -31,16 +31,18 @@ class CompressibleNeoHookean : public LargeStrainMechModel {
             lambda[i] = bulk_modulus[i] - (2.0 / 3.0) * mu[i];
         }
 
-        // Compute reference tangent at F=I
-        kapparef_mat = Matrix<double, 9, 9>::Zero();
-        for (int mat_idx = 0; mat_idx < n_mat; ++mat_idx) {
-            Matrix3d             F_identity = Matrix3d::Identity();
-            Matrix3d             S          = compute_S(F_identity, mat_idx, 0);
-            Matrix<double, 6, 6> C_mandel   = compute_material_tangent(F_identity, mat_idx);
-            Matrix<double, 9, 9> A          = compute_spatial_tangent(F_identity, S, C_mandel);
-            kapparef_mat += A;
+        if (kapparef_mat.isZero()) {
+            // Compute reference tangent at F=I
+            kapparef_mat = Matrix<double, 9, 9>::Zero();
+            for (int mat_idx = 0; mat_idx < n_mat; ++mat_idx) {
+                Matrix3d             F_identity = Matrix3d::Identity();
+                Matrix3d             S          = compute_S(F_identity, mat_idx, 0);
+                Matrix<double, 6, 6> C_mandel   = compute_material_tangent(F_identity, mat_idx);
+                Matrix<double, 9, 9> A          = compute_spatial_tangent(F_identity, S, C_mandel);
+                kapparef_mat += A;
+            }
+            kapparef_mat /= static_cast<double>(n_mat);
         }
-        kapparef_mat /= static_cast<double>(n_mat);
     }
 
     Matrix3d compute_S(const Matrix3d &F, int mat_index, ptrdiff_t element_idx) override

diff --git a/include/material_models/GBDiffusion.h b/include/material_models/GBDiffusion.h
@@ -95,7 +95,7 @@ class GBDiffusion : public ThermalModel, public LinearModel<1, 3> {
             throw std::runtime_error("Error in GBDiffusion initialization: " + std::string(e.what()));
         }
 
-        kapparef_mat = Matrix3d::Zero(3, 3);
+        Matrix3d kappa_temp = Matrix3d::Zero();
         Matrix3d phase_kappa;
         phase_stiffness = new Matrix<double, 8, 8>[n_mat];
 
@@ -112,12 +112,14 @@ class GBDiffusion : public ThermalModel, public LinearModel<1, 3> {
             } else {
                 throw std::runtime_error("GBDiffusion: Unknown material index");
             }
-            kapparef_mat += phase_kappa;
+            kappa_temp += phase_kappa;
             for (int p = 0; p < n_gp; ++p) {
                 phase_stiffness[i] += B_int[p].transpose() * phase_kappa * B_int[p] * v_e / n_gp;
             }
         }
-        kapparef_mat /= n_mat;
+        if (kapparef_mat.isZero()) {
+            kapparef_mat = kappa_temp / n_mat;
+        }
     }
     ~GBDiffusion()
     {

diff --git a/include/material_models/J2Plasticity.h b/include/material_models/J2Plasticity.h
@@ -22,12 +22,14 @@ class J2Plasticity : public SmallStrainMechModel {
         }
         n_mat = bulk_modulus.size();
 
-        const double Kbar      = std::accumulate(bulk_modulus.begin(), bulk_modulus.end(), 0.0) / static_cast<double>(n_mat);
-        const double Gbar      = std::accumulate(shear_modulus.begin(), shear_modulus.end(), 0.0) / static_cast<double>(n_mat);
-        const double lambdabar = Kbar - 2.0 * Gbar / 3.0;
-        kapparef_mat.setZero();
-        kapparef_mat.topLeftCorner(3, 3).setConstant(lambdabar);
-        kapparef_mat.diagonal().array() += 2.0 * Gbar;
+        if (kapparef_mat.isZero()) {
+            const double Kbar      = std::accumulate(bulk_modulus.begin(), bulk_modulus.end(), 0.0) / static_cast<double>(n_mat);
+            const double Gbar      = std::accumulate(shear_modulus.begin(), shear_modulus.end(), 0.0) / static_cast<double>(n_mat);
+            const double lambdabar = Kbar - 2.0 * Gbar / 3.0;
+            kapparef_mat.setZero();
+            kapparef_mat.topLeftCorner(3, 3).setConstant(lambdabar);
+            kapparef_mat.diagonal().array() += 2.0 * Gbar;
+        }
 
         // Allocate the member matrices/vectors for performance optimization
         sqrt_two_over_three = sqrt(2.0 / 3.0);

diff --git a/include/material_models/LinearElastic.h b/include/material_models/LinearElastic.h
@@ -24,16 +24,18 @@ class LinearElasticIsotropic : public SmallStrainMechModel, public LinearModel<3
             lambda[i] = bulk_modulus[i] - (2.0 / 3.0) * mu[i];
         }
 
-        double lambda_ref = (*max_element(lambda.begin(), lambda.end()) +
-                             *min_element(lambda.begin(), lambda.end())) /
+        if (kapparef_mat.isZero()) {
+            double lambda_ref = (*max_element(lambda.begin(), lambda.end()) +
+                                 *min_element(lambda.begin(), lambda.end())) /
+                                2;
+            double mu_ref = (*max_element(mu.begin(), mu.end()) +
+                             *min_element(mu.begin(), mu.end())) /
                             2;
-        double mu_ref = (*max_element(mu.begin(), mu.end()) +
-                         *min_element(mu.begin(), mu.end())) /
-                        2;
 
-        kapparef_mat = Matrix<double, 6, 6>::Zero();
-        kapparef_mat.topLeftCorner(3, 3).setConstant(lambda_ref);
-        kapparef_mat += 2 * mu_ref * Matrix<double, 6, 6>::Identity();
+            kapparef_mat = Matrix<double, 6, 6>::Zero();
+            kapparef_mat.topLeftCorner(3, 3).setConstant(lambda_ref);
+            kapparef_mat += 2 * mu_ref * Matrix<double, 6, 6>::Identity();
+        }
 
         phase_stiffness = new Matrix<double, 24, 24>[n_mat];
         Matrix<double, 6, 6> phase_kappa;
@@ -105,7 +107,7 @@ class LinearElasticTriclinic : public SmallStrainMechModel, public LinearModel<3
 
         // Assemble stiffness matrices for each material
         C_mats.resize(n_mat);
-        kapparef_mat = Matrix<double, 6, 6>::Zero();
+        Matrix<double, 6, 6> kappa_temp = Matrix<double, 6, 6>::Zero();
 
         for (size_t i = 0; i < n_mat; ++i) {
             Matrix<double, 6, 6> C_i = Matrix<double, 6, 6>::Zero();
@@ -122,10 +124,12 @@ class LinearElasticTriclinic : public SmallStrainMechModel, public LinearModel<3
             C_i = C_i.selfadjointView<Eigen::Upper>();
 
             C_mats[i] = C_i;
-            kapparef_mat += C_i;
+            kappa_temp += C_i;
         }
 
-        kapparef_mat /= n_mat;
+        if (kapparef_mat.isZero()) {
+            kapparef_mat = kappa_temp / n_mat;
+        }
 
         // Compute phase stiffness matrices
         phase_stiffness = new Matrix<double, 24, 24>[n_mat];

diff --git a/include/material_models/LinearThermal.h b/include/material_models/LinearThermal.h
@@ -16,10 +16,12 @@ class LinearThermalIsotropic : public ThermalModel, public LinearModel<1, 3> {
         }
         n_mat = conductivity.size();
 
-        double kappa_ref = (*max_element(conductivity.begin(), conductivity.end()) +
-                            *min_element(conductivity.begin(), conductivity.end())) /
-                           2;
-        kapparef_mat = kappa_ref * Matrix3d::Identity();
+        if (kapparef_mat.isZero()) {
+            double kappa_ref = (*max_element(conductivity.begin(), conductivity.end()) +
+                                *min_element(conductivity.begin(), conductivity.end())) /
+                               2;
+            kapparef_mat = kappa_ref * Matrix3d::Identity();
+        }
 
         Matrix3d phase_kappa;
         phase_stiffness = new Matrix<double, 8, 8>[n_mat];

diff --git a/include/material_models/PseudoPlastic.h b/include/material_models/PseudoPlastic.h
@@ -32,12 +32,14 @@ class PseudoPlastic : public SmallStrainMechModel {
         }
         n_mat = bulk_modulus.size();
 
-        const double Kbar      = std::accumulate(bulk_modulus.begin(), bulk_modulus.end(), 0.0) / static_cast<double>(n_mat);
-        const double Gbar      = std::accumulate(shear_modulus.begin(), shear_modulus.end(), 0.0) / static_cast<double>(n_mat);
-        const double lambdabar = Kbar - 2.0 * Gbar / 3.0;
-        kapparef_mat.setZero();
-        kapparef_mat.topLeftCorner(3, 3).setConstant(lambdabar);
-        kapparef_mat.diagonal().array() += 2.0 * Gbar;
+        if (kapparef_mat.isZero()) {
+            const double Kbar      = std::accumulate(bulk_modulus.begin(), bulk_modulus.end(), 0.0) / static_cast<double>(n_mat);
+            const double Gbar      = std::accumulate(shear_modulus.begin(), shear_modulus.end(), 0.0) / static_cast<double>(n_mat);
+            const double lambdabar = Kbar - 2.0 * Gbar / 3.0;
+            kapparef_mat.setZero();
+            kapparef_mat.topLeftCorner(3, 3).setConstant(lambdabar);
+            kapparef_mat.diagonal().array() += 2.0 * Gbar;
+        }
     }
 
     void initializeInternalVariables(ptrdiff_t num_elements, int num_gauss_points) override

diff --git a/include/material_models/SaintVenantKirchhoff.h b/include/material_models/SaintVenantKirchhoff.h
@@ -29,16 +29,18 @@ class SaintVenantKirchhoff : public LargeStrainMechModel {
             lambda[i] = bulk_modulus[i] - (2.0 / 3.0) * mu[i];
         }
 
-        // Compute reference tangent at F=I
-        kapparef_mat = Matrix<double, 9, 9>::Zero();
-        for (int mat_idx = 0; mat_idx < n_mat; ++mat_idx) {
-            Matrix3d             F_identity = Matrix3d::Identity();
-            Matrix3d             S          = compute_S(F_identity, mat_idx, 0);
-            Matrix<double, 6, 6> C_mandel   = compute_material_tangent(F_identity, mat_idx);
-            Matrix<double, 9, 9> A          = compute_spatial_tangent(F_identity, S, C_mandel);
-            kapparef_mat += A;
+        if (kapparef_mat.isZero()) {
+            // Compute reference tangent at F=I
+            kapparef_mat = Matrix<double, 9, 9>::Zero();
+            for (int mat_idx = 0; mat_idx < n_mat; ++mat_idx) {
+                Matrix3d             F_identity = Matrix3d::Identity();
+                Matrix3d             S          = compute_S(F_identity, mat_idx, 0);
+                Matrix<double, 6, 6> C_mandel   = compute_material_tangent(F_identity, mat_idx);
+                Matrix<double, 9, 9> A          = compute_spatial_tangent(F_identity, S, C_mandel);
+                kapparef_mat += A;
+            }
+            kapparef_mat /= static_cast<double>(n_mat);
         }
-        kapparef_mat /= static_cast<double>(n_mat);
     }
 
     Matrix3d compute_S(const Matrix3d &F, int mat_index, ptrdiff_t element_idx) override

diff --git a/include/matmodel.h b/include/matmodel.h
@@ -85,6 +85,29 @@ Matmodel<howmany, n_str>::Matmodel(const Reader &reader)
     eps.resize(n_str * n_gp);
     g0.resize(n_str * n_gp);
     sigma.resize(n_str * n_gp);
+
+    kapparef_mat.setZero();
+
+    // Check for optional user-defined reference material
+    if (reader.materialProperties.contains("reference_material")) {
+        auto ref_mat = reader.materialProperties["reference_material"].get<vector<vector<double>>>();
+        if (ref_mat.size() != n_str || ref_mat[0].size() != n_str) {
+            throw std::invalid_argument("reference_material must be " + std::to_string(n_str) + "x" + std::to_string(n_str));
+        }
+        for (int i = 0; i < n_str; ++i) {
+            for (int j = 0; j < n_str; ++j) {
+                kapparef_mat(i, j) = ref_mat[i][j];
+            }
+        }
+        // Verify SPD using Cholesky decomposition
+        Eigen::LLT<Matrix<double, n_str, n_str>> llt(kapparef_mat);
+        if (llt.info() != Eigen::Success) {
+            throw std::invalid_argument("reference_material must be symmetric positive definite");
+        }
+        if (reader.world_rank == 0) {
+            cout << "# Using user-defined reference material for fundamental solution." << endl;
+        }
+    }
 }
 
 template <int howmany, int n_str>

diff --git a/include/mixedBCs.h b/include/mixedBCs.h
@@ -153,7 +153,8 @@ struct MixedBCController {
   protected:
     MixedBC  mbc_local;
     size_t   step_idx = 0;
-    VectorXd g0_vec; // current macro strain (size n_str)
+    VectorXd g0_vec;      // current macro strain (size n_str)
+    VectorXd g0_vec_prev; // previous time step for extrapolation
 
     // call from user code after v_u update each iteration
     template <typename SolverType>
@@ -198,10 +199,18 @@ struct MixedBCController {
                 // Small strain: zero initialization
                 g0_vec = VectorXd::Zero(n_str);
             }
+            g0_vec_prev = g0_vec; // Initialize history
+        } else {
+            // Linear extrapolation for stress-controlled components (t > 0)
+            VectorXd delta = g0_vec - g0_vec_prev;
+            g0_vec_prev    = g0_vec; // Store current before updating
+
+            for (int i = 0; i < mbc_local.idx_F.size(); ++i) {
+                g0_vec(mbc_local.idx_F(i)) += delta(mbc_local.idx_F(i));
+            }
         }
 
-        // Update ONLY the strain-controlled components from the current time step
-        // Stress-controlled components remain from previous converged state
+        // Update strain-controlled components from prescribed path
         if (mbc_local.idx_E.size()) {
             VectorXd prescribed = mbc_local.F_E_path.row(t).transpose();
             for (int i = 0; i < mbc_local.idx_E.size(); ++i) {

diff --git a/include/solver.h b/include/solver.h
@@ -31,6 +31,7 @@ class Solver : private MixedBCController<howmany> {
     unsigned short *ms;  // Micro-structure
     double         *v_r; //!< Residual vector
     double         *v_u;
+    double         *v_u_prev; //!< Previous displacement for extrapolation
     double         *buffer_padding;
 
     RealArray      v_r_real; // can't do the "classname()" intialization here, and Map doesn't have a default constructor
@@ -44,7 +45,8 @@ class Solver : private MixedBCController<howmany> {
     void iterateCubes(F f);
 
     void         solve();
-    virtual void internalSolve() {}; // important to have "{}" here, otherwise we get an error about undefined reference to vtable
+    void         extrapolateDisplacement(); //!< Linear extrapolation for next time step
+    virtual void internalSolve() {};        // important to have "{}" here, otherwise we get an error about undefined reference to vtable
 
     template <int padding, typename F>
     void compute_residual_basic(RealArray &r_matrix, RealArray &u_matrix, F f);
@@ -120,6 +122,7 @@ Solver<howmany, n_str>::Solver(Reader &reader, Matmodel<howmany, n_str> *mat)
 
       v_u(fftw_alloc_real((local_n0 + 1) * n_y * n_z * howmany)),
       v_u_real(v_u, n_z * howmany, local_n0 * n_y, OuterStride<>(n_z * howmany)),
+      v_u_prev(fftw_alloc_real(local_n0 * n_y * n_z * howmany)),
 
       rhat((std::complex<double> *) v_r, local_n1 * n_x * (n_z / 2 + 1) * howmany), // actual initialization is below
       buffer_padding(fftw_alloc_real(n_y * (n_z + 2) * howmany))
@@ -128,6 +131,7 @@ Solver<howmany, n_str>::Solver(Reader &reader, Matmodel<howmany, n_str> *mat)
     for (ptrdiff_t i = local_n0 * n_y * n_z * howmany; i < (local_n0 + 1) * n_y * n_z * howmany; i++) {
         this->v_u[i] = 0;
     }
+    std::memset(v_u_prev, 0, local_n0 * n_y * n_z * howmany * sizeof(double));
 
     matmodel->initializeInternalVariables(local_n0 * n_y * n_z, matmodel->n_gp);
 
@@ -291,6 +295,17 @@ void Solver<howmany, n_str>::solve()
     matmodel->updateInternalVariables();
 }
 
+template <int howmany, int n_str>
+void Solver<howmany, n_str>::extrapolateDisplacement()
+{
+    const size_t n = local_n0 * n_y * n_z * howmany;
+    for (size_t i = 0; i < n; ++i) {
+        const double delta = v_u[i] - v_u_prev[i];
+        v_u_prev[i]        = v_u[i];
+        v_u[i] += delta; // v_u = v_u + (v_u - v_u_prev)
+    }
+}
+
 template <int howmany, int n_str>
 template <int padding, typename F>
 void Solver<howmany, n_str>::iterateCubes(F f)
@@ -598,6 +613,7 @@ 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>
@@ -678,6 +694,10 @@ Solver<howmany, n_str>::~Solver()
         fftw_free(v_u);
         v_u = nullptr;
     }
+    if (v_u_prev) {
+        fftw_free(v_u_prev);
+        v_u_prev = nullptr;
+    }
     if (buffer_padding) {
         fftw_free(buffer_padding);
         buffer_padding = nullptr;

diff --git a/src/main.cpp b/src/main.cpp
@@ -15,7 +15,17 @@ void runSolver(Reader &reader)
         Matmodel<howmany, n_str> *matmodel = createMatmodel<howmany, n_str>(reader);
         Solver<howmany, n_str>   *solver   = createSolver(reader, matmodel);
 
+        if (reader.world_rank == 0) {
+            printf("\n╔════════════════════════════════════════════════════════════ Load case %zu/%zu: %zu time steps ════════════════════════════════════════════════════════════╗\n",
+                   load_path_idx + 1, reader.load_cases.size(), reader.load_cases[load_path_idx].n_steps);
+        }
+
         for (size_t time_step_idx = 0; time_step_idx < reader.load_cases[load_path_idx].n_steps; ++time_step_idx) {
+            if (reader.world_rank == 0) {
+                printf("║   ▶ Time step %zu/%zu (load case %zu/%zu) ◀ \n",
+                       time_step_idx + 1, reader.load_cases[load_path_idx].n_steps,
+                       load_path_idx + 1, reader.load_cases.size());
+            }
             if (reader.load_cases[load_path_idx].mixed) {
                 solver->enableMixedBC(reader.load_cases[load_path_idx].mbc, time_step_idx);
             } else {
@@ -27,6 +37,9 @@ void runSolver(Reader &reader)
         }
         delete solver;
         delete matmodel;
+        if (reader.world_rank == 0) {
+            printf("╚══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════╝\n");
+        }
     }
 }