September 2021 - petsc-users - mailman3.cels.anl.gov

How to combine different element types into a single DMPlex?
by Eric Chamberland 01 Aug '24

01 Aug '24

Hi, while playing with DMPlexBuildFromCellListParallel, I noticed we have to specify "numCorners" which is a fixed value, then gives a fixed number of nodes for a series of elements. How can I then add, for example, triangles and quadrangles into a DMPlex? Thanks, Eric -- Eric Chamberland, ing., M. Ing Professionnel de recherche GIREF/Université Laval (418) 656-2131 poste 41 22 42

4 13

Nodes coordinates in distributed dmplex
by Alejandro D Otero 16 Jul '22

16 Jul '22

Hi, sorry if this is an obvious question, but I cannot figure out how to recover finite element nodes coordinates once I have distributed a mesh stored as a dmplex. I am using petsc4py as interface to petsc rutines. I first created a dmplex using: dm.createFromCellList() In a sequential run I got the coordinates with: Coords = dm.getCoordinates() which gave a sequential vector with the coordinates of the mesh nodes. When I distribute the mesh with: dm.distribute() each mpi process has it own dm but the indexing of the vector resulting from getCoordinates() or getCoordinatesLocal() seems not consistent with the local numbering of the cells and nodes. Which is the correct way of doing this in PETSc philosophy? Thanks in advance, Alejandro

4 9

DMView and DMLoad
by Berend van Wachem 09 Mar '22

09 Mar '22

Dear PETSc-team, We are trying to save and load distributed DMPlex and its associated physical fields (created with DMCreateGlobalVector) (Uvelocity, VVelocity, ...) in HDF5_XDMF format. To achieve this, we do the following: 1) save in the same xdmf.h5 file: DMView( DM , H5_XDMF_Viewer ); VecView( UVelocity, H5_XDMF_Viewer ); 2) load the dm: DMPlexCreateFromfile(PETSC_COMM_WORLD, Filename, PETSC_TRUE, DM); 3) load the physical field: VecLoad( UVelocity, H5_XDMF_Viewer ); There are no errors in the execution, but the loaded DM is distributed differently to the original one, which results in the incorrect placement of the values of the physical fields (UVelocity etc.) in the domain. This approach is used to restart the simulation with the last saved DM. Is there something we are missing, or there exists alternative routes to this goal? Can we somehow get the IS of the redistribution, so we can re-distribute the vector data as well? Many thanks, best regards, Berend.

5 13

Question about PCFieldSplit
by Jorti, Zakariae 31 Jan '22

31 Jan '22

Hello, I have a DMStag grid with one dof on each edge and face center. I want to use a PCFieldSplit preconditioner on a Jacobian matrix that I assume is already split but I am not sure how to determine the fields. In the DMStag examples (ex2.c and ex3.c), the function PCFieldSplitSetDetectSaddlePoint is used to determine those fields based on zero diagonal entries. In my case, I have a Jacobian matrix that does not have zero diagonal entries. Can I use that PCFieldSplitSetDetectSaddlePoint in this case? If not, how should I do? Should I do like this example (https://www.mcs.anl.gov/petsc/petsc-master/src/ksp/ksp/tutorials/ex43.c.html) const PetscInt Bfields[1] = {0},Efields[1] = {1}; KSPGetPC(ksp,&pc); PCFieldSplitSetBlockSize(pc,2); PCFieldSplitSetFields(pc,"B",1,Bfields,Bfields); PCFieldSplitSetFields(pc,"E",1,Efields,Efields); where my B unknowns are defined on face centers and E unknowns are defined on edge centers? One last thing, I do not know which field comes first. Is it the one defined for face dofs or edge dofs. Thank you. Best regards, Zakariae

4 14

PC shell destroy
by Alfredo J Duarte Gomez 19 Nov '21

19 Nov '21

Good afternoon PETSC team, I am currently developing an application for PETSC in which I use my own preconditioner with a PCSHELL. I have successfully set all the functions and the performance of the preconditioner is good. I am using this PCSHELL within a TS object, and it is imperative that the objects in the PCSHELL context are freed every time since the memory requirements of those are large. I have set up the Preconditioner before the TS starts with the following block of code: ---------------------------------------------------------------------------------------------------- ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr); ierr = ShellPCCreate(&shell);CHKERRQ(ierr); ierr = PCShellSetApply(pc,MatrixFreePreconditioner);CHKERRQ(ierr); ierr = PCShellSetContext(pc,shell);CHKERRQ(ierr); ierr = PCShellSetDestroy(pc,ShellPCDestroy);CHKERRQ(ierr); ierr = PCShellSetName(pc,"MyPreconditioner");CHKERRQ(ierr); ierr = ShellPCSetUp(pc,da,0.0,dt,u,user);CHKERRQ(ierr); ierr = TSSetPreStep(ts,PreStep);CHKERRQ(ierr); ------------------------------------------------------------------------------------------------ The shell context is then updated by using the following code within the TSPreStep function: --------------------------------------------------------------------------- ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); // Get necessary objects from TS context TSGetTime(ts,&time); TSGetApplicationContext(ts,&user); TSGetSolution(ts,&X); TSGetTimeStep(ts,&dt); TSGetStepNumber(ts, &stepi); TSGetDM(ts,&da); tdt = time+dt; // Update preconditioner context with current values ierr = ShellPCSetUp(pc,da,tdt,dt,X,user);CHKERRQ(ierr); --------------------------------------------------------------------------- I have set up the necessary code in the function ShellPCDestroy to free the objects within this context, however I am unsure when/if this function is called automatically. Do I have to free the context myself after every step? How would I call the function myself? I am running out of memory after a few steps, and I think this shell context is the culprit. In addition to that, is it possible to get what is called the "ashift" dF/dU + a*dF/dU_t in this function from the TS object? https://petsc.org/release/docs/manualpages/TS/TSSetIJacobian.html I need it as an input for my preconditioner (currrently hardcoded for TSBEULER where ashift is always 1/dt). Thank you, -Alfredo -- Alfredo Duarte Graduate Research Assistant The University of Texas at Austin

3 3

Is it possible to keep track of original elements # after a call to DMPlexDistribute ?
by Eric Chamberland 11 Nov '21

11 Nov '21

Hi, I want to use DMPlexDistribute from PETSc for computing overlapping and play with the different partitioners supported. However, after calling DMPlexDistribute, I noticed the elements are renumbered and then the original number is lost. What would be the best way to keep track of the element renumbering? a) Adding an optional parameter to let the user retrieve a vector or "IS" giving the old number? b) Adding a DMLabel (seems a wrong good solution) c) Other idea? Of course, I don't want to loose performances with the need of this "mapping"... Thanks, Eric -- Eric Chamberland, ing., M. Ing Professionnel de recherche GIREF/Université Laval (418) 656-2131 poste 41 22 42

2 31

Disconnected domains and Poisson equation
by Marco Cisternino 06 Oct '21

06 Oct '21

Good morning, I want to solve the Poisson equation on a 3D domain with 2 non-connected sub-domains. I am using FGMRES+GAMG and I have no problem if the two sub-domains see a Dirichlet boundary condition each. On the same domain I would like to solve the Poisson equation imposing periodic boundary condition in one direction and homogenous Neumann boundary conditions in the other two directions. The two sub-domains are symmetric with respect to the separation between them and the operator discretization and the right hand side are symmetric as well. It would be nice to have the same solution in both the sub-domains. Setting the null space to the constant, the solver converges to a solution having the same gradients in both sub-domains but different values. Am I doing some wrong with the null space? I'm not setting a block matrix (one block for each sub-domain), should I? I tested the null space against the matrix using MatNullSpaceTest and the answer is true. Can I do something more to have a symmetric solution as outcome of the solver? Thank you in advance for any comments and hints. Best regards, Marco Cisternino

4 13

Re: [petsc-users] (percent time in this phase)
by Karthikeyan Chockalingam - STFC UKRI 05 Oct '21

05 Oct '21

When comparing the MatSolve data for GPU MatSolve 341 1.0 1.3009e+02 1.6 2.96e+11 1.0 0.0e+00 0.0e+00 0.0e+00 36 34 0 0 0 36 34 0 0 0 4536 4538 1 5.34e+02 0 0.00e+00 100 MatCUSPARSSolAnl 2 1.0 3.2338e+01 1.5 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 9 0 0 0 0 9 0 0 0 0 0 0 0 0.00e+00 0 0.00e+00 0 and CPU MatSolve 352 1.0 1.3553e+02 1.0 1.02e+11 1.0 0.0e+00 0.0e+00 0.0e+00 35 34 0 0 0 35 34 0 0 0 4489 the time spent is almost the same for this preconditioner. Look like MatCUSPARSSolAnl is called only twice (since I am running on two cores) mpirun -n 2 ./ex45 -log_summary -da_grid_x 511 -da_grid_y 511 -da_grid_z 511 -dm_mat_type mpiaijcusparse -dm_vec_type mpicuda -ksp_type cg -pc_type bjacobi -ksp_monitor So would it be fair to assume MatCUSPARSSolAnl is not accounted for in MatSolve and it is an exclusive event? KSPSolve (53%) + PCSetup (16%) + DMCreateMat (23%) + MatCUSPARSSolAnl (9%) ~ 100 % Best, Karthik. From: Matthew Knepley <knepley(a)gmail.com> Date: Wednesday, 29 September 2021 at 16:29 To: "Chockalingam, Karthikeyan (STFC,DL,HC)" <karthikeyan.chockalingam(a)stfc.ac.uk> Cc: Barry Smith <bsmith(a)petsc.dev>, "petsc-users(a)mcs.anl.gov" <petsc-users(a)mcs.anl.gov> Subject: Re: [petsc-users] %T (percent time in this phase) On Wed, Sep 29, 2021 at 10:18 AM Karthikeyan Chockalingam - STFC UKRI <karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]>> wrote: Thank you! Just to summarize KSPSolve (53%) + PCSetup (16%) + DMCreateMat (23%) + MatCUSPARSSolAnl (9%) ~ 100 % You didn’t happen to mention how MatCUSPARSSolAnl is accounted for? Am I right in accounting for it as above? I am not sure.I thought it might be the GPU part of MatSolve(). I will have to look in the code. I am not as familiar with the GPU part. MatCUSPARSSolAnl 2 1.0 3.2338e+01 1.5 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 9 0 0 0 0 9 0 0 0 0 0 0 0 0.00e+00 0 0.00e+00 0 Finally, I believe the vector events, VecNorn, VecTDot, VecAXPY, and VecAYPX are mutually exclusive? Yes. Thanks, Matt Best, Karthik. From: Matthew Knepley <knepley(a)gmail.com<mailto:[email protected]>> Date: Wednesday, 29 September 2021 at 11:58 To: "Chockalingam, Karthikeyan (STFC,DL,HC)" <karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]>> Cc: Barry Smith <bsmith(a)petsc.dev<mailto:[email protected]>>, "petsc-users(a)mcs.anl.gov<mailto:[email protected]>" <petsc-users(a)mcs.anl.gov<mailto:[email protected]>> Subject: Re: [petsc-users] %T (percent time in this phase) On Wed, Sep 29, 2021 at 6:24 AM Karthikeyan Chockalingam - STFC UKRI <karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]>> wrote: Thank you Mathew. Now, it is all making sense to me. From data file ksp_ex45_N511_gpu_2.txt KSPSolve (53%) + KSPSetup (0%) = PCSetup (16%) + PCApply (37%). However, you said “So an iteration would mostly consist of MatMult + PCApply, with some vector work” 1) You do one solve, but 2 KSPSetUp()s. You must be running on more than one process and using Block-Jacobi . Half the time is spent in the solve (53%) KSPSetUp 2 1.0 5.3149e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 2.4e+01 0 0 0 0 1 0 0 0 0 1 0 0 0 0.00e+00 0 0.00e+00 0 KSPSolve 1 1.0 1.5837e+02 1.1 8.63e+11 1.0 6.8e+02 2.1e+06 4.4e+03 53100100100 95 53100100100 96 10881 11730 1022 6.40e+03 1021 8.17e-03 100 2) The preconditioner look like BJacobi-ILU. The setup time is 16%, which is all setup of the individual blocks, and this is all used by the numerical ILU factorization. PCSetUp 2 1.0 4.9623e+01 1.3 1.45e+09 1.0 0.0e+00 0.0e+00 0.0e+00 16 0 0 0 0 16 0 0 0 0 58 0 2 6.93e+03 0 0.00e+00 0 PCSetUpOnBlocks 1 1.0 4.9274e+01 1.3 1.45e+09 1.0 0.0e+00 0.0e+00 0.0e+00 15 0 0 0 0 15 0 0 0 0 59 0 2 6.93e+03 0 0.00e+00 0 MatLUFactorNum 1 1.0 4.6126e+01 1.3 1.45e+09 1.0 0.0e+00 0.0e+00 0.0e+00 14 0 0 0 0 14 0 0 0 0 63 0 2 6.93e+03 0 0.00e+00 0 MatILUFactorSym 1 1.0 2.5110e+00 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 1 0 0 0 0 1 0 0 0 0 0 0 0 0.00e+00 0 0.00e+00 0 3) The preconditioner application takes 37% of the time, which is all solving the factors and recorded in MatSolve(). Matrix multiplication takes 4%. PCApply 341 1.0 1.3068e+02 1.6 2.96e+11 1.0 0.0e+00 0.0e+00 0.0e+00 37 34 0 0 0 37 34 0 0 0 4516 4523 1 5.34e+02 0 0.00e+00 100 MatSolve 341 1.0 1.3009e+02 1.6 2.96e+11 1.0 0.0e+00 0.0e+00 0.0e+00 36 34 0 0 0 36 34 0 0 0 4536 4538 1 5.34e+02 0 0.00e+00 100 MatMult 341 1.0 1.0774e+01 1.1 2.96e+11 1.0 6.9e+02 2.1e+06 2.0e+00 4 34100100 0 4 34100100 0 54801 66441 2 5.86e+03 0 0.00e+00 100 4) The significant vector time is all in norms (11%) since they are really slow on the GPU. VecNorm 342 1.0 6.2261e+01129.9 4.57e+10 1.0 0.0e+00 0.0e+00 6.8e+02 11 5 0 0 15 11 5 0 0 15 1466 196884 0 0.00e+00 342 2.74e-03 100 VecTDot 680 1.0 1.7107e+00 1.3 9.09e+10 1.0 0.0e+00 0.0e+00 1.4e+03 1 10 0 0 29 1 10 0 0 29 106079 133922 0 0.00e+00 680 5.44e-03 100 VecAXPY 681 1.0 3.2036e+00 1.7 9.10e+10 1.0 0.0e+00 0.0e+00 0.0e+00 1 11 0 0 0 1 11 0 0 0 56728 58367 682 5.34e+02 0 0.00e+00 100 VecAYPX 339 1.0 2.6502e+00 1.8 4.53e+10 1.0 0.0e+00 0.0e+00 0.0e+00 1 5 0 0 0 1 5 0 0 0 34136 34153 339 2.71e-03 0 0.00e+00 100 So the solve time is: 53% ~ 37% + 4% + 11% and the setup time is about 16%. I was wrong about the SetUp time being included, as it is outside the event: https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/interface/itfunc.c#L… It looks like the remainder of the time (23%) is spent preallocating the matrix. Thanks, Matt The MalMult event is 4 %. How does this event figure into the above equation; if preconditioning (MatMult + PCApply) is included in KSPSolve? Best, Karthik. From: Matthew Knepley <knepley(a)gmail.com<mailto:[email protected]>> Date: Wednesday, 29 September 2021 at 10:58 To: "Chockalingam, Karthikeyan (STFC,DL,HC)" <karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]>> Cc: Barry Smith <bsmith(a)petsc.dev<mailto:[email protected]>>, "petsc-users(a)mcs.anl.gov<mailto:[email protected]>" <petsc-users(a)mcs.anl.gov<mailto:[email protected]>> Subject: Re: [petsc-users] %T (percent time in this phase) On Wed, Sep 29, 2021 at 5:52 AM Karthikeyan Chockalingam - STFC UKRI <karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]>> wrote: That was helpful. I would like to provide some additional details of my run on cpus and gpus. Please find the following attachments: 1. graph.pdf a plot showing overall time and various petsc events. 2. ksp_ex45_N511_cpu_6.txt data file of the log_summary 3. ksp_ex45_N511_gpu_2.txt data file of the log_summary I used the following petsc options for cpu mpirun -n 6 ./ex45 -log_summary -da_grid_x 511 -da_grid_y 511 -da_grid_z 511 -dm_mat_type mpiaij -dm_vec_type mpi -ksp_type cg -pc_type bjacobi -ksp_monitor and for gpus mpirun -n 1 ./ex45 -log_summary -da_grid_x 511 -da_grid_y 511 -da_grid_z 511 -dm_mat_type mpiaijcusparse -dm_vec_type mpicuda -ksp_type cg -pc_type bjacobi -ksp_monitor to run the following problem https://petsc.org/release/src/ksp/ksp/tutorials/ex45.c.html From the above code, I see is there no individual function called KSPSetUp(), so I gather KSPSetDM, KSPSetComputeInitialGuess, KSPSetComputeRHS, kSPSetComputeOperators all are timed together as KSPSetUp. For this example, is KSPSetUp time and KSPSolve time mutually exclusive? No, KSPSetUp() will be contained in KSPSolve() if it is called automatically. In your response you said that “PCSetUp() time may be in KSPSetUp() or it maybe in PCApply() it depends on how much of the preconditioner construction can take place early, so depends exactly on the preconditioner used.” I don’t see a explicit call to PCSetUp() or PCApply() in ex45; so for this particular preconditioner (bjacobi) how can I tell how they are timed? They are all inside KSPSolve(). If you have a preconditioned linear solve, the oreconditioning happens during the iteration. So an iteration would mostly consist of MatMult + PCApply, with some vector work. I am hoping to time KSP solving and preconditioning mutually exclusively. I am not sure that concept makes sense here. See above. Thanks, Matt Kind regards, Karthik. From: Barry Smith <bsmith(a)petsc.dev<mailto:[email protected]>> Date: Tuesday, 28 September 2021 at 19:19 To: "Chockalingam, Karthikeyan (STFC,DL,HC)" <karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]>> Cc: "petsc-users(a)mcs.anl.gov<mailto:[email protected]>" <petsc-users(a)mcs.anl.gov<mailto:[email protected]>> Subject: Re: [petsc-users] %T (percent time in this phase) On Sep 28, 2021, at 12:11 PM, Karthikeyan Chockalingam - STFC UKRI <karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]>> wrote: Thanks for Barry for your response. I was just benchmarking the problem with various preconditioner on cpu and gpu. I understand, it is not possible to get mutually exclusive timing. However, can you tell if KSPSolve time includes both PCSetup and PCApply? And if KSPSolve and KSPSetup are mutually exclusive? Likewise for PCSetUp and PCApply. If you do not call KSPSetUp() separately from KSPSolve() then its time is included with KSPSolve(). PCSetUp() time may be in KSPSetUp() or it maybe in PCApply() it depends on how much of the preconditioner construction can take place early, so depends exactly on the preconditioner used. So yes the answer is not totally satisfying. The one thing I would recommend is to not call KSPSetUp() directly and then KSPSolve() will always include the total time of the solve plus all setup time. PCApply will contain all the time to apply the preconditioner but may also include some setup time. Barry Best, Karthik. From: Barry Smith <bsmith(a)petsc.dev<mailto:[email protected]>> Date: Tuesday, 28 September 2021 at 16:56 To: "Chockalingam, Karthikeyan (STFC,DL,HC)" <karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]>> Cc: "petsc-users(a)mcs.anl.gov<mailto:[email protected]>" <petsc-users(a)mcs.anl.gov<mailto:[email protected]>> Subject: Re: [petsc-users] %T (percent time in this phase) On Sep 28, 2021, at 10:55 AM, Karthikeyan Chockalingam - STFC UKRI <karthikeyan.chockalingam(a)stfc.ac.uk<mailto:[email protected]>> wrote: Hello, I ran ex45 in the KPS tutorial, which is a 3D finite-difference Poisson problem. I noticed from the output from using the flag -log_summary that for various events their respective %T (percent time in this phase) do not add up to 100 but rather exceeds 100. So, I gather there is some overlap among these events. I am primarily looking at the events KSPSetUp, KSPSolve, PCSetUp and PCSolve. Is it possible to get a mutually exclusive %T or Time for these individual events? I have attached the log_summary output file from my run for your reference. For nested solvers it is tricky to get the times to be mutually exclusive because some parts of the building of the preconditioner is for some preconditioners delayed until the solve has started. It looks like you are using the default preconditioner options which for this example are taking more or less no time since so many iterations are needed. It is best to use -pc_type mg to use geometric multigrid on this problem. Barry Thanks! Karthik. This email and any attachments are intended solely for the use of the named recipients. If you are not the intended recipient you must not use, disclose, copy or distribute this email or any of its attachments and should notify the sender immediately and delete this email from your system. UK Research and Innovation (UKRI) has taken every reasonable precaution to minimise risk of this email or any attachments containing viruses or malware but the recipient should carry out its own virus and malware checks before opening the attachments. UKRI does not accept any liability for any losses or damages which the recipient may sustain due to presence of any viruses. <ksp_ex45_N511_cpu_6.txt> -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/> -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/> -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>

3 18

SLEPc: smallest eigenvalues
by Varun Hiremath 05 Oct '21

05 Oct '21

Hi All, I am trying to compute the smallest eigenvalues of a generalized system A*x= lambda*B*x. I don't explicitly know the matrix A (so I am using a shell matrix with a custom matmult function) however, the matrix B is explicitly known so I compute inv(B)*A within the shell matrix and solve inv(B)*A*x = lambda*x. To compute the smallest eigenvalues it is recommended to solve the inverted system, but since matrix A is not explicitly known I can't invert the system. Moreover, the size of the system can be really big, and with the default Krylov solver, it is extremely slow. So is there a better way for me to compute the smallest eigenvalues of this system? Thanks, Varun

2 20

Convergence rate for spatially varying Helmholtz system
by Ramakrishnan Thirumalaisamy 01 Oct '21

01 Oct '21

Hi all, I am trying to solve the Helmholtz equation for temperature T: (C I + Div D grad) T = f in IBAMR, in which C is the spatially varying diagonal entries, and D is the spatially varying diffusion coefficient. I use a matrix-free solver with matrix-based PETSc preconditioner. For the matrix-free solver, I use gmres solver and for the matrix based preconditioner, I use Richardson ksp + Jacobi as a preconditioner. As the simulation progresses, the iterations start to increase. To understand the cause, I set D to be zero, which results in a diagonal system: C T = f. This should result in convergence within a single iteration, but I get convergence in 3 iterations. Residual norms for temperature_ solve. 0 KSP preconditioned resid norm 4.590811647875e-02 true resid norm 2.406067589273e+09 ||r(i)||/||b|| 4.455533946945e-05 1 KSP preconditioned resid norm 2.347767895880e-06 true resid norm 1.210763896685e+05 ||r(i)||/||b|| 2.242081505717e-09 2 KSP preconditioned resid norm 1.245406571896e-10 true resid norm 6.328828824310e+00 ||r(i)||/||b|| 1.171966730978e-13 Linear temperature_ solve converged due to CONVERGED_RTOL iterations 2 To verify that I am indeed solving a diagonal system I printed the PETSc matrix from the preconditioner and viewed it in Matlab. It indeed shows it to be a diagonal system. Attached is the plot of the spy command on the printed matrix. The matrix in binary form is also attached. My understanding is that because the C coefficient is varying in 4 orders of magnitude, i.e., Max(C)/Min(C) ~ 10^4, the matrix is poorly scaled. When I rescale my matrix by 1/C then the system converges in 1 iteration as expected. Is my understanding correct, and that scaling 1/C should be done even for a diagonal system? When D is non-zero, then scaling by 1/C seems to be very inconvenient as D is stored as side-centered data for the matrix free solver. In the case that I do not scale my equations by 1/C, is there some solver setting that improves the convergence rate? (With D as non-zero, I have also tried gmres as the ksp solver in the matrix-based preconditioner to get better performance, but it didn't matter much.) Thanks, Ramakrishnan Thirumalaisamy San Diego State University.

4 10