search

fprimeau / EARTHSS220

Julia jupyter notebooks for UCI's EarthSS 220 Course.
4 stars 2 forks source link

DIP+DOP+POP+DO2 is too big #7

Open briochemc opened 5 years ago

briochemc commented 5 years ago

@fprimeau As we discussed might happen, I can't even factorize the matrix on my laptop. FWIW, this is the sparsity pattern, which looks OK (but, yes quite big, with 10'000'000 non-zeros...):

Screen Shot 2019-05-28 at 10 28 36 AM

I'll try to draft a block-newton algorithm and see if it works on this as a test.

briochemc commented 5 years ago

I'm making some progress but it's not easy...

Of course, another solution is to use a Newton-Krylov solver. So when I'm tired of trying the block-Newton approach, I'll try to see if it's easy to just make the current solver use a Krylov iterative solver for each backslash by adding a Krylov flag.

fprimeau commented 5 years ago

Hi Benoit, At the very least, we can separate the DO2 equation from the P-cycle equations and solve the two separate systems sequentially.

Do you know if the sparse factorization has the options for incomplete factorization available in Julia? iLU is an easy way to construct a preconditioner that sometimes works.

Best, François.

Sent from my iPhone

On May 27, 2019, at 7:43 PM, Benoit Pasquier notifications@github.com wrote:

I'm making some progress but it's not easy...

Of course, another solution is to use a Newton-Krylov solver. So when I'm tired of trying the block-Newton approach, I'll try to see if it's easy to just make the current solver use a Krylov iterative solver for each backslash by adding a Krylov flag.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

briochemc commented 5 years ago

Hi François,

I'm not sure (never used any of these iterative solvers before) but there is IncompleteLU.jl for what you are suggesting, and IterativeSolvers.jl with a bunch of potential solutions, too. Worth trying out I guess, because my block-Newton efforts are going nowhere 😅

Cheers, Benoit

briochemc commented 5 years ago

So just to be sure, you recommend GMRES + ILU preconditing?

fprimeau commented 5 years ago

Hi Benoit, GMRES can work, but you should have a look at CTK’s nsoli.m solver because implementation in the context of a Newton solver is a bit more involved. And yes ILU is a viable option for preconditioning.
François.

Sent from my iPhone

On May 27, 2019, at 9:34 PM, Benoit Pasquier notifications@github.com wrote:

So just to be sure, you recommend GMRES + ILU preconditing?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

fprimeau commented 5 years ago

Hi Benoit, Another preconditioning idea is to use a tanime-stepper. Specifically you would use implicit Euler backward scheme for the advection-diffusion but treat all the other tendency terms explicitly using an Euler forward scheme. That way you only have to factor one block, and only once, not at every Newton iteration. The tunable bobs to make the preconditioner effective would be the time step size and the number of time steps to take. You can see that in the limit of many very small time steps this preconditioner is equivalent to the inverse of the Jacobean because it will produce the steady state equilibrium. What we hope for this is that we can take only a few time steps with a reasonably big dt and have the Newton Krylov solver converge fairly quickly.

I have been wanting to try such a preconditioner for a while now but I never find the time to implement and test it.
Best, François.

Sent from my iPhone

On May 27, 2019, at 9:34 PM, Benoit Pasquier notifications@github.com wrote:

So just to be sure, you recommend GMRES + ILU preconditing?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

fprimeau commented 5 years ago

Bobs —> nobs Tanime—> time.

Sent from my iPhone

On May 27, 2019, at 9:34 PM, Benoit Pasquier notifications@github.com wrote:

So just to be sure, you recommend GMRES + ILU preconditing?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

briochemc commented 5 years ago

Hi François,

I have a few questions:

Anyway, I would like to chat about all these points! 🙂 Cheers, Benoit

fprimeau commented 5 years ago

Hi Benoit, You definitely need to use an iterative solver. ILU on it’s own will not work. My point was on;y that gares on it’s own as a replacement for the factorization is complicated to implement. It is better to think holistically of the Newton and Krylov solver.
The nsoli.m implementation is good and you should use that.

We can chat on the phone tomorrow afternoon for me at 4:30pm CA time if that works for you.
Best, Francois. Sent from my iPad

On May 29, 2019, at 6:16 PM, Benoit Pasquier notifications@github.com wrote:

Hi François,

I have a few questions:

After mentioning the use of an incomplete LU, you said

GMRES can work, but

Do you mean we can use the incomplete LU without an iterative solver like GMRES? I.e., directly solve for A * x = b via

Af = ilu(A) x = Af \ b (I do not think that's what you meant but I want to check with you.)

If I understood correctly, you only meant for me to

have a look at CTK’s nsoli.m solver because implementation in the context of a Newton solver is a bit more involved

I just did, but I'm not sure I caught all the "bit more involved" parts. Maybe we can chat about it? (before I try to implement it in AIBECS?)

I am unsure I understood your time-stepper preconditioner concept. Is it similar to the preconditoner used in the Katiwala 2008 paper? I'd like to chat about this.

I have an idea for what I think would be a cheap preconditioner that I want to pass through you. Although I am unsure it would be very efficient. Assuming we geologically restored every tracer, we could use the inverse (the LU factors) of

because its L factor would be the block-diagonal matrix of the L factors of each small matrix, while the big U factor would be the block-diagonal matrix of the U factors. Is that a good or a dumb idea?

Anyway, I would like to chat about all these points! 🙂 Cheers, Benoit

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

briochemc commented 5 years ago

Hi François,

gmres on it’s own as a replacement for the factorization is complicated to implement.

But I would not implement it myself — I would use IterativeSolvers.jl, which has a restarted GMRES implemented already. I think I just need to remove my Shamanskii solver's update_jacobian and δx .= ∇ₓF(x) \ -F(x) (in AIBECS' solver code that was translated from CTK's nsold.m) and use

δx = gmres!(δx, A, F(x), pl = Q)

instead, where gmres! is taken from IterativeSolvers.jl.

FWIW, there is also Isaac.jl, which seems to already be doing the translation of nsoli.m from MATLAB to Julia.

And yes, let's chat tomorrow 🙂

Cheers, Benoit