Commit 8ae27fe5 authored by Christopher Subich's avatar Christopher Subich

Simplify GMRES dot products; fix convergence errors

This commit amends Cheb_2dmg::resid_dot (dot-product of
residual vectors in the context of the 2D GMRES/multigrid
iteration) to perform the dot-product in the simplest way,
as a sum over grid elements.

Before this, the code (which is left in the file but
disabled by #if 1 / #else) tried to be clever and
evenly-weight contributions from the interior, the boundary
conditions, and the mean-pressure error (if the problem was
indeterminate).  This evidently caused the long-standing
problem of the GMRES inner iteration being far more
optimistic about error (as reported by LAPACK) than would be
computed by directly measuring the residual error.

This problem was ultimately discovered when the
doc_map_iwave case broke.  Even though this broken residual
code has been in SPINS since before it was
version-controlled via GIT, the map_iwave case would still
converge (and relatively quickly) when the multigrid
iteration had a relatively large coarse-solve size.

Git bisect pinned the broken documentation case on the
change that reduced the coase-grid size, which reduced the
GMRES convergence rates enough that the errors caused by
*this* problem broke convergence.

Since this change replaces the error norm calculations used
by GMRES, it should be treated with caution and mapped-grid
cases should have a mark-1-eyeball inspection of results for
a while.

--

Additionally, for debugging utility this change also
includes crude residual/basis output for the 2D GMRES
solver.  All current uses are commented out.
parent e6a00da4
......@@ -74,6 +74,10 @@ template <class resid_type, class basis_type> class GMRES_Interface {
virtual void resid_scale(resid_type & a, double c) = 0;
virtual void basis_scale(basis_type & a, double c) = 0;
/* For debugging purposes, resid and basis outputs */
virtual void resid_write(resid_type & a, int seq) {};
virtual void basis_write(basis_type & a, int seq) {};
virtual bool noisy() const { return false; };
};
......@@ -278,6 +282,8 @@ template <class Controller> class GMRES_Solver {
/* And scale for assignment to the residual vector */
ops->resid_copy(resid_vec(0),start_r);
ops->resid_scale(resid_vec(0),1/starting_norm);
//ops->resid_write(resid_vec(0),0);
int my_it = 1;
......@@ -286,6 +292,7 @@ template <class Controller> class GMRES_Solver {
/* Find an x vector that approximately solves the prior basis */
timing_push("gmres_precond");
ops->precondition(resid_vec(my_it-1), basis_vec(my_it-1));
//ops->basis_write(basis_vec(my_it-1),my_it-1);
timing_pop();
timing_push("gmres_matmul");
/* Apply the operator the current basis vector */
......@@ -308,6 +315,9 @@ template <class Controller> class GMRES_Solver {
double norm = sqrt(ops->resid_dot(resid_vec(my_it),resid_vec(my_it)));
// fprintf(stderr,"norm %g\n",norm);
ops->resid_scale(resid_vec(my_it),1/norm);
//ops->resid_write(resid_vec(my_it),my_it);
//cout << resid_vec[my_it];
hess(my_it+1,my_it) = norm;
timing_pop();
......@@ -316,10 +326,12 @@ template <class Controller> class GMRES_Solver {
timing_push("gmres_hess");
{
/* LAPACK overwrites the matrix, so copy */
// std::cerr.precision(8);
// std::cerr.width(8);
// std::cerr.setf( std::ios::scientific );
// std::cerr << hess;
#if 0
std::cerr.precision(5);
std::cerr.width(5);
std::cerr.setf( std::ios::scientific );
std::cerr << hess;
#endif
hess_copy = hess;
int M = my_it+1, N = my_it;
rhs_vec = 0; rhs_vec(1) = 1;
......
......@@ -128,10 +128,23 @@ void Cheb_2dmg::free_basis(ubox * & f) {
}
double Cheb_2dmg::resid_dot(fbox * & l, fbox * & r) {
// double s = pssum(sum(*l->gridbox * *r->gridbox),my_comm);
#if 1
double s = pssum(sum(*l->gridbox * *r->gridbox),my_comm);
// return pssum(sum(*l->gridbox * *r->gridbox),my_comm);
// if (indefinite_problem)
// s = s + l->mean*r->mean/(szx*szz);
if (indefinite_problem)
s = s + l->mean*r->mean;
#else
/* The following code attempts to create a dot product for residuals that more-evenly balances
RHS error (from nabla^2 u = f), boundary condition error, and error attributable to the
mean-pressure condition (<p> = 0), if this is an indeterminate problem.
It also appears to be wrong. Actually enabling this code causes GMRES itself to become
very confused, where error measures reported by LAPACK within the GMRES inner iteration
no longer match a re-calculation from the original RHS and the estimated solution.
Whether this is due to some simple bug or a more subtle, mathematical error is
currently unknown (Feb 2018). Further investigation will wait for the simple calculation
(above) actually causing problems in use. */
double s = 0;
int imin, imax;
if (type_x == CHEBY) {
......@@ -173,6 +186,7 @@ double Cheb_2dmg::resid_dot(fbox * & l, fbox * & r) {
// indefinite problem
if (indefinite_problem)
s = s + l->mean*r->mean;
#endif
return pssum(s,my_comm);
}
......@@ -999,6 +1013,18 @@ void Cheb_2dmg::matrix_multiply(ubox * & lhs, fbox * & rhs) {
bool Cheb_2dmg::noisy() const {
return false && master(my_comm);
//return indefinite_problem && master(my_comm);
}
void Cheb_2dmg::resid_write(fbox * & vec,int seq) {
write_array(*(vec->gridbox),"resid",seq);
write_reader(*(vec->gridbox),"resid",true);
fprintf(stderr,"Residual %d mean: %.5e\n",seq,vec->mean);
}
void Cheb_2dmg::basis_write(ubox * & vec,int seq) {
write_array(*(vec->gridbox),"basis",seq);
write_reader(*(vec->gridbox),"basis",true);
fprintf(stderr,"Basis %d sigma: %.5e\n",seq,vec->sigma);
}
......
......@@ -132,6 +132,9 @@ class Cheb_2dmg : public GMRES_Interface < fbox * , ubox *> {
void resid_scale(fbox * & vec, double scale);
void basis_scale(ubox * & vec, double scale);
void resid_write(fbox * & vec, int seq);
void basis_write(ubox * & vec, int seq);
bool noisy() const;
private:
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment