Sophie,
Great, everything looks good.
So the new version takes about 7 times longer, due to the relatively modest
increase (about 25 percent) in the number of iterations from the poorer
preconditioner convergence and the rest from the much slower matrix-vector
product due to using matrix free instead of matrix based precondtioner. Both of
these are expected.
The matrix is taking about 10% of the memory it used to require, also
expected.
I noticed in the logging the memory for the vectors
Vector 85 85 82303208 0.
Matrix 15 15 8744032 0.
is substantial/huge, with the much smaller matrix memory the vector memory
dominates.
It indicates 85 vectors are used. This is a large number, there are some
needed for the TS (maybe 5?) and some needed for the KSP solve (maybe about 37)
but I am not sure why there are so many. Perhaps this number could be reduced.
Are there are lot of vectors created in the Xolotyl code? I would it could run
with about 45 vectors.
Barry
> On Sep 15, 2020, at 5:12 PM, Blondel, Sophie <[email protected]> wrote:
>
> Hi Barry,
>
> I fixed everything and re-ran the 4 cases in 1D. They took more time than
> before because I used the Kokkos serial backend on the Xolotl side instead of
> the CUDA one previously (long story short, I tried to update CUDA and messed
> up the whole installation). Step 4 looks much better than prevously, I was
> even able to remove MatSetOptions(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE)
> from the code and it ran without throwing errors. The log files are attached.
>
> Cheers,
>
> Sophie
> From: Barry Smith <[email protected]>
> Sent: Friday, September 11, 2020 18:03
> To: Blondel, Sophie <[email protected]>
> Cc: [email protected] <[email protected]>;
> [email protected]
> <[email protected]>
> Subject: Re: [petsc-users] Matrix Free Method questions
>
>
>
>> On Sep 11, 2020, at 7:45 AM, Blondel, Sophie <[email protected]
>> <mailto:[email protected]>> wrote:
>>
>> Thank you Barry,
>>
>> Step 3 worked after I moved MatSetOption at the beginning of
>> computeJacobian(). Attached is the updated log which is pretty similar to
>> what I had before. Step 4 still uses many more iterations.
>>
>> I checked the Jacobian using -ksp_view_pmat ascii (on a simpler case), I can
>> see the difference between step 3 and 4 is that the contribution from the
>> reactions is not included in the step 4 Jacobian (as expected from the fact
>> that I removed their setting from the code).
>>
>> Looking back at one of your previous email, you wrote "This routine should
>> only compute the elements of the Jacobian needed for this reduced matrix
>> Jacobian, so the diagonals and the diffusion/convection terms. ", does it
>> mean that I should still include the contributions from the reactions that
>> affect the pure diagonal terms?
>
> Yes, you need to leave in everything that affects the diagonal otherwise
> the "Jacobi" preconditioner will not reflect the true Jacobi preconditioner
> and likely perform poorly.
>
> Barry
>
>>
>> Cheers,
>>
>> Sophie
>> From: Barry Smith <[email protected] <mailto:[email protected]>>
>> Sent: Thursday, September 10, 2020 17:04
>> To: Blondel, Sophie <[email protected] <mailto:[email protected]>>
>> Cc: [email protected] <mailto:[email protected]>
>> <[email protected] <mailto:[email protected]>>;
>> [email protected]
>> <mailto:[email protected]>
>> <[email protected]
>> <mailto:[email protected]>>
>> Subject: Re: [petsc-users] Matrix Free Method questions
>>
>>
>>
>>> On Sep 10, 2020, at 2:46 PM, Blondel, Sophie <[email protected]
>>> <mailto:[email protected]>> wrote:
>>>
>>> Hi Barry,
>>>
>>> Going through the different changes again to understand what was going
>>> wrong with the last step, I discovered that my changes from 2 to 3
>>> (keeping only the pure diagonal for the reaction Jacobian setup and adding
>>> MatSetOptions(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);) were wrong: the
>>> sparsity of the matrix was correct but then the RHSJacobian method was
>>> wrong. I updated it
>>
>> I'm not sure what you mean here. My hope was that in step 3 you won't
>> need to change RHSJacobian at all (that is just for step 4).
>>
>>> but now when I run step 3 again I get the following error:
>>>
>>> [2]PETSC ERROR: --------------------- Error Message
>>> --------------------------------------------------------------
>>> [2]PETSC ERROR: Argument out of range
>>> [2]PETSC ERROR: Inserting a new nonzero at global row/column (310400,
>>> 316825) into matrix
>>> [2]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html
>>> <https://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting.
>>> [2]PETSC ERROR: Petsc Development GIT revision: v3.13.4-885-gf58a62b032
>>> GIT Date: 2020-09-01 13:07:58 -0500
>>> [2]PETSC ERROR: Unknown Name on a 20200902 named iguazu by bqo Thu Sep 10
>>> 15:38:58 2020
>>> [2]PETSC ERROR: Configure options
>>> PETSC_DIR=/home2/bqo/libraries/petsc-barry PETSC_ARCH=20200902
>>> --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif77 --with-debugging=no
>>> --with-shared-libraries --download-fblaslapack=1
>>> [2]PETSC ERROR: #1 MatSetValues_MPIAIJ() line 606 in
>>> /home2/bqo/libraries/petsc-barry/src/mat/impls/aij/mpi/mpiaij.c
>>> [2]PETSC ERROR: #2 MatSetValues() line 1392 in
>>> /home2/bqo/libraries/petsc-barry/src/mat/interface/matrix.c
>>> [2]PETSC ERROR: #3 MatSetValuesLocal() line 2207 in
>>> /home2/bqo/libraries/petsc-barry/src/mat/interface/matrix.c
>>> [2]PETSC ERROR: #4 MatSetValuesStencil() line 1595 in
>>> /home2/bqo/libraries/petsc-barry/src/mat/interface/matrix.c
>>> PetscSolverExpHandler::computeJacobian: MatSetValuesStencil (reactions)
>>> failed.
>>>
>>> Because the RHSJacobian method is trying to update the elements
>>> corresponding to the reactions. I'm not sure I understood correctly what
>>> step 3 was supposed to be.
>>
>> In step the three the RHSJacobian was suppose to be unchanged, only the
>> option to ignore the "unneeded" Jacobian entries inside MatSetValues (set
>> with MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);) was needed
>> (plus changing the DMDASetBlockFillsXXX argument).
>>
>> The error message Inserting a new nonzero at global row/column (310400,
>> 316825) into matrix indicates that somehow the MatOption
>> MAT_NEW_NONZERO_LOCATION_ERR is in control instead of the option
>> MAT_NEW_NONZERO_LOCATIONS, when it is setting values the Jacobian values.
>>
>> The MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS_ERR,PETSC_TRUE);) is
>> normally called inside the DMCreateMatrix() so I am not sure how they could
>> be getting called in the wrong order but it seems somehow it is
>>
>> When do you call MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);)
>> in the code? You can call it at the beginning of computeJacobian().
>>
>> If this still doesn't work and you get the same error you can run in the
>> debugger on one process and put a breakpoint for MatSetOptions() to found
>> out how the MAT_NEW_NONZERO_LOCATIONS_ERR comes in late to upset the apple
>> cart. You should see MatSetOption() called at least twice and the last one
>> should have the MAT_NEW_NONZERO_LOCATION flag.
>>
>> Barry
>>
>>
>>
>>
>>>
>>> Cheers,
>>>
>>> Sophie
>>>
>>>
>>> From: Barry Smith <[email protected] <mailto:[email protected]>>
>>> Sent: Friday, September 4, 2020 01:06
>>> To: Blondel, Sophie <[email protected] <mailto:[email protected]>>
>>> Cc: [email protected] <mailto:[email protected]>
>>> <[email protected] <mailto:[email protected]>>;
>>> [email protected]
>>> <mailto:[email protected]>
>>> <[email protected]
>>> <mailto:[email protected]>>
>>> Subject: Re: [petsc-users] Matrix Free Method questions
>>>
>>>
>>> Sophie,
>>>
>>> Thanks. I have started looking through the logs
>>>
>>> The change to matrix-free multiple (from 1 to 2) which reduces the
>>> accuracy of the multiply to about half the digits is not surprising.
>>>
>>> * It roughly doubles the time since doing the matrix-free product
>>> requires a function evaluation
>>>
>>> * It increases the iteration count, but not significantly since the
>>> reduced precision of the multiple induces some additional linear iterations
>>>
>>> The change from 2 to 3 (not storing the entire matrix)
>>>
>>> * number of nonzeros goes from 49459966 to 1558766 = 3.15 percent so
>>> it succeds in not storing the unneeded part of the matrix
>>>
>>> * the number of MatMult_MF goes from 2331 to 2418. I don't understand
>>> this, I expected it to be identical because it should be using the same
>>> preconditioner in 3 as in 2 and thus get the same convergence. Could
>>> possibility be due to the variability in convergence due to different runs
>>> with the matrix-free preconditioner preconditioner and not related to not
>>> storing the entire matrix.
>>>
>>> * the KSPSolve() time goes from 3.8774e+0 to 3.7855e+02 a trivial
>>> difference which is what I would expect
>>>
>>> * the SNESSolve time goes from 5.0047e+02 to 4.3275e+02 about a 14
>>> percent drop which is reasonable because 3 doesn't spend as much time
>>> inserting matrix values (it still computes them but doesn't insert the ones
>>> we don't want for the preconditioner).
>>>
>>> The change from 3 to 4
>>>
>>> * something goes seriously wrong here. The total number of linear solve
>>> iterations goes from 2282 to 97403 so something has gone seriously wrong
>>> with the preconditioner, but since the preconditioner operations are the
>>> same it seems something has gone wrong with the new reduced preconditioner.
>>>
>>> I think there is an error in computing the reduced matrix entries, that is
>>> the new compute Jacobian code is not computing the entries it needs to
>>> correctly.
>>>
>>> To debug this you can run case 3 and case 4 for a single time step with
>>> -ksp_view_pmat binary This should create a binary file with the initial
>>> Jacobian matrices in each. You can use Matlab or Python to do the
>>> difference in the matrices and see how possibly the new Jacobian
>>> computation code is not producing the correct values in some locations.
>>>
>>> Good luck,
>>>
>>> Barry
>>>
>>>
>>>
>>>
>>>> On Sep 3, 2020, at 12:26 PM, Blondel, Sophie <[email protected]
>>>> <mailto:[email protected]>> wrote:
>>>>
>>>> Hi Barry,
>>>>
>>>> Attached are the log files for the 1D case, for each of the 4 steps. I
>>>> don't know how I did it yesterday but the differences between steps look
>>>> better today, except for step 4 that takes many more iterations and
>>>> smaller time steps.
>>>>
>>>> Cheers,
>>>>
>>>> Sophie
>>>>
>>>> De : Barry Smith <[email protected] <mailto:[email protected]>>
>>>> Envoyé : mercredi 2 septembre 2020 15:53
>>>> À : Blondel, Sophie <[email protected] <mailto:[email protected]>>
>>>> Cc : [email protected] <mailto:[email protected]>
>>>> <[email protected] <mailto:[email protected]>>;
>>>> [email protected]
>>>> <mailto:[email protected]>
>>>> <[email protected]
>>>> <mailto:[email protected]>>
>>>> Objet : Re: [petsc-users] Matrix Free Method questions
>>>>
>>>>
>>>>
>>>>> On Sep 2, 2020, at 1:44 PM, Blondel, Sophie <[email protected]
>>>>> <mailto:[email protected]>> wrote:
>>>>>
>>>>> Thank you Barry,
>>>>>
>>>>> The code ran with your branch but it's much slower than running with the
>>>>> full Jacobian and Jacobi PC subtype (around 10 times slower). It is using
>>>>> less memory as expected. I tried step 2 as well and it's even slower.
>>>>
>>>> Sophie,
>>>>
>>>> That is puzzling. It should be using the same matrix in the solver so
>>>> should be the same speed and the setup time should be a bit better since
>>>> it does not form the full Jacobian. (We'll get to this later)
>>>>
>>>>> The TS iteration for step 1 are the same as with full Jacobian. Let me
>>>>> know what I can look at to check if I've done something wrong.
>>>>
>>>> We need to see if the KSP iterations are pretty similar for four
>>>> approaches (1) original code with Jacobi PC subtype (2) matrix free with
>>>> Jacobi PC (just add -snes_mf_operator to case 1) (3) the new code with the
>>>> MatSetOption() to not store the entire Jacobian also with the
>>>> -snes_mf_operator and (4) the new code that doesn't compute the unneeded
>>>> part of the Jacobian also with the -snes_mf_operator
>>>>
>>>> You could run each case with same 20 timesteps and -ts_monitor
>>>> -ksp_monitor and -ts_view and send the four output files around.
>>>>
>>>> Once we are sure the four cases are behaving as expected then you can get
>>>> timings for them but let's not do that until we confirm the similar
>>>> results. There could easily be a flaw in my reasoning or the PETSc code
>>>> somewhere that affects the correctness so its best to check that first.
>>>>
>>>>
>>>> Barry
>>>>
>>>>>
>>>>> Cheers,
>>>>>
>>>>> Sophie
>>>>> De : Barry Smith <[email protected] <mailto:[email protected]>>
>>>>> Envoyé : mardi 1 septembre 2020 14:12
>>>>> À : Blondel, Sophie <[email protected] <mailto:[email protected]>>
>>>>> Cc : [email protected] <mailto:[email protected]>
>>>>> <[email protected] <mailto:[email protected]>>;
>>>>> [email protected]
>>>>> <mailto:[email protected]>
>>>>> <[email protected]
>>>>> <mailto:[email protected]>>
>>>>> Objet : Re: [petsc-users] Matrix Free Method questions
>>>>>
>>>>>
>>>>> Sophie,
>>>>>
>>>>> Sorry, looks like an old bug in PETSc that was undetected due to lack
>>>>> of use. The code is trying to use the first of the two matrices to
>>>>> determine the preconditioner which won't work in your case since it is
>>>>> matrix-free. It should be using the second matrix.
>>>>>
>>>>> I hope the branch barry/2020-09-01/fix-fieldsplit-mf resolves this
>>>>> issue for you.
>>>>>
>>>>> Barry
>>>>>
>>>>>
>>>>>> On Sep 1, 2020, at 12:45 PM, Blondel, Sophie <[email protected]
>>>>>> <mailto:[email protected]>> wrote:
>>>>>>
>>>>>> Hi Barry,
>>>>>>
>>>>>> I'm working through step 1) but I think I am doing something wrong. I'm
>>>>>> using DMDASetBlockFillsSparse to set the non-zeros only for the
>>>>>> diffusing clusters (small He clusters here, from size 1 to 7) and all
>>>>>> the diagonal entries. Then I added a few lines in the code:
>>>>>> Mat mat;
>>>>>> DMCreateMatrix(da, &mat);
>>>>>> MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
>>>>>>
>>>>>> When I try to run with the following options: -snes_mf_operator -ts_dt
>>>>>> 1.0e-12 -ts_adapt_time_step_increase_delay 2 -snes_force_iteration
>>>>>> -pc_fieldsplit_detect_coupling -pc_type fieldsplit -fieldsplit_0_pc_type
>>>>>> jacobi -fieldsplit_1_pc_type redundant -ts_max_time 1000.0
>>>>>> -ts_adapt_dt_max 2.0e-3 -ts_adapt_wnormtype INFINITY
>>>>>> -ts_exact_final_time stepover -ts_max_snes_failures -1 -ts_monitor
>>>>>> -ts_max_steps 20
>>>>>>
>>>>>> I get an error:
>>>>>> [0]PETSC ERROR: --------------------- Error Message
>>>>>> --------------------------------------------------------------
>>>>>> [0]PETSC ERROR: No support for this operation for this object type
>>>>>> [0]PETSC ERROR: Matrix type mffd does not have a find off block diagonal
>>>>>> entries defined
>>>>>> [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html
>>>>>> <https://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble
>>>>>> shooting.
>>>>>> [0]PETSC ERROR: Petsc Development GIT revision: v3.13.4-851-gde18fec8da
>>>>>> GIT Date: 2020-08-28 16:47:50 +0000
>>>>>> [0]PETSC ERROR: Unknown Name on a 20200828 named sophie-Precision-5530
>>>>>> by sophie Tue Sep 1 10:58:44 2020
>>>>>> [0]PETSC ERROR: Configure options PETSC_DIR=/home/sophie/Code/petsc
>>>>>> PETSC_ARCH=20200828 --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif77
>>>>>> --with-debugging=no --with-shared-libraries
>>>>>> [0]PETSC ERROR: #1 MatFindOffBlockDiagonalEntries() line 9847 in
>>>>>> /home/sophie/Code/petsc/src/mat/interface/matrix.c
>>>>>> [0]PETSC ERROR: #2 PCFieldSplitSetDefaults() line 504 in
>>>>>> /home/sophie/Code/petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c
>>>>>> [0]PETSC ERROR: #3 PCSetUp_FieldSplit() line 606 in
>>>>>> /home/sophie/Code/petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c
>>>>>> [0]PETSC ERROR: #4 PCSetUp() line 1009 in
>>>>>> /home/sophie/Code/petsc/src/ksp/pc/interface/precon.c
>>>>>> [0]PETSC ERROR: #5 KSPSetUp() line 406 in
>>>>>> /home/sophie/Code/petsc/src/ksp/ksp/interface/itfunc.c
>>>>>> [0]PETSC ERROR: #6 KSPSolve_Private() line 658 in
>>>>>> /home/sophie/Code/petsc/src/ksp/ksp/interface/itfunc.c
>>>>>> [0]PETSC ERROR: #7 KSPSolve() line 889 in
>>>>>> /home/sophie/Code/petsc/src/ksp/ksp/interface/itfunc.c
>>>>>> [0]PETSC ERROR: #8 SNESSolve_NEWTONLS() line 225 in
>>>>>> /home/sophie/Code/petsc/src/snes/impls/ls/ls.c
>>>>>> [0]PETSC ERROR: #9 SNESSolve() line 4524 in
>>>>>> /home/sophie/Code/petsc/src/snes/interface/snes.c
>>>>>> [0]PETSC ERROR: #10 TSStep_ARKIMEX() line 811 in
>>>>>> /home/sophie/Code/petsc/src/ts/impls/arkimex/arkimex.c
>>>>>> [0]PETSC ERROR: #11 TSStep() line 3731 in
>>>>>> /home/sophie/Code/petsc/src/ts/interface/ts.c
>>>>>> [0]PETSC ERROR: #12 TSSolve() line 4128 in
>>>>>> /home/sophie/Code/petsc/src/ts/interface/ts.c
>>>>>> PetscSolver::solve: TSSolve failed.
>>>>>>
>>>>>> Cheers,
>>>>>>
>>>>>> Sophie
>>>>>> De : Barry Smith <[email protected] <mailto:[email protected]>>
>>>>>> Envoyé : lundi 31 août 2020 14:50
>>>>>> À : Blondel, Sophie <[email protected] <mailto:[email protected]>>
>>>>>> Cc : [email protected] <mailto:[email protected]>
>>>>>> <[email protected] <mailto:[email protected]>>;
>>>>>> [email protected]
>>>>>> <mailto:[email protected]>
>>>>>> <[email protected]
>>>>>> <mailto:[email protected]>>
>>>>>> Objet : Re: [petsc-users] Matrix Free Method questions
>>>>>>
>>>>>>
>>>>>> Sophie,
>>>>>>
>>>>>> Thanks.
>>>>>>
>>>>>> The factor of 4 is lot, the 1.5 not so bad.
>>>>>>
>>>>>> You will definitely want to retain the full matrix assembly codes for
>>>>>> speed and to verify a reduced matrix version.
>>>>>>
>>>>>> It is worth trying a "reduced matrix version" with matrix-free
>>>>>> multiply based on these numbers. This reduced matrix Jacobian will only
>>>>>> have the diagonals and all the terms connected to the cluster sizes that
>>>>>> move. In other words you will be building just the part of the Jacobian
>>>>>> needed for the new preconditioner (PC subtype for Jacobi) and doing the
>>>>>> matrix-vector product matrix free. (SOR requires all the Jacobian
>>>>>> entries).
>>>>>>
>>>>>> Fortunately this is hopefully pretty straightforward for this code.
>>>>>> You will not have to change the structure of the main code at all.
>>>>>>
>>>>>> Step 1) create a new "sparse matrix" that will be passed to
>>>>>> DMDASetBlockFillsSparse(). This new "sparse matrix" needs to retain all
>>>>>> the diagonal entries and also all the entries that are associated with
>>>>>> the variables that diffuse. If I remember correctly these are just the
>>>>>> smallest cluster size, plain Helium?
>>>>>>
>>>>>> Call MatSetOptions(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
>>>>>>
>>>>>> Then you would run the code with -snes_mf_operator and the new PC
>>>>>> subtype for Jacobi.
>>>>>>
>>>>>> A test that the new reduced Jacobian is correct will be that you get
>>>>>> almost the same iterations as the runs you just make using the PC
>>>>>> subtype of Jacobi. Hopefully not slower and using a great deal less
>>>>>> memory. The iterations will not be identical because of the matrix-free
>>>>>> multiple.
>>>>>>
>>>>>> Step 2) create a new version of the Jacobian computation routine. This
>>>>>> routine should only compute the elements of the Jacobian needed for this
>>>>>> reduced matrix Jacobian, so the diagonals and the diffusion/convection
>>>>>> terms.
>>>>>>
>>>>>> Again run with with -snes_mf_operator and the new PC subtype for
>>>>>> Jacobi and you should again get the same convergence history.
>>>>>>
>>>>>> I made two steps because it makes it easier to validate and debug to
>>>>>> get the same results as before. The first step cheats in that it still
>>>>>> computes the full Jacobian but ignores the entries that we don't need to
>>>>>> store for the preconditioner. The second step is more efficient because
>>>>>> it only computes the Jacobian entries needed for the preconditioner but
>>>>>> it requires you going through the Jacobian code and making sure only the
>>>>>> needed parts are computed.
>>>>>>
>>>>>>
>>>>>> If you have any questions please let me know.
>>>>>>
>>>>>> Barry
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>> On Aug 31, 2020, at 1:13 PM, Blondel, Sophie <[email protected]
>>>>>>> <mailto:[email protected]>> wrote:
>>>>>>>
>>>>>>> Hi Barry,
>>>>>>>
>>>>>>> I ran the 2 cases to look at the effect of the Jacobi pre-conditionner:
>>>>>>> 1D with 200 grid points and 7759 DOF per grid point (for the PSI
>>>>>>> application), for 20 TS: the factor between SOR and Jacobi is ~4 (976
>>>>>>> MatMult for SOR and 4162 MatMult for Jacobi)
>>>>>>> 2D with 63x63 grid points and 4124 DOF per grid point (for the NE
>>>>>>> application), for 20 TS: the factor is 1.5 (6657 for SOR, 10379 for
>>>>>>> Jacobi)
>>>>>>> Cheers,
>>>>>>>
>>>>>>> Sophie
>>>>>>> De : Barry Smith <[email protected] <mailto:[email protected]>>
>>>>>>> Envoyé : vendredi 28 août 2020 18:31
>>>>>>> À : Blondel, Sophie <[email protected] <mailto:[email protected]>>
>>>>>>> Cc : [email protected] <mailto:[email protected]>
>>>>>>> <[email protected] <mailto:[email protected]>>;
>>>>>>> [email protected]
>>>>>>> <mailto:[email protected]>
>>>>>>> <[email protected]
>>>>>>> <mailto:[email protected]>>
>>>>>>> Objet : Re: [petsc-users] Matrix Free Method questions
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> On Aug 28, 2020, at 4:11 PM, Blondel, Sophie <[email protected]
>>>>>>>> <mailto:[email protected]>> wrote:
>>>>>>>>
>>>>>>>> Thank you Jed and Barry,
>>>>>>>>
>>>>>>>> First, attached are the logs from the benchmark runs I did without
>>>>>>>> (log_std.txt) and with MF method (log_mf.txt). It took me some trouble
>>>>>>>> to get the -log_view to work because I'm using push and pop for the
>>>>>>>> options which means that PETSc is initialized with no argument so the
>>>>>>>> command line argument was not taken into account, but I guess this is
>>>>>>>> for a separate discussion.
>>>>>>>>
>>>>>>>> To answer questions about the current per-conditioners:
>>>>>>>> I used the same pre-conditioner options as listed in my previous email
>>>>>>>> when I added the -snes_mf option; I did try to remove all the PC
>>>>>>>> related options at one point with the MF method but didn't see a
>>>>>>>> change in runtime so I put them back in
>>>>>>>> this benchmark is for a 1D DMDA using 20 grid points; when running in
>>>>>>>> 2D or 3D I switch the PC options to: -pc_type fieldsplit
>>>>>>>> -fieldsplit_0_pc_type sor -fieldsplit_1_pc_type gamg
>>>>>>>> -fieldsplit_1_ksp_type gmres -ksp_type fgmres
>>>>>>>> -fieldsplit_1_pc_gamg_threshold -1
>>>>>>>> I haven't tried a Jacobi PC instead of SOR, I will run a set of more
>>>>>>>> realistic runs (1D and 2D) without MF but with Jacobi and report on it
>>>>>>>> next week. When you say "iterations" do you mean what is given by
>>>>>>>> -ksp_monitor?
>>>>>>>
>>>>>>> Yes, the number of MatMult is a good enough surrogate.
>>>>>>>
>>>>>>> So using matrix-free (which means no preconditioning) has
>>>>>>>
>>>>>>> 35846/160
>>>>>>>
>>>>>>> ans =
>>>>>>>
>>>>>>> 224.0375
>>>>>>>
>>>>>>> or 224 as many iterations. So even for this modest 1d problem
>>>>>>> preconditioning is doing a great deal.
>>>>>>>
>>>>>>> Barry
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>>
>>>>>>>> Cheers,
>>>>>>>>
>>>>>>>> Sophie
>>>>>>>> De : Barry Smith <[email protected] <mailto:[email protected]>>
>>>>>>>> Envoyé : vendredi 28 août 2020 12:12
>>>>>>>> À : Blondel, Sophie <[email protected] <mailto:[email protected]>>
>>>>>>>> Cc : [email protected] <mailto:[email protected]>
>>>>>>>> <[email protected] <mailto:[email protected]>>;
>>>>>>>> [email protected]
>>>>>>>> <mailto:[email protected]>
>>>>>>>> <[email protected]
>>>>>>>> <mailto:[email protected]>>
>>>>>>>> Objet : Re: [petsc-users] Matrix Free Method questions
>>>>>>>>
>>>>>>>> [External Email]
>>>>>>>>
>>>>>>>> Sophie,
>>>>>>>>
>>>>>>>> This is exactly what i would expect. If you run with -ksp_monitor
>>>>>>>> you will see the -snes_mf run takes many more iterations.
>>>>>>>>
>>>>>>>> I am puzzled that the argument -pc_type fieldsplit did not stop the
>>>>>>>> run since this is under normal circumstances not a viable
>>>>>>>> preconditioner with -snes_mf. Did you also remove the -pc_type
>>>>>>>> fieldsplit argument?
>>>>>>>>
>>>>>>>> In order to see how one can avoid forming the entire matrix and use
>>>>>>>> matrix-free to do the matrix-vector but still have an effective
>>>>>>>> preconditioner let's look at what the current preconditioner options
>>>>>>>> do.
>>>>>>>>
>>>>>>>>> -pc_fieldsplit_detect_coupling
>>>>>>>>
>>>>>>>> creates two sub-preconditioners, the first for all the variables and
>>>>>>>> the second for those that are coupled by the matrix to variables in
>>>>>>>> neighboring cells Since only the smallest cluster sizes have
>>>>>>>> diffusion/advection this second set contains only the cluster size one
>>>>>>>> variables.
>>>>>>>>
>>>>>>>>> -fieldsplit_0_pc_type sor
>>>>>>>>
>>>>>>>> Runs SOR on all the variables; you can think of this as running SOR on
>>>>>>>> the reactions, it is a pretty good preconditioner for the reactions
>>>>>>>> since the reactions are local, per cell.
>>>>>>>>
>>>>>>>>> -fieldsplit_1_pc_type redundant
>>>>>>>>
>>>>>>>>
>>>>>>>> This runs the default preconditioner (ILU) on just the variables that
>>>>>>>> diffuse, i.e. the elliptic part. For smallish problems this is fine,
>>>>>>>> for larger problems and 2d and 3d presumably you have also
>>>>>>>> -redundant_pc_type gamg to use algebraic multigrid for the diffusion.
>>>>>>>> This part of the matrix will always need to be formed and used in the
>>>>>>>> preconditioner. It is very important since the diffusion is what
>>>>>>>> brings in most of the ill-conditioning for larger problems into the
>>>>>>>> linear system. Note that it only needs the matrix entries for the
>>>>>>>> cluster size of 1 so it is very small compared to the entire sparse
>>>>>>>> matrix.
>>>>>>>>
>>>>>>>> ----
>>>>>>>> The first preconditioner SOR requires ALL the matrix entries which
>>>>>>>> are almost all (except for the diffusion terms) the coupling between
>>>>>>>> different size clusters within a cell. Especially each cell has its
>>>>>>>> own sparse matrix of the size of total number of clusters, it is
>>>>>>>> sparse but not super sparse.
>>>>>>>>
>>>>>>>> So the to significantly lower memory usage we need to remove the SOR
>>>>>>>> and the storing of all the matrix entries but still have an efficient
>>>>>>>> preconditioner for the "reaction" terms.
>>>>>>>>
>>>>>>>> The simplest thing would be to use Jacobi instead of SOR for the
>>>>>>>> first subpreconditioner since it only requires the diagonal entries in
>>>>>>>> the matrix. But Jacobi is a worse preconditioner than SOR (since it
>>>>>>>> totally ignores the matrix coupling) and sometimes can be much worse.
>>>>>>>>
>>>>>>>> Before anyone writes additional code we need to know if doing
>>>>>>>> something along these lines does not ruin the convergence that.
>>>>>>>>
>>>>>>>> Have you used the same options as before but with
>>>>>>>> -fieldsplit_0_pc_type jacobi ? (Not using any matrix free). We need to
>>>>>>>> get an idea of how many more linear iterations it requires (not time,
>>>>>>>> comparing time won't be helpful for this exercise.) We also need this
>>>>>>>> information for realistic size problems in 2 or 3 dimensions that you
>>>>>>>> really want to run; for small problems this approach will work ok and
>>>>>>>> give misleading information about what happens for large problems.
>>>>>>>>
>>>>>>>> I suspect the iteration counts will shot up. Can you run some cases
>>>>>>>> and see how the iteration counts change?
>>>>>>>>
>>>>>>>> Based on that we can decide if we still retain "good convergence" by
>>>>>>>> changing the SOR to Jacobi and then change the code to make this
>>>>>>>> change efficient (basically by skipping the explicit computation of
>>>>>>>> the reaction Jacobian terms and using matrix-free on the outside of
>>>>>>>> the PCFIELDSPLIT.)
>>>>>>>>
>>>>>>>> Barry
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>> On Aug 28, 2020, at 9:49 AM, Blondel, Sophie via petsc-users
>>>>>>>>> <[email protected] <mailto:[email protected]>> wrote:
>>>>>>>>>
>>>>>>>>> Hi everyone,
>>>>>>>>>
>>>>>>>>> I have been using PETSc for a few years with a fully implicit TS
>>>>>>>>> ARKIMEX method and am now exploring the matrix free method option.
>>>>>>>>> Here is the list of PETSc options I typically use: -ts_dt 1.0e-12
>>>>>>>>> -ts_adapt_time_step_increase_delay 5 -snes_force_iteration
>>>>>>>>> -ts_max_time 1000.0 -ts_adapt_dt_max 2.0e-3 -ts_adapt_wnormtype
>>>>>>>>> INFINITY -ts_exact_final_time stepover -fieldsplit_0_pc_type sor
>>>>>>>>> -ts_max_snes_failures -1 -pc_fieldsplit_detect_coupling -ts_monitor
>>>>>>>>> -pc_type fieldsplit -fieldsplit_1_pc_type redundant -ts_max_steps 100
>>>>>>>>>
>>>>>>>>> I started to compare the performance of the code without changing
>>>>>>>>> anything of the executable and simply adding "-snes_mf", I see a
>>>>>>>>> reduction of memory usage as expected and a benchmark that would
>>>>>>>>> usually take ~5min to run now takes ~50min. Reading the documentation
>>>>>>>>> I saw that there are a few option to play with the matrix free method
>>>>>>>>> like -snes_mf_err, -snes_mf_umin, or switching to -snes_mf_type wp. I
>>>>>>>>> used and modified the values of each of these options separately but
>>>>>>>>> never saw a sizable change in runtime, is it expected?
>>>>>>>>>
>>>>>>>>> And are there other ways to make the matrix free method faster? I saw
>>>>>>>>> in the documentation that you can define your own per-conditioner for
>>>>>>>>> instance. Let me know if you need additional information about the
>>>>>>>>> PETSc setup in the application I use.
>>>>>>>>>
>>>>>>>>> Best,
>>>>>>>>>
>>>>>>>>> Sophie
>>>>>>>>
>>>>>>>> <log_mf.txt><log_std.txt>
>>>>
>>>> <log_1D_1.txt><log_1D_2.txt><log_1D_3.txt><log_1D_4.txt>
>>
>> <log_1D_3_bis.txt>
>
> <log_1D_1.txt><log_1D_2.txt><log_1D_3.txt><log_1D_4.txt>
