GMRES not converging for Step-29

[David Fernández]

Hello all,

I am new to deal.ii. I am experimenting with step-29 example and I was trying to solve it using an iterative solver instead of the direct solver.  Because the system matrix is asymmetric I decided to try out GMRES. The solver did not converge, but since I just recently starting using deal.ii. I am assuming there could be a bug in my code.  I am only changing the solve method in the "UltrasoundProblem" class as shown bellow, leaving the rest of the file unchanged.

Thanks in advance for the advice.

Best regards,

David

------------------------------------------------------------------------------------------------------
 template <int dim>
  void UltrasoundProblem<dim>::solve ()
  {
    deallog << "Solving linear system... ";
    Timer timer;
    timer.start ();
    SolverControl control(1000, 1.e-4); //, 1.e-10, true, true );
    PreconditionSOR<SparseMatrix<double> > preconditioner;
    preconditioner.initialize (system_matrix);
    SolverGMRES<Vector <double>> solver_gmres (control );
    solver_gmres.solve(system_matrix, solution, system_rhs, preconditioner);

       std::cout << ">>> Using iterative solver!!!" << std::endl;

    timer.stop ();
    deallog << "done ("
            << timer ()
            << "s)"
            << std::endl;
  }

------------------------------------------------------------------------------------------------------


OUPUT:

DEAL::Generating grid... done (0.0722460s)
DEAL::  Number of active cells:  25600
DEAL::Setting up system... done (0.0713560s)
DEAL::  Number of degrees of freedom: 51842
DEAL::Assembling system matrix... done (0.307943s)
DEAL:GMRES::Solving linear system... Starting value 10.7267
DEAL:GMRES::Failure step 1000 value 1.73422


----------------------------------------------------
Exception on processing:

--------------------------------------------------------
An error occurred in line <828> of file </usr/local/deal.II/deal.II_8.1_mt/include/deal.II/lac/solver_gmres.h> in function
    void dealii::SolverGMRES<VECTOR>::solve(const MATRIX&, VECTOR&, const VECTOR&, const PRECONDITIONER&) [with MATRIX = dealii::SparseMatrix<double>; PRECONDITIONER = dealii::PreconditionSOR<dealii::SparseMatrix<double> >; VECTOR = dealii::Vector<double>]
The violated condition was:
    false
The name and call sequence of the exception was:
    SolverControl::NoConvergence (this->control().last_step(), this->control().last_value())
Additional Information:
Iterative method reported convergence failure in step 1000 with residual 1.73422
--------------------------------------------------------

Aborting!
----------------------------------------------------

[Uwe Köcher]

David,

your code is okay, but the parameters lead the solver to abort...

I've tried it with the similar code (please see this as debug code, some of the parameters are poor

and the final solution is not precise!):
   

    PreconditionSOR<SparseMatrix<double>> preconditioner;
    preconditioner.initialize(system_matrix);
    
    SolverControl sc(20000, 1e0, true, true);

    SolverGMRES<Vector<double>>::AdditionalData additional_data;
    additional_data.max_n_tmp_vectors = 1000;
    additional_data.force_re_orthogonalization = true;
    SolverGMRES<Vector<double>> A_direct(sc,additional_data);
    A_direct.solve(system_matrix, solution, system_rhs, preconditioner);

I think you have to find out, what is a good strategy for the numerical solution with GMRES for

such a problem.

You have a lot of possibilities for better preconditioners and the parameter values for GMRES,

which bring the performance to the solver...

[David Fernández]

Thank you Uwe Köcher for you reply.

I tried today playing around with this problem and setting up a GMRES solver tweaking the vector space size and the number of iterations with a realistic tolerance.
Also, I tried experimenting with BiCGStab.
In both cases I tried using the Identity, Diagonal (Jacobi), SOR and SparsILU preconditioners to no avail. The way I use the SparseILU preconditioner is the following:

    SparseILU<double> preconditioner_ilu;
    preconditioner_ilu.initialize (system_matrix, SparseILU<double>::AdditionalData());

So I printed the system (system_matrix, system_rhs and the solution obtained with the direct solver) and read the files from MatLab. With this, I was able to solve the system using both BiCGStab and GMRES from MatLab with an ILU preconditioner with these settings:

      [L,U] = ilu(A,struct('type','ilutp','droptol',1e-2));
      [ x_iter,flag,relres,iter,resvec] = gmres(A,b,[],1e-9,100, L, U);

I am wondering, has anyone been able to solve this problem with an iterative solver from deal.ii?
I would like to have a solution using an iterative solver for this problem.

Regards,

[Bruno Turcksin]

David,

first, solvers in deal.II uses an absolute tolerance, I would recommend that you use a relative tolerance instead something like:

SolverControl control(1000,1e-6*system_rhs.l2_norm());

I was able to solve the problem by increasing the vector space size and using:

SparseILU<double>::AdditionalData additional_data(0,100);
SparseILU<double> preconditioner_ilu;
preconditioner_ilu.initialize (system_matrix, additional_data);

The ILU that you used with Matlab and the ILU in deal.II do not work the same way. With Matlab, you drop the terms in the matrix that are small and then, you do an exact LU decomposition. In deal.II, you control the size of the inverse of the matrix ("how full" do you allow the inverse to be). 

Best,

Bruno

[David Fernández]

Thanks Bruno that did the trick!!!

[Guido Kanschat]

...and instead of computing a relative tolerance in this complicated way, you can use ReductionControl.

Best,

Guido

[Jason McIntosh]

I'm trying to do the same thing David was trying to do over a year ago, but I'm not figuring it out.  I tried to follow this thread and changed the original step29 code for  UltrasoundProblem<dim>::solve() from

        SparseDirectUMFPACK  A_direct;
        A_direct.initialize(system_matrix);
        A_direct.vmult(solution, system_rhs);

to

        SolverControl      solver_control(1000, 1e-6*system_rhs.l2_norm());
        SparseILU<double>::AdditionalData additional_data(0, 100);
        SparseILU<double> preconditioner_ilu;
        preconditioner_ilu.initialize(system_matrix, additional_data);
        SolverGMRES<Vector<double>> solver(solver_control);
        solver.solve(system_matrix, solution, system_rhs, preconditioner_ilu);


but it's not converging.  Suggestions?

Jason

[Bruno Turcksin]

Jason,

you can increase the maximum number of iterations (change the 1000). The best way to speed up the convergence of GMRES is to do less restart but this will require more memory. To do that you can do (https://dealii.org/developer/doxygen/deal.II/structSolverGMRES_1_1AdditionalData.html):
unsigned int restart = 35;
SolverGMRES< Vector<double> >::AdditionalData gmres_additional_data(restart+2);
SolverGMRES< Vector<double> > solver(solver_control, gmres_additional_data);

Best,

Bruno

[Jason McIntosh]

Hi Bruno,

Thanks for the reply, but it's still not converging and unfortunatley I don't know these classes or algorithms well enough to know what I should try next. 

I added the gmres_additional_data object (see code below).  Prior to adding it, the code reported a "Starting value" of

1.0677e+24

before, and aborted with a value of

2.949e+16

after 1000 Iterations.  With the addition of gmres_additional_data, it aborts at step 1000 with a value of

4.416e+15

Maybe you could call that an improvement, but I think that the iteration is stuck bouncing around some point.

When I included the gmres_additional_data class, the solver reported:

DEAL:GMRES::Re-orthogonalization enabled at step 30

which doesn't happen if gmres_additional_data isn't used. 

Any more suggestions to try?

thanks, Jason


       double cc = 1e-6*system_rhs.l2_norm();
        std::cout << "convergence condition = " << cc << std::endl;

        SolverControl      solver_control(1000, cc);

        SparseILU<double>::AdditionalData additional_data(0, 100);
        SparseILU<double> preconditioner_ilu;
        preconditioner_ilu.initialize(system_matrix, additional_data);

        unsigned int restart = 35;
        SolverGMRES< Vector<double> >::AdditionalData gmres_additional_data(restart + 2);
        SolverGMRES<Vector<double>> solver(solver_control, gmres_additional_data);

        solver.solve(system_matrix, solution, system_rhs, preconditioner_ilu);

[Wolfgang Bangerth]

On 03/31/2016 02:11 PM, Jason McIntosh wrote:
>
> I added the gmres_additional_data object (see code below). Prior to
> adding it, the code reported a "Starting value" of
>
> 1.0677e+24
>
> before, and aborted with a value of
>
> 2.949e+16
>
> after 1000 Iterations. With the addition of gmres_additional_data, it
> aborts at step 1000 with a value of
>
> 4.416e+15

Both of these numbers are essentially zero compared to the starting
residual. You have successfully solved your linear system, assuming that
you started with a zero vector and the 1.06e24 reflects only the norm of
the right hand side vector.

What is the value of your cc variable in the code snippet you show?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Jason McIntosh]

The convergence condition (cc in the code) is 288300.  Quite a bit smaller than the e15 values that we're converging to.

I compiled the code in Windows using VisualStudio.  Maybe that's the problem.  Timo gave me a link to a linux based VirtualBox image with deal II installed.  I'll get that installed and try the code with that.

Jason

[Wolfgang Bangerth]

On 03/31/2016 04:34 PM, Jason McIntosh wrote:
> The convergence condition (cc in the code) is 288300. Quite a bit smaller
> than the e15 values that we're converging to.

If you start with a residual of 1e24, you can't expect to converge down to a
residual of 2.8e5. Think about the fact that you will never get a relative
accuracy less than 1e-16 from floating point arithmetic, which would be around
1e8 in the absolute best case.

[Jason McIntosh]

>If you start with a residual of 1e24, you can't expect to converge down to a residual of 2.8e5.

Yes, that makes sense.

Unfortunately the step-29 code is not behaving the same using the VisualStudio compiler and the linux compiler.  I installed a Ubuntu image that Timo pointed me to.  When I use the UMFPACK solver, I get a good result from the linux code.  When I use the GMRES solver with linux, it converges, but when I plot the results using VisIt, the resulting solution is obviously bad.  This is the output from the linux code.  Apparently it's a debug build:


[100%] Run step-29 with Debug configuration
DEAL::Generating grid... done (0.117687s)

DEAL::  Number of active cells:  25600

DEAL::Setting up system... done (0.0688160s)

DEAL::  Number of degrees of freedom: 51842

DEAL::Assembling system matrix... done (0.470934s)
convergence condition = 203572
DEAL:GMRES::Solving linear system... Starting value 2.36079e+09
DEAL:GMRES::Convergence step 23 value 94124.6
DEAL::done (2.86510s)
DEAL::Generating output... done (0.410455s)
[100%] Built target run

This is the output from the Windows code.  This is a release build, but the step-29.cc code is identical to the linux code:

DEAL::Generating grid... done (1.24801s)
DEAL::  Number of active cells:  102400
DEAL::Setting up system... done (0.936006s)
DEAL::  Number of degrees of freedom: 206082
DEAL::Assembling system matrix... done (1.87201s)
convergence condition = 288300
DEAL:GMRES::Solving linear system... Starting value 1.06777e+24

DEAL:GMRES::Re-orthogonalization enabled at step 30

DEAL:GMRES::Failure step 1000 value 4.41606e+15


One extra complication here is that the image Timo gave me had deal II 8.3 installed, while I'm using 8.4 with Windows. 

It looks like I need to get UMFPACK working in Windows and hope it solves the problem correctly.

Jason

[Bruno Turcksin]

Jason,

there is one more thing that you can try. SparseILU does an inexact LU
decomposition because it throws away some entries in the matrix to
keep the memory requirement low. If you increase the number of
off-diagonal elements allowed, you will get a better preconditioner
and in the limit where the matrix can be full, you should get an exact
LU decomposition. So replace SparseILU<double>::AdditionalData
additional_data(0, 100); with SparseILU<double>::AdditionalData
additional_data(0, 500); and see if it works better.


Best,


Bruno

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/-NflTyZpOd4/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Jason McIntosh]

Thanks Bruno.  I'll give it a try and get back.

I wanted to give an update to my previous post that indicated that the identical code was giving different results between the 8.3 linux build and the 8.4 windows build.  I checked the step-29.prm parameter file and found that the 8.3 code on Timo's image has slightly different parameter values than what I got with the 8.4 code.  Using the identical parameter file on both builds produced the exact same numerical results with GMRES.  They didn't  converge, but at least I know that the code is behaving the same.  Big relief for me.

Jason

[Jason McIntosh]

Using the same code and parameter files, I'm getting the same results between linux and windows.  However, I still can't get the GMRES implementation to produce a good result.  I can get it to reduce it's error by 7 or 8 orders of magnitude, but the resulting intensity plots look bad.  When I use the UMFPACK solver in linux, I get a nice output that looks like the step-29 tutorial documentation.

My ultimate problem is going to look a lot like step-29, but I can't get a good solution to run just using the deal ii solvers.  It seems that UMFPACK (which requires LAPACK and BLAS) will have to be used.  I now have to embark down the path of getting these libraries to compile with the VisualStudio compiler.  Ultimately speed will be a concern for me, so I'm assuming that UMFPACK and an optimized BLAS library will be the best way to go anyway. 

Jason

[Jason McIntosh]

I can't get the GMRES solver to converge.  The attached image "GMRES solution.jpg" shows what I get with GMRES using the code below, and "UMFPACK solution.jpg" shows what the original step-29 code using the UMFPACK solver generates.

 

                                SolverControl      solver_control(1000, convergence_condition);

 

                                SparseILU<double>::AdditionalData additional_data(0, 100);

                                SparseILU<double> preconditioner_ilu;

                                preconditioner_ilu.initialize(system_matrix, additional_data);

 

                                unsigned int restart = 35;

                                SolverGMRES< Vector<double> >::AdditionalData gmres_additional_data(restart + 2);

                                SolverGMRES<Vector<double>> solver(solver_control, gmres_additional_data);

 

                                solver.solve(system_matrix, solution, system_rhs, preconditioner_ilu);

 

The GMRES solver does reduce the error by over 6 orders of magnitude  (1.88354e+09/371.922 = 5e6) but that's as good as I got with the set of parameters I was working with.  (Note that I changed the frequency in the *.prm file from 3e7 to 2e7 to make the solution a bit smoother.  This helped GMRES to at least generate a solution that didn't look like numerical noise.)

 

                DEAL::  Number of degrees of freedom: 206082

                DEAL::Assembling system matrix... done (1.87201s)

                convergence condition = 400

                DEAL:GMRES::Solving linear system... Starting value 1.88354e+09

                DEAL:GMRES::Convergence step 17 value 371.922

                DEAL::done (59.2804s)

                output_file = solution

                DEAL::Generating output... done (28.3922s)

 

 

So far what you've seen is the solutions with the parameter file having the variable "Number of refinements" set to 6, which is what I got with the 8.3 code.  But the 8.3 code I got from Timo's VirtualBox image had this set to 5.  When I set it to 5 and run it with the same GMRES code as shown above, I can reduce the error from 7.44224e+14 to 40.  Pretty good, but the solution is mostly noise (see "GMRES solution with 5 refinement steps.jpg").

 

                DEAL:GMRES::Solving linear system... Starting value 7.44224e+14

                DEAL:GMRES::Re-orthogonalization enabled at step 15

                DEAL:GMRES::Convergence step 210 value 39.9634

 

 

To try to understand what could be going on I  exported the sparsity pattern and plotted it in "sparsity pattern.jpg".  It seems well formed, at least as far as I can tell.  I looked at what was going on at about the center of the matrix and decided to see how many non-zero entries there were in the matrix and how "wide" the triangular spread is.  I generated the data below for the 6 refinement steps case.  Note that the size of the matrix is 206082x206082. 

 

All rows have 18 non-zero entries.  Some of the rows are quite tight (just 25 indexes separating the non-zero entries) , while some are quite wide (over 35k entries separating the non-zero entries.  Does this mean anything?

 

At this point I'm lost as to what I should do and I need some suggestions.  I'll attach the code and parameter file I'm using.  I believe I've tried all of the suggestions that have been given to me so far, but if anyone wants to try to outline a range of parameters or approaches to try, I'll do so.

 

Thanks for reading, Jason

 

1st column = row #  

2nd column = # of non-zero entries in row

3rd column = max non-zero column index - min non-zero column index

                      (this is basically how wide the sparse matrix is)

100000      18      73

100001      18      73

100002      18   35657

100003      18   35657

100004      18   35641

100005      18   35641

100006      18      41

100007      18      41

100008      18      25

100009      18      25

100010      18   35641

100011      18   35641

100012      18   37001

100013      18   37001

100014      18      25

100015      18      25

100016      18    1385

100017      18    1385

100018      18      41

100019      18      41

100020      18      25

100021      18      25

100022      18      73

100023      18      73

100024      18      57

100025      18      57

100026      18      25

100027      18      25

100028      18    1385

100029      18    1385

100030      18      57

100031      18      57

100032      18    1417

100033      18    1417

100034      18     137

100035      18     137

100036      18      57

100037      18      57

100038      18     105

100039      18     105

100040      18      25

[Jason McIntosh]

PS:

So I guess what I need is for someone to get step-29 to work with one of the deal II iterative solvers, or point me in a direction to get it solving myself.  I don't know what to try at this point.

[Bruno Turcksin]

Jason,

you have three possibilities:
- use a physics-based preconditioner
- use an algebraic multigrid preconditioner. This is the most general
approach but you will need Trilinos or PETSc and I am not sure that
these libraries support windows
- use a geometric multigrid preconditioner. Timo and Martin know a
lot more than me about that. Any input on this, guys ?

Best,

Bruno

[Timo Heister]

I would first try to run GMRES with something like ILU on a very
coarse mesh. This should work and give correct results. You only need
to think about better preconditioners if the number of iterations
needed explodes for h->0. If you already have issues on a coarse mesh,
the issue is something else: maybe the system is not invertible,
scaling of the equations is inconsistent, etc..

--
Timo Heister
http://www.math.clemson.edu/~heister/

[Jason McIntosh]

Thanks Bruno.  Unfortunately I don't know how to make a better perconditioner.  I can look into it though.  Any guidance you can suggest would be appreciated.  (I'm guessing that I'll be purchasing some books...  I want to learn this stuff, but I've got a finite amount of time I can devote to it.)

One thing that's bothering me is the GMRES starting values.  You originally suggested to David F. to use a "relative tolerance" of 1e-6*system_rhs.l2_norm().  However, the GMRES algorithm is printing out error values that aren't tracking the system_rhs.l2_norm() value.  Below is a little table I put together. 

There's a strong dependence on the extra data passed to SparseILU<double>::AdditionalData, the frequency, and the # of refinements.  But none of these changes are strongly effecting system_rhs.l2_norm().  Does this indicate some other hidden problem?

Jason


# of refinements    system_rhs.l2_norm()    GMRES starting value    extra data     omega (freq)
6                               2.98e11                                    1.88e9                          100               2e7
5                                2.07e11                                    7.44e14                       100              2e7
5                                2.07e11                                    2.09e9                         200               2e7
5                                2.03e11                                    7.50e6                         200                3e7
5                                2.03e11                                    2.36e9                         100                3e7
5                                2.03e11                                    5.080e35                        0                3e7

I'll paste in the solver code again for your reference:

        SolverControl      solver_control(1000, convergence_condition);

        SparseILU<double>::AdditionalData additional_data(0, 200*0);

[Bruno Turcksin]

Jason,

2016-04-04 16:55 GMT-04:00 Jason McIntosh <jason...@gmail.com>:
> There's a strong dependence on the extra data passed to
> SparseILU<double>::AdditionalData, the frequency, and the # of refinements.
> But none of these changes are strongly effecting system_rhs.l2_norm(). Does
> this indicate some other hidden problem?

This is perfectly normal. The system rhs depends on the equation, the
mesh, and the finite elements that you are using. Refining the mesh
has a very strong effect to the norm of the rhs. The reason I said to
use a relative tolerance is the following: image that you are
computing the position of the moon relative to the Earth in meters,
all your numbers are going to be huge and it wouldn't make sense to
put the tolerance at 1e-9. The tolerance should be adapted to the
order of magnitude of your problem. The starting value on the other
hand also depends on the preconditioner that you are using that's why
when you get a better ILU, this number decrease.

Best,

Bruno

[Jason McIntosh]

>I would first try to run GMRES with something like ILU on a very
>coarse mesh. This should work and give correct results. You only need
 

I ran the UMF and GMF code for different #'s of refinements.  The table below gives the size of the system for each refinement level and the GMRES starting and ending error values.  It seems that GMRES is doing a good job of driving the error down, but the solution only looked like the UMF solution for the 3 refinement case.  I've attached images of what the solution looked like 3,4, and 5 refinements. 

UMF gets a good looking solution for all levels of refinement, but it's way off for the smallest system (3 refinements). 

While GMRES gets a solution that looks a lot like UFM for the 3 refinement case, it generates noisy junk for higher refinement levels.

Now what?

Jason

extra data = 100
omega=3e7

# of refinements       # of active cells           #DOF                GMRES start        GMRES stop
3                                    1600                        3362                1.50e0                     1.21e-31
4                                    6400                      13122                4.05e12                    2.32e3
5                                  25600                      51242                2.36e9                      1.95e3
6                                 102400                    206082               1.30e16                    3.10e7

[Jason McIntosh]

I don't like the idea of separating the real and imaginary parts of this equation.  Just makes me feel uneasy to have two equations that are VERY similar to one another, but yet completely independent.  So I had the idea of adding some physical damping to the system.  This can be done by making the speed of sound complex, or

c = co - i *α

where co is the "normal" speed of sound and  α is an attenuation, typically α << co.  (Sorry for the bad equation formatting, I don't know how to insert an equation into this editor.)  So the standard wave equation changes from

−ω2U−c2ΔU=0

to

−ω2U−(co − i∗α)2ΔU=0

when U is decomposed into it's real and imaginary parts v=real(U) w=imag(U) we get

−ω2v−(co2−α2)Δv − 2coαΔw=0

−ω2w−(co2−α2)Δw+2coαΔv=0

So now we've got a loosely coupled set of equations instead of equations that
 are only coupled through the boundary condition.

My problem now is how do I modify the step-29 code to take this coupling into account.

I somewhat blindly change the code to that shown below, but it hasn't made 
any difference with the solver.  Did I do it right?


double alpha = c / 100.0; // attenuation embedded into a complex speed of sound
					        // speed of sound = c-i*alpha

if (fe.system_to_component_index(i).first ==
	fe.system_to_component_index(j).first)
{
// this part is unchanged
	for (unsigned int q_point = 0; q_point<n_q_points; ++q_point)
		cell_matrix(i, j) += (((fe_values.shape_value(i, q_point) *
			fe_values.shape_value(j, q_point)) *
			(-omega * omega)
			+
			(fe_values.shape_grad(i, q_point) *
				fe_values.shape_grad(j, q_point)) *
			(c * c - alpha * alpha)) *
			fe_values.JxW(q_point));

}
else
{
// this is the part i added... am I doing something wrong???
	double coeff = 2*c*alpha;
	if (fe.system_to_component_index(i).first == 0) // real part?
		coeff = -coeff;

	unsigned int component_i = fe.system_to_component_index(i).first;
	unsigned int component_j = fe.system_to_component_index(j).first;
						
	for (unsigned int q_point = 0; q_point<n_q_points; ++q_point)
		cell_matrix(i, j) += 
			(fe_values.shape_grad(i, q_point)[component_i] *
				fe_values.shape_grad(j, q_point)[component_j]) *
			coeff *
			fe_values.JxW(q_point);
}

[Jason McIntosh]

Putting equations into that last post didn't seem to work out very well.  I hope this is readable.

I'm still struggling with getting step-29 to converge using the iterative solver.  I thought I'd smooth the solution and couple the real and imaginary parts of the problem together by adding an attenuation term to the fundamental equations.  I did this by making the speed of sound complex:

c=co - i*alpha

when this is done, the wave equation:

-w^2 * U  - (co - i*alpha)^2  DELTA U = 0

breaks into the real and imaginary equations

-w^2 * V  - (co^2 - alpha^2) DELTA V    -  2*co*alpha* DELTA W = 0
-w^2 * W  - (co^2 - alpha^2) DELTA W  + 2*co*alpha* DELTA V = 0

where w is the radian frequency, V=real(U) and W=imag(U).

Now the V and W equations are coupled.  But I don't know how to add this to the step-29 program.  What I tried is pasted in below.  If someone could tell me if it's right or correct it, that would be great.

Jason

double alpha = c / 100.0; // attention embedded into a complex speed of sound

                                         // speed of sound = c-i*alpha

if (fe.system_to_component_index(i).first ==
    fe.system_to_component_index(j).first)
{

// this is the old code

    for (unsigned int q_point = 0; q_point<n_q_points; ++q_point)
        cell_matrix(i, j) += (((fe_values.shape_value(i, q_point) *
            fe_values.shape_value(j, q_point)) *
            (-omega * omega)
            +
            (fe_values.shape_grad(i, q_point) *
                fe_values.shape_grad(j, q_point)) *
            (c * c - alpha * alpha)) *
            fe_values.JxW(q_point));

}
else
{

// and this is what I added

[Wolfgang Bangerth]

On 04/04/2016 05:34 PM, Jason McIntosh wrote:
> >I would first try to run GMRES with something like ILU on a very
> >coarse mesh. This should work and give correct results. You only need
>
> I ran the UMF and GMF code for different #'s of refinements. The table below
> gives the size of the system for each refinement level and the GMRES starting
> and ending error values. It seems that GMRES is doing a good job of driving
> the error down, but the solution only looked like the UMF solution for the 3
> refinement case. I've attached images of what the solution looked like 3,4,
> and 5 refinements.
>
> UMF gets a good looking solution for all levels of refinement, but it's way
> off for the smallest system (3 refinements).
>
> While GMRES gets a solution that looks a lot like UFM for the 3 refinement
> case, it generates noisy junk for higher refinement levels.
>
> Now what?

I would start by trying to make the problem simpler. What happens if you
remove the coupling between the real and imaginary part by making all parts of
the boundary Dirichlet, for example?

In this case, you end up with a block matrix of the form
[ A 0 ]
[ 0 A ]
for which the perfect preconditioner is going to be
[ A^-1 0 ]
[ 0 A^-1 ]
In fact, this preconditioner would also be a pretty good choice for the
coupled problem
[ A B ]
[ -B A ]
where B is induced by the boundary integrals.

The question is how you implement the preconditioner above. To understand how
to do something like this, you'll need to understand how we build block
matrices and block preconditioners. Take, for example, a look at step-22 (and
in particular the "Results" section). Or, if you're using the newest deal.II
version, you could use the facilities of the linear operator classes:
https://www.dealii.org/8.4.0/doxygen/deal.II/group__LAOperators.html
Combine this with a BlockLinearOperator to make a preconditioner out of the
individual blocks.

What this all means is that you only need to be able to solve with the matrix
A, for which you'll have to have a preconditioner for A. This should be
substantially simpler than finding a preconditioner for the entire,
non-symmetric matrix.

As an aside, my take on building preconditioners for coupled problems is in
lecture 38 here:
http://www.math.tamu.edu/~bangerth/videos.html

[Jason McIntosh]

Thanks Wolfgang.  I'll look into those links and try to understand the system better. 

If purely real Dirichlet boundary conditions are used, then the solution will be purely real.  I suppose it would be worthwhile to see if that kind of system would converge.

I have a lot of concern about a complex system converging with an iterative solution, when the real and imaginary parts of the equations have been separated.  When I think of the equations that are being solved, I see a  "real system" and a completely independent "imaginary system" trying to be solved.  So the iterator is going to be trying to strongly drive the solution for these two decoupled systems.  However, the boundary conditions will be coupling these two systems, and so these relatively few boundary condition equations will be propagating information up into the decoupled real and imaginary systems.  This propagation will be a relatively "weak force" on the numerical solution, partly because their relatively few equations contributing a relative small contribution to the error, but largely because their influences will be slowly propagated up through the large A systems.  Worse, when the real or imaginary system change, the other (imaginary or real) system must change with it.  But again there's no strong coupling in the equations to get them to move together.  Their change must be coupled through the boundary conditions and then propagated back through the other system. 

[ A 0 ]    real system
[ 0 A ]    imaginary system
[ B B ]    boundary conditions

So in my mind, if the A system is on the order of 100,000 equations, then there must be something on the order of 100,000 iterations for a change in the boundary conditions to fully propagate through the A system.  This just doesn't seem like a system that's ever going to converge, so I would expect that a completely complex system where the real and imaginary parts are combined into a single A system would have a much better chance of converging. 

I have to apologize for ignorantly handwaving my way through a problem that you're an expert in.  I'm doing my best to understand the problem.

Do you see any merit to trying to implement a complex system matrix and complex solver?  Do these iterative solver approaches work for complex systems?  If that's the best solution, I'm willing to put the effort into coding it.

Jason

[Wolfgang Bangerth]

On 04/07/2016 11:22 AM, Jason McIntosh wrote:
>
> If purely real Dirichlet boundary conditions are used, then the solution
> will be purely real. I suppose it would be worthwhile to see if that
> kind of system would converge.

Think of it this way: if you can't solve this simpler system
successfully, the more complex coupled problem won't be solvable using
this scheme either. So start simple.

But that is not how Krylov space methods work. Even without
preconditioning, information propagates from every node to all of its
neighbors in each iteration just through the matrix-vector product, and
then they do global scalings that affect all other nodes as well. With
most preconditioners, information propagates through the entire domain
in each iteration.

> Do you see any merit to trying to implement a complex system matrix and
> complex solver?

No. Any solver that works on the complex system should, with suitable
modifications, also work on the coupled problem of separate real and
imaginary parts; and vice versa.

There is nothing wrong with implementing such a solver. But you still
need a preconditioner, and that seems to me what you don't currently
have here.

[Jason McIntosh]

I finally got the GMRES solver to converge to a good solution for step-29!

 

 I'll paste in the parameters and the code that I used below.  What finally got it to converge was to set the max_n_tmp_vectors of SolverGMRES::AdditionalData to 500 instead of 35.  Someone had recommended using this parameter, but had suggested a value of 35.  I got frustrated and cranked it up to 500 and I finally GMRES to converge to a good solution.  I don't know what this does yet, but obviously I need to look into it more.

 

Thanks for everyone's help!

 

Next question: 

The parameters below are requiring around 350 iteration steps.  This is going to be too slow for most of my problems.  However, my solutions will be run over a frequency range, usually about 200 frequencies, with the solution for the previous frequency being very similar to the solution for the next frequency.  Is it possible to use the solution from a previous frequency to help speed the system convergence for the next frequency?  Note that once all of the frequency dependent data is put into the problem, the dependence on omega w is going to be more complex than just the w^2 that shows up in the Helmholtz problem.

 

Jason

 

 

Focal distance        = 0.3

 Number of refinements = 5

 c     = 1.5e5

omega = 3.0e7

convergence_condition = .01

nmax_iterations= 1000

additional_diagional_terms = 200

max_n_tmp_vectors = 500

 

SolverControl      solver_control(nmax_iterations, convergence_condition);

SparseILU<double>::AdditionalData additional_data(0, additional_diagional_terms);

SparseILU<double> preconditioner_ilu;

preconditioner_ilu.initialize(system_matrix, additional_data);

SolverGMRES< Vector<double> >::AdditionalData gmres_additional_data(max_n_tmp_vectors + 2);

[Bruno Turcksin]

Jason,

2016-04-08 13:13 GMT-04:00 Jason McIntosh <jason...@gmail.com>:
> I'll paste in the parameters and the code that I used below. What finally
> got it to converge was to set the max_n_tmp_vectors of
> SolverGMRES::AdditionalData to 500 instead of 35. Someone had recommended
> using this parameter, but had suggested a value of 35. I got frustrated and
> cranked it up to 500 and I finally GMRES to converge to a good solution. I
> don't know what this does yet, but obviously I need to look into it more.

This is the number of iterations that GMRES is doing before it
restarts. The bigger the number the better GMRES converges but the
more memory you will need (in exact exact arithmetic, GMRES is
guaranteed to converge to the exact solution if you don't restart it).
With this number to 500, you keep up to 502 vectors with the size of
your solution.

> The parameters below are requiring around 350 iteration steps. This is
> going to be too slow for most of my problems. However, my solutions will be
> run over a frequency range, usually about 200 frequencies, with the solution
> for the previous frequency being very similar to the solution for the next
> frequency. Is it possible to use the solution from a previous frequency to
> help speed the system convergence for the next frequency?

When you call solver.solve(system_matrix, solution, system_rhs,
preconditioner_ilu); solution is used as an initial guess. If the
solution doesn't change too much, GMRES should converge very quickly.

Best,

Bruno