Cleanup output, function for computing matrix norms.

LLNL · Oct 25, 2024 · 457cdb4 · 457cdb4
1 parent 54c9ee4
commit 457cdb4
Show file tree

Hide file tree

Showing 2 changed files with 36 additions and 17 deletions.
diff --git a/src/serac/numerics/equation_solver.cpp b/src/serac/numerics/equation_solver.cpp
@@ -119,8 +119,10 @@ class NewtonSolver : public mfem::NewtonSolver {
       }
 
       if (norm != norm) {
-        mfem::out << "Initial residual for Newton iteration is undefined/nan." << std::endl;
-        mfem::out << "Newton: No convergence!\n";
+        if (print_options.first_and_last) {
+          mfem::out << "Initial residual for Newton iteration is undefined/nan." << std::endl;
+          mfem::out << "Newton: No convergence!\n";
+        }
         throw SolveException("Initial residual for Newton iteration is undefined/nan.");
         return;
       }
@@ -542,12 +544,14 @@ class TrustRegion : public mfem::NewtonSolver {
         } else {
           mfem::out << ", norm goal = " << std::setw(13) << norm_goal;
         }
-        mfem::out << '\n';
+        mfem::out << " :: ";
       }
 
       if (norm != norm) {
-        mfem::out << "Initial residual for trust-region iteration is undefined/nan." << std::endl;
-        mfem::out << "Newton: No convergence!\n";
+        if (print_options.first_and_last) {
+          mfem::out << "Initial residual for trust-region iteration is undefined/nan." << std::endl;
+          mfem::out << "Newton: No convergence!\n";
+        }
         throw SolveException("Initial residual for trust-region iteration is undefined/nan.");
         return;
       }

diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp
@@ -27,6 +27,15 @@
 
 namespace serac {
 
+inline double matrixNorm(std::unique_ptr<mfem::HypreParMatrix>& K, const FiniteElementState& x, int iters)
+{
+  mfem::HypreParMatrix* H = K.get();
+  hypre_ParCSRMatrix * Hhypre = static_cast<hypre_ParCSRMatrix *>(*H);
+  double Hfronorm;
+  hypre_ParCSRMatrixNormFro(Hhypre, &Hfronorm);
+  return Hfronorm;
+}
+
 namespace solid_mechanics {
 
 namespace detail {
@@ -872,9 +881,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
         dx_dX += du_dX;
       }
 
-      auto flux = dot(stress); // * det(dx_dX);
-
-      return serac::tuple{material_.density * d2u_dt2, flux};
+      return serac::tuple{material_.density * d2u_dt2, stress};
     }
   };
 
@@ -1163,6 +1170,14 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
           J_.reset();
           J_ = assemble(drdu);
 
+          //auto Jt =  std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());
+          //Jt  = std::unique_ptr<mfem::HypreParMatrix>(mfem::Add(0.5, *J_, -0.5, *Jt));
+          //double n1 = matrixNorm(J_, displacement_, 10);
+          //double n2 = matrixNorm(Jt, displacement_, 10);
+          //if (mpi_rank_==0) {
+          //  std::cout << "approx norms = " << n1 << ", " << n2 << " :: ";
+          //}
+
           if (symmetrize) {
             auto J_T = std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());
             J_  = std::unique_ptr<mfem::HypreParMatrix>(mfem::Add(0.5, *J_, 0.5, *J_T));
@@ -1622,7 +1637,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
   GeometricNonlinearities geom_nonlin_;
 
   /// @brief A flag denoting whether to symmetrize linear operators passed to the solvers
-  bool symmetrize = true;
+  bool symmetrize = false; //true;
 
   /// @brief A flag denoting whether to compute the warm start for improved robustness
   bool use_warm_start_;
@@ -1645,9 +1660,9 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
       }...};
 
   /// @brief Solve the Quasi-static Newton system
-  virtual void quasiStaticSolve(double dt, double a, double b)
+  virtual void quasiStaticSolve(double dt, double a, double b, int level=0)
   {
-    if (b < 1e-5) {
+    if (level >= 6) {
       std::cout << "Too many boundary condition cutbacks, try increasing the number of load steps " << std::endl;
       return;
     }
@@ -1660,18 +1675,18 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
       warmStartDisplacement(dt, b);
       nonlin_solver_->solve(displacement_);
     } catch (const std::exception& e) {
-      std::cout << "caught: " << e.what() << std::endl;
+      if (mpi_rank_==0) mfem::out << "caught: " << e.what() << std::endl;
       displacement_ -= du_;
       solver_success = false;
-      quasiStaticSolve(dt, 1.0, 0.5*b);
-      quasiStaticSolve(dt, 1.0, 1.0);
+      quasiStaticSolve(dt, 1.0, 0.5*b, level+1);
+      quasiStaticSolve(dt, 1.0, 1.0, level+1);
     }
 
     if (solver_success) {
       if (b==1.0) {
-        std::cout << "final solve succeeded for time " << time_ << " dt = " << dt << std::endl;
+        if (mpi_rank_==0) mfem::out << "final solve succeeded for time " << time_ << " dt = " << dt << std::endl;
       } else {
-        std::cout << "substep solve succeeded for time " << time_ << " dt = " << dt << std::endl;
+        if (mpi_rank_==0) mfem::out << "substep solve succeeded for time " << time_ << " dt = " << dt << std::endl;
       }
     } 
   }
@@ -1775,7 +1790,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
   {
     SERAC_MARK_FUNCTION;
 
-    std::cout << "Solving with displacement factor = " << displacement_scale_factor << std::endl;
+    if (mpi_rank_==0) mfem::out << "Solving with displacement factor = " << displacement_scale_factor << std::endl;
 
     du_ = 0.0;
     for (auto& bc : bcs_.essentials()) {