diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h index 2ea6a21b716258097209ce17bf6045ccf2bd59c0..0fd8e68dd137a407e194e2da9e9d2a4e7c52ccda 100644 --- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h +++ b/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(); diff --git a/MathLib/LinAlg/BLAS.cpp b/MathLib/LinAlg/BLAS.cpp index 87550f22621bb9cf1493ab22f42e2a52eba642a6..51f2add2370cb4fd1169c593476234958b4013f8 100644 --- a/MathLib/LinAlg/BLAS.cpp +++ b/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 diff --git a/MathLib/LinAlg/BLAS.h b/MathLib/LinAlg/BLAS.h index a042b045446fb65f9bc7a48c3c5b896f8fa3cd94..5d4fae356e8ecb72fae0311a02c7bf84fc0e14e4 100644 --- a/MathLib/LinAlg/BLAS.h +++ b/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 diff --git a/MathLib/LinAlg/Lis/LisOption.h b/MathLib/LinAlg/Lis/LisOption.h index 74776c98bd7693e32087386dcb8bf8fec5a186d9..c67421c268a6fcbf003277a1ed9d1432854c609f 100644 --- a/MathLib/LinAlg/Lis/LisOption.h +++ b/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"; diff --git a/NumLib/Function/Interpolation.h b/NumLib/Function/Interpolation.h index b687665798bda6b8d2a765d0a72962176fe4e3dd..655e4a29bd0c2a8fba6105ceb8b4e3a016648464 100644 --- a/NumLib/Function/Interpolation.h +++ b/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...); diff --git a/NumLib/ODESolver/NonlinearSolver-impl.h b/NumLib/ODESolver/NonlinearSolver-impl.h index 6a5346f645f206f18125d2715bbaa74b25bc8ddf..5876ff11b79793b4b3c5b4b170d8b6f742c6aee8 100644 --- a/NumLib/ODESolver/NonlinearSolver-impl.h +++ b/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); diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h index 8e14d4cc9cda7eb8df156b58da2e73ba824bea6b..bda093cc0c37af58db4819599268498fab15c5a2 100644 --- a/ProcessLib/Process.h +++ b/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