> On Sep 11, 2020, at 7:45 AM, Blondel, Sophie <[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>
