Merge pull request #1191 from chleh/tes-for-pr-side

Misc changes

Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h

@@ -397,6 +397,9 @@ loop(ProjectData& project)
         t        = ts.current();
         timestep = ts.steps();
 
+        INFO("=== timestep #%u (t=%gs, dt=%gs) ==============================",
+             timestep, t, delta_t);
+
         // TODO use process name
         unsigned pcs_idx = 0;
         for (auto p = project.processesBegin(); p != project.processesEnd();

MathLib/LinAlg/BLAS.cpp

@@ -323,6 +323,11 @@ void matMultAdd(EigenMatrix const& A, EigenVector const& v1, EigenVector const&
     v3.getRawVector() = v2.getRawVector() + A.getRawMatrix()*v1.getRawVector();
 }
 
+void finalizeAssembly(EigenMatrix& x)
+{
+    x.getRawMatrix().makeCompressed();
+}
+
 } // namespace BLAS
 
 } // namespace MathLib

MathLib/LinAlg/BLAS.h

@@ -226,6 +226,8 @@ void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y);
 void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
                 EigenVector const& v2, EigenVector& v3);
 
+void finalizeAssembly(EigenMatrix& A);
+
 } // namespace BLAS
 
 } // namespace MathLib

MathLib/LinAlg/Lis/LisOption.h

@@ -42,16 +42,15 @@ struct LisOption
 {
 	LisOption(BaseLib::ConfigTree const* const options)
 	{
-		if (options)
-		{
+		if (options) {
 			ignoreOtherLinearSolvers(*options, "lis");
 			if (auto s = options->getConfParamOptional<std::string>("lis")) {
-				if (!s->empty()) _option_string += " " + *s;
+				if (!s->empty()) {
+					_option_string += " " + *s;
+					INFO("Lis options: \"%s\"", _option_string.c_str());
+				}
 			}
 		}
-#ifndef NDEBUG
-		INFO("Lis options: \"%s\"", _option_string.c_str());
-#endif
 	}
 
 	std::string _option_string = "-initxzeros 0";

NumLib/Function/Interpolation.h

@@ -79,7 +79,7 @@ void shapeFunctionInterpolate(
     auto const num_nodes = shape_matrix_N.size();
 
     assert(num_nodes*num_nodal_dof == nodal_values.size());
-    (void) num_nodal_dof; // no warnings when not in debug build
+    (void) num_nodal_dof; (void) num_nodes; // no warnings when not in debug build
 
     detail::shapeFunctionInterpolate<0>(nodal_values, shape_matrix_N, interpolated_value,
                                         interpolated_values...);

NumLib/ODESolver/NonlinearSolver-impl.h

@@ -47,7 +47,8 @@ solve(Vector &x)
 
     BLAS::copy(x, x_new); // set initial guess, TODO save the copy
 
-    for (unsigned iteration=1; iteration<_maxiter; ++iteration)
+    unsigned iteration=1;
+    for (; iteration<=_maxiter; ++iteration)
     {
         sys.preIteration(iteration, x);
 
@@ -65,7 +66,7 @@ solve(Vector &x)
 
         if (!iteration_succeeded)
         {
-            ERR("The linear solver failed.");
+            ERR("Picard: The linear solver failed.");
         }
         else
         {
@@ -76,12 +77,15 @@ solve(Vector &x)
                 // Although currently it is not.
                 break;
             case IterationResult::FAILURE:
+                ERR("Picard: The postIteration() hook reported a non-recoverable error.");
                 iteration_succeeded = false;
                 // Copy new solution to x.
                 // Thereby the failed solution can be used by the caller for debugging purposes.
                 BLAS::copy(x_new, x);
                 break;
             case IterationResult::REPEAT_ITERATION:
+                INFO("Picard: The postIteration() hook decided that this iteration"
+                     " has to be repeated.");
                 continue; // That throws the iteration result away.
             }
         }
@@ -92,26 +96,33 @@ solve(Vector &x)
             break;
         }
 
+        auto const norm_x = BLAS::norm2(x);
         // x is used as delta_x in order to compute the error.
         BLAS::aypx(x, -1.0, x_new); // x = _x_new - x
-        auto const error = BLAS::norm2(x);
-        // INFO("  picard iteration %u error: %e", iteration, error);
+        auto const error_dx = BLAS::norm2(x);
+        INFO("Picard: Iteration #%u |dx|=%.4e, |x|=%.4e, |dx|/|x|=%.4e,"
+             " tolerance(dx)=%.4e",
+             iteration, error_dx, norm_x, error_dx/norm_x, _tol);
 
         // Update x s.t. in the next iteration we will compute the right delta x
         BLAS::copy(x_new, x);
 
-        if (error < _tol) {
+        if (error_dx < _tol) {
             error_norms_met = true;
             break;
         }
 
         if (sys.isLinear()) {
-            // INFO("  picard linear system. not looping");
             error_norms_met = true;
             break;
         }
     }
 
+    if (iteration > _maxiter) {
+        ERR("Picard: Could not solve the given nonlinear system within %u iterations",
+            _maxiter);
+    }
+
     MathLib::GlobalMatrixProvider<Matrix>::provider.releaseMatrix(A);
     MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(rhs);
     MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(x_new);
@@ -153,7 +164,8 @@ solve(Vector &x)
     BLAS::copy(x, minus_delta_x);
     minus_delta_x.setZero();
 
-    for (unsigned iteration=1; iteration<_maxiter; ++iteration)
+    unsigned iteration=1;
+    for (; iteration<_maxiter; ++iteration)
     {
         sys.preIteration(iteration, x);
 
@@ -172,7 +184,7 @@ solve(Vector &x)
 
         if (!iteration_succeeded)
         {
-            ERR("The linear solver failed.");
+            ERR("Newton: The linear solver failed.");
         }
         else
         {
@@ -187,9 +199,12 @@ solve(Vector &x)
             case IterationResult::SUCCESS:
                 break;
             case IterationResult::FAILURE:
+                ERR("Newton: The postIteration() hook reported a non-recoverable error.");
                 iteration_succeeded = false;
                 break;
             case IterationResult::REPEAT_ITERATION:
+                INFO("Newton: The postIteration() hook decided that this iteration"
+                     " has to be repeated.");
                 // TODO introduce some onDestroy hook.
                 MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(x_new);
                 continue; // That throws the iteration result away.
@@ -208,7 +223,10 @@ solve(Vector &x)
         }
 
         auto const error_dx = BLAS::norm2(minus_delta_x);
-        DBUG("error of -delta_x %g and of residual %g,", error_dx, error_res);
+        auto const norm_x   = BLAS::norm2(x);
+        INFO("Newton: Iteration #%u |dx|=%.4e, |r|=%.4e, |x|=%.4e, |dx|/|x|=%.4e,"
+             " tolerance(dx)=%.4e",
+             iteration, error_dx, error_res, norm_x, error_dx/norm_x, _tol);
 
         if (error_dx < _tol) {
             error_norms_met = true;
@@ -216,12 +234,16 @@ solve(Vector &x)
         }
 
         if (sys.isLinear()) {
-            // INFO("  newton linear system. not looping");
             error_norms_met = true;
             break;
         }
     }
 
+    if (iteration > _maxiter) {
+        ERR("Newton: Could not solve the given nonlinear system within %u iterations",
+            _maxiter);
+    }
+
     MathLib::GlobalMatrixProvider<Matrix>::provider.releaseMatrix(J);
     MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(res);
     MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(minus_delta_x);

ProcessLib/Process.h

@@ -224,7 +224,7 @@ private:
 		_local_to_global_index_map.reset(
 		    new AssemblerLib::LocalToGlobalIndexMap(
 		        std::move(all_mesh_subsets),
-		        AssemblerLib::ComponentOrder::BY_COMPONENT));
+		        AssemblerLib::ComponentOrder::BY_LOCATION));
 	}
 
 	/// Sets the initial condition values in the solution vector x for a given