mjayadharan / mmmfe-st-dd Goto Github PK

I believe the test for deciding when to resize in GMRES in wrong (lines 1281 .... in darcy_vtdd.cc)
As written, the test is activated at the first two iterations (which is OK, but wasteful, and not the intended behavior). If the max number of iterations is chosen too small, the code will eventually crash at the last iteration, probably when doing

	Q_side[side][k_counter+1]=q[side];

at line 1445 (not sure).

There are actually two different issues here:

• the variable temp_array_size is declared as an int, but then the comparison
  temp_array_size<k_counter-2 is between an int and an unsigned int and may be wrong, so temp_array_size should be an unsigned int;
• furthermore, I believe the test should be temp_array_size<k_counter + 2 as otherwise the same line will crash (can you check this ?)

I wrote a small test program to reproduce the behavior. You can compile it as is, or with -D MODIFIED to see the fix I propose. You can see the results in the 2 enclosed files
out-mod.txt
out-orig.txt

#include <iostream>
using namespace std;

int main(int argc, char *argv[])
{
    unsigned int maxiter = 27;
#ifndef MODIFIED
    int temp_array_size = maxiter/4;
#else
    int temp_array_size = maxiter/4;
#endif    

    unsigned int k_counter = 0;
    std::cout << "Initial size " << temp_array_size << std::endl;
    while (k_counter < maxiter)
	  {
	      std::cout << "Doing iteration " << k_counter << std::endl;
#ifndef MODIFIED
	      if(temp_array_size<k_counter-2){
#else
	      if(temp_array_size < k_counter + 2){
#endif
		  std::cout << "\nResizing " << std::endl; 
		  temp_array_size*=2;
		  std::cout << "New size = " << temp_array_size << std::endl;
	      }
	      k_counter++;
	  }
    std::cout << "\nDone" << std::endl;

	  return 0;
}

I am confused by the computation of the GMRES residual, and the stopping test:

The variable e_all_iter is initialized as r_norm which is the actual norm of the initial residual (line 1273)

Then, during the iterations, one has (lines 1461 - 1465)

              combined_error_iter=fabs(Beta[k_counter+1])/r_norm;
              e_all_iter[k_counter+1]=(combined_error_iter);

so that e_all_iter is the decrease in the residual (but then the first value is wrong ?), and this is the value that gets printed.

Finally, the stopping test, and the value reported (lines 1471 - 1473) is combined_error_iter/e_all_iter[0], so there is one extra division by the initial residual.

I think what should be done is to set e_all_iter to the actual value of the residual, and when it is needed to divide by e_all_iter[0] (for deciding where to stop). Maybe also report the "target" residual when the iterations are started. Something like:
"initial residual is " << r_norm << " and stopping criteria is residual < " << tolerance * r_norm

When the pvtu files are created, the names of the corresponding .vtu files includes the directory name, and then the .pvtu and the .vtu files are written in the same directory.

If the .pvtu files are in the same directory as the .vtu files, the directory name should not be in the file name (another solution would be to write the .pvtu files in the parent directory, but I don't think that's a good idea).