January 2024 - petsc-users - mailman3.cels.anl.gov

[Gmsh] Access both default sets and region names
by Noam T. 12 Jan '24

12 Jan '24

Would it be feasible to have an option (e.g. new flag along the lines of -dm_plex_gmsh_...) that allows the user to access both the default sets (Cell / Face / Vertex) together with user-defined gorups (those under $PhysicalNames, available when using -dm_plex_gmsh_use_regions)? That is, with a *.msh file containing $PhysicalNames 2 2 100 "my_surface" 3 200 "my_vol" the return of DMGetLabelName(dm, n, name) would be (order may differ) n = 0, name = "celltype" n = 1, name = "depth" n = 2, name = "Cell Sets" n = 3, name = "my_vol" n = 4, name = "Face Sets" n = 5, name = "my_surface" ... I poked into src/dm/impls/plex/plexgmsh.c and have managed to print all the labels after changing a couple of variable values, so perhaps it is doable. The changes made are not a solution, simply naively set some variables to skip checking for the use_regions flag, so it understandably crashes soon after. Thanks, Noam

2 7

Re: [petsc-users] KSP number of iterations different from size of residual array
by Barry Smith 11 Jan '24

11 Jan '24

Kevin A couple of different things are at play here producing the unexpected results. I have created a merge request https://gitlab.com/petsc/petsc/-/merge_requests/7179 clarifying why the results obtained from the KSPGetResidualHistory() and KSPGetIterationNumber() can be different in the docs. I also fixed a couple of locations of KSPLogResidual() (in gmres and fgmres) that resulted in extra incorrect logging of the history. In summary, with the "standard" textbook Krylov methods, one expects numIts = nEntries - 1, but this need not be the case for advanced Krylov methods (like those with inner iterations or pipelining) or under exceptional circumstances like the use of CG in trust region methods. Barry > On Jan 11, 2024, at 1:20 PM, Kevin G. Wang <kevinw3(a)vt.edu> wrote: > > Hi Barry, > > Thanks for your help!! > > I have checked that in KSPSetResidualHistory, "reset" is set to PETSC_TRUE. I did a few quick tests after reading your message. There seems to be some patterns between "numIts" (given by KSPGetIterationNumber) and "nEntries" (given by KSPGetResidualHistory): > > 1. With gmres or fgmres as the solver: > - If the number of iterations (until error tolerance is met) is small, like 20 - 30, indeed as you said, numIts = nEntries - 1. > - if the number of iterations is large, this is no longer true. I have a case where nEntries = 372, numIts = 360. > 2. With bcgsl, it looks like numIts = 2*(nEntries - 1). > 3. With ibcgs, nEntries = 0, while numIts is nonzero. > > In all these tests, I have set the preconditioner to "none". > > My code (where the KSP functions are called) is here: https://github.com/kevinwgy/m2c/blob/main/LinearSystemSolver.cpp > > I am using PETSc 3.12.4. > > Thanks! > Kevin > > > On Thu, Jan 11, 2024 at 12:26 PM Barry Smith <bsmith(a)petsc.dev <mailto:[email protected]>> wrote: >> >> Trying again. >> >> Normally, numIts would be one less than nEntries since the initial residual is computed (and stored in the history) before any iterations. >> >> Is this what you are seeing or are you seeing other values for the two? >> >> I've started a run of the PETSc test suite that compares the two values for inconsistencies for all tests to see if I can find any problems. >> >> Barry >> >> Also note the importance of the reset value in KSPSetResidualHistory() which means the values will not match when reset is PETSC_FALSE. >> >>> On Jan 10, 2024, at 7:09 PM, Kevin G. Wang <kevinw3(a)vt.edu <mailto:[email protected]>> wrote: >>> >>> Hello everyone! >>> >>> I am writing a code that uses PETSc/KSP to solve linear systems. I just realized that after running "KSPSolve(...)", the number of iterations given by >>> >>> KSPGetIterationNumber(ksp, &numIts) >>> >>> is *different* from the size of the residual history given by >>> >>> KSPGetResidualHistory(ksp, NULL, &nEntries); >>> >>> That is, "numIts" is not equal to "nEntries". Is this expected, or a bug in my code? (I thought they should be the same...) >>> >>> I have tried several pairs of solvers and preconditioners (e.g., fgmres & bjacobi, ibcgs & bjacobi). This issue happens to all of them. >>> >>> Thanks! >>> Kevin >>> >>> -- >>> Kevin G. Wang, Ph.D. >>> Associate Professor >>> Kevin T. Crofton Department of Aerospace and Ocean Engineering >>> Virginia Tech >>> 1600 Innovation Dr., VTSS Rm 224H, Blacksburg, VA 24061 >>> Office: (540) 231-7547 | Mobile: (650) 862-2663 >>> URL: https://www.aoe.vt.edu/people/faculty/wang.html >>> Codes: https://github.com/kevinwgy >> > > > -- > Kevin G. Wang, Ph.D. > Associate Professor > Kevin T. Crofton Department of Aerospace and Ocean Engineering > Virginia Tech > 1600 Innovation Dr., VTSS Rm 224H, Blacksburg, VA 24061 > Office: (540) 231-7547 | Mobile: (650) 862-2663 > URL: https://www.aoe.vt.edu/people/faculty/wang.html > Codes: https://github.com/kevinwgy

1 0

KSP number of iterations different from size of residual array
by Kevin G. Wang 11 Jan '24

11 Jan '24

Hello everyone! I am writing a code that uses PETSc/KSP to solve linear systems. I just realized that after running "KSPSolve(...)", the number of iterations given by KSPGetIterationNumber(ksp, &numIts) is *different* from the size of the residual history given by KSPGetResidualHistory(ksp, NULL, &nEntries); That is, "numIts" is not equal to "nEntries". Is this expected, or a bug in my code? (I thought they should be the same...) I have tried several pairs of solvers and preconditioners (e.g., fgmres & bjacobi, ibcgs & bjacobi). This issue happens to all of them. Thanks! Kevin -- Kevin G. Wang, Ph.D. Associate Professor Kevin T. Crofton Department of Aerospace and Ocean Engineering Virginia Tech 1600 Innovation Dr., VTSS Rm 224H, Blacksburg, VA 24061 Office: (540) 231-7547 | Mobile: (650) 862-2663 URL: https://www.aoe.vt.edu/people/faculty/wang.html Codes: https://github.com/kevinwgy

2 3

Solving reduced system of equations with global preconditioner
by johnny 09 Jan '24

09 Jan '24

Dear Mr. Johnny, I trust this email finds you in good health. I am writing to provide information on the efficient resolution of reduced systems of equations through the incorporation of a global preconditioner. To optimize the convergence of iterative solvers, consider the utilization of a global preconditioner, particularly when dealing with intricate and interconnected equation systems. This strategy can significantly enhance the computational efficiency and accelerate the solution process, especially in scenarios involving large-scale equations. I encourage you to explore various preconditioning techniques and tailor them to the specifics of your problem. Experimentation is essential in identifying the most suitable approach for your application. Should you have any queries or require further clarification on this subject, please do not hesitate to contact me. https://howisdom.com/inside-of-a-blueberry/ Sincerely.

1 0

Solving reduced system of equations with global preconditioner
by johnny 09 Jan '24

09 Jan '24

Hi Mike, I hope this message finds you well. I wanted to reach out regarding the topic of solving reduced systems of equations with a global preconditioner. This approach proves to be highly effective in improving the efficiency of numerical methods for solving complex systems. To delve into this further, consider implementing a global preconditioner to enhance the convergence of your iterative solver when dealing with reduced systems. By doing so, you can significantly accelerate the overall solution process, especially in scenarios with large-scale and interconnected equations. Feel free to explore various preconditioning strategies and adapt them to your specific problem[.](https://artcoolz.com/good-morning-texts.html) Experimentation is key to finding the most suitable approach for your application. If you have any questions or need further assistance on this topic, please don't hesitate to reach out. Best regards, Johnny.

1 0

M2 macs
by Sanjay Govindjee 09 Jan '24

09 Jan '24

I was wondering if anyone has build experience with PETSc + FORTRAN on an M2-based MAC? In particular, I am looking for compiler recommendations. -sanjay

3 2

SNES seems not use my matrix-free operation
by Yi Hu 08 Jan '24

08 Jan '24

Dear PETSc Experts, I am implementing a matrix-free jacobian for my SNES solver in Fortran. (command line option -snes_type newtonls -ksp_type gmres) In the main program, I define my residual and jacobian and matrix-free jacobian like the following, … call DMDASNESSetFunctionLocal(DM_mech, INSERT_VALUES, formResidual, PETSC_NULL_SNES, err_PETSc) call DMSNESSetJacobianLocal(DM_mech, formJacobian, PETSC_NULL_SNES, err_PETSc) … subroutine formJacobian(residual_subdomain,F,Jac_pre,Jac,dummy,err_PETSc) #include <petsc/finclude/petscmat.h> use petscmat implicit None DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & residual_subdomain !< DMDA info (needs to be named "in" for macros like XRANGE to work) real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: & F !< deformation gradient field Mat :: Jac, Jac_pre PetscObject :: dummy PetscErrorCode :: err_PETSc PetscInt :: N_dof ! global number of DoF, maybe only a placeholder N_dof = 9*product(cells(1:2))*cells3 print*, 'in my jac' call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N_dof,N_dof,0,Jac,err_PETSc) CHKERRQ(err_PETSc) call MatShellSetOperation(Jac,MATOP_MULT,GK_op,err_PETSc) CHKERRQ(err_PETSc) print*, 'in my jac' ! for jac preconditioner call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N_dof,N_dof,0,Jac_pre,err_PETSc) CHKERRQ(err_PETSc) call MatShellSetOperation(Jac_pre,MATOP_MULT,GK_op,err_PETSc) CHKERRQ(err_PETSc) print*, 'in my jac' end subroutine formJacobian subroutine GK_op(Jac,dF,output,err_PETSc) real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: & dF real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(out) :: & output real(pREAL), dimension(3,3) :: & deltaF_aim = 0.0_pREAL Mat :: Jac PetscErrorCode :: err_PETSc integer :: i, j, k, e … a lot of calculations … print*, 'in GK op' end subroutine GK_op The first question is that: it seems I still need to explicitly define the interface of MatCreateShell() and MatShellSetOperation() to properly use them, even though I include them via “use petscmat”. It is a little bit strange to me, since some examples do not perform this step. Then the main issue is that I can build my own Jacobian from my call back function formJacobian, and confirm my Jacobian is a shell matrix (by MatView). However, my customized operator GK_op is not called when solving the nonlinear system (not print my “in GK op”). When I try to monitor my SNES, it gives me some conventional output not mentioning my matrix-free operations. So I guess my customized MATOP_MULT may be not associated with Jacobian. Or my configuration is somehow wrong. Could you help me solve this issue? Thanks, Yi ------------------------------------------------- Stay up to date and follow us on LinkedIn, Twitter and YouTube. Max-Planck-Institut für Eisenforschung GmbH Max-Planck-Straße 1 D-40237 Düsseldorf Handelsregister B 2533 Amtsgericht Düsseldorf Geschäftsführung Prof. Dr. Gerhard Dehm Prof. Dr. Jörg Neugebauer Prof. Dr. Dierk Raabe Dr. Kai de Weldige Ust.-Id.-Nr.: DE 11 93 58 514 Steuernummer: 105 5891 1000 Please consider that invitations and e-mails of our institute are only valid if they end with …(a)mpie.de. If you are not sure of the validity please contact rco(a)mpie.de Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails aus unserem Haus nur mit der Endung …(a)mpie.de gültig sind. In Zweifelsfällen wenden Sie sich bitte an rco(a)mpie.de -------------------------------------------------

2 1

Question about MATMPIAIJ and column indexes ordering
by Deuse, Mathieu 08 Jan '24

08 Jan '24

Hello, I have a piece of code which generates a matrix in CSR format, but the without sorting the column indexes in increasing order within each row. This seems not to be 100% compatible with the MATMPIAIJ format: the documentation of MatCreateMPIAIJWithArrays indeed mentions 'row-major ordering'. For example, consider the 2x2 matrix (1 2; 3 4), which in my code could be stored as i=[0, 2, 4], j=[1, 0, 0, 1], v=[2, 1, 3, 4]. I can generate the matrix as follows (on 1 proc): MatCreateMPIAIJWithArrays(PETSC_COMM_SELF, 2, 2, 2, 2, i, j, v, &matrix). This appears to work fine, and I can then use the matrix in a KSP for example. However, if I try to update the entry values (same order and values v=[2, 1, 3, 4]) with MatUpdateMPIAIJWithArray(matrix, v), it seems that PETSc does not memorize the order of the column indexes and the matrix that I get now is (2 1; 3 4). I get the same result with MatUpdateMPIAIJWithArrays(matrix, 2, 2, 2, 2, i, j, v). On the other hand, if the column indexes are sorted within each row (i=[0, 2, 4], j=[0, 1, 0, 1], v=[1, 2, 3, 4]), then it works fine. I have attached a minimal working example (C++). Can I safely rely on MatCreateMPIAIJWithArrays working fine with unsorted column indexes as long as I do not use MatUpdateMPIAIJWithArray(s)? Or should I do the sorting myself before calling MatCreateMPIAIJWithArrays? (or alternatively use another matrix format). Thanks in advance for the help. Kind regards, Mathieu Deuse

2 2

Re: [petsc-users] petsc-users Digest, Vol 181, Issue 9
by Khaled Nabil Shar Abdelaziz 07 Jan '24

07 Jan '24

Issue resolved, stopped using the vector in the struct and used the one in the arguments. Same issue was taking place in the Jacobian function which resulted in error (-8) SNES_DIVERGED_LOCAL_MIN. Adjusting the jacobian function resolved that issue as well. -----Original Message----- From: petsc-users <petsc-users-bounces(a)mcs.anl.gov> On Behalf Of petsc-users-request(a)mcs.anl.gov Sent: Sunday, January 7, 2024 3:52 PM To: petsc-users(a)mcs.anl.gov Subject: petsc-users Digest, Vol 181, Issue 9 ---- External Email: Use caution with attachments, links, or sharing data ---- Send petsc-users mailing list submissions to petsc-users(a)mcs.anl.gov To subscribe or unsubscribe via the World Wide Web, visit https://lists.mcs.anl.gov/mailman/listinfo/petsc-users or, via email, send a message with subject or body 'help' to petsc-users-request(a)mcs.anl.gov You can reach the person managing the list at petsc-users-owner(a)mcs.anl.gov When replying, please edit your Subject line so it is more specific than "Re: Contents of petsc-users digest..." Today's Topics: 1. SNES solve residual vector - Fortran (Khaled Nabil Shar Abdelaziz) ---------------------------------------------------------------------- Message: 1 Date: Sun, 7 Jan 2024 20:51:33 +0000 From: Khaled Nabil Shar Abdelaziz <kabdelaz(a)purdue.edu> To: "petsc-users(a)mcs.anl.gov" <petsc-users(a)mcs.anl.gov> Subject: [petsc-users] SNES solve residual vector - Fortran Message-ID: <PH7PR22MB363758580C91A9958802A9C6BD642(a)PH7PR22MB3637.namprd22.prod.outlook.com> Content-Type: text/plain; charset="us-ascii" Hello everyone, I am running into an issue with SNES solver where the residual in the first iteration gets calculated fine, but on the 2nd iteration, even though it is getting calculated locally, it doesn't seem to be "assigned" to the residual vector correctly. The residual vector here is referred to as petsc_PVector which is defined in a type (struct). It seems like the first iteration, it behaves as expected and that petsc_PVector gets updated, 2nd iteration however seems like there are 2 versions of it; one that is being passed into FormFunctionU, and another that exists in the type (struct). This result in the residual on the 2nd iteration to be exactly = 0 and convergence. Any insights, suggestions, or experiences you can share regarding this issue would be immensely appreciated. I'm particularly interested in understanding why the residual vector behaves differently in the second iteration. Best, Khaled The following is a snapshot of how I initialize and use SNES: subroutine set_function_jacobian_to_solver_u(nlU) type(variable_node_u),target,intent(inout)::nlU PetscErrorCode::ierr type(variable_node_part),pointer::var_nd_part integer::ipart ! external FormFunctionU,FormJacobianU do ipart=1,nlU%variable_node%npart var_nd_part=>nlU%variable_node%parts(ipart) PetscCallA(SNESSetFunction(var_nd_part%petsc_snes,var_nd_part%petsc_PVector,FormFunctionU,nlU,ierr)) PetscCallA(SNESSetJacobian(var_nd_part%petsc_snes,var_nd_part%petsc_KMatrix,var_nd_part%petsc_KMatrix,FormJacobianU,nlU,ierr)) PetscCallA(SNESSetFromOptions(var_nd_part%petsc_snes,ierr)) end do var_nd_part=>null() end subroutine set_function_jacobian_to_solver_u subroutine FormFunctionU(petsc_snes,petsc_vector,petsc_PVector,nlU,ierr) type(tSNES)::petsc_snes type(tVec)::petsc_vector type(tVec)::petsc_PVector type(PetscScalar),pointer::vec_ptr_a(:),vec_ptr_b(:) PetscErrorCode::ierr type(variable_node_u),target::nlU type(PetscScalar),pointer::val_petsc_PVector(:) type(parameter),pointer::para_data type(mesh),pointer::mesh_data type(element),pointer::eles(:) integer::k,i,flag,pid PetscInt::its para_data=>nlU%para_data mesh_data=>nlU%mesh_data eles=>nlU%eles call MPI_Comm_rank(MPI_COMM_WORLD, pid, ierr) call update_data_by_petsc_vector(nlU%variable_node,petsc_vector) call set_petsc_PVector_to_zero(nlU%variable_node) call update_petsc_PVector(nlU%variable_node) flag = 0 do k=1,mesh_data%nel call get_element_value(nlU%data_e, eles(k)%nnd_e, mesh_data%elements(:,k), nlU%nvar_total, nlU%data) call cal_element_ResidualU(nlU%PVector_e, eles(k), nlU%data_e, para_data) call set_petsc_PVector_values_by_element(nlU%variable_node, nlU%PVector_e, eles(k)%nnd_e, mesh_data%elements(:,k)) end do call update_petsc_PVector(nlU%variable_node) para_data=>null() mesh_data=>null() eles=>null() return end subroutine FormFunctionU call SNESSolve(nlU%variable_node%petsc_snes,PETSC_NULL_VEC, nlU%variable_node%petsc_vector,error_petsc)

1 0

SNES solve residual vector - Fortran
by Khaled Nabil Shar Abdelaziz 07 Jan '24

07 Jan '24

Hello everyone, I am running into an issue with SNES solver where the residual in the first iteration gets calculated fine, but on the 2nd iteration, even though it is getting calculated locally, it doesn't seem to be "assigned" to the residual vector correctly. The residual vector here is referred to as petsc_PVector which is defined in a type (struct). It seems like the first iteration, it behaves as expected and that petsc_PVector gets updated, 2nd iteration however seems like there are 2 versions of it; one that is being passed into FormFunctionU, and another that exists in the type (struct). This result in the residual on the 2nd iteration to be exactly = 0 and convergence. Any insights, suggestions, or experiences you can share regarding this issue would be immensely appreciated. I'm particularly interested in understanding why the residual vector behaves differently in the second iteration. Best, Khaled The following is a snapshot of how I initialize and use SNES: subroutine set_function_jacobian_to_solver_u(nlU) type(variable_node_u),target,intent(inout)::nlU PetscErrorCode::ierr type(variable_node_part),pointer::var_nd_part integer::ipart ! external FormFunctionU,FormJacobianU do ipart=1,nlU%variable_node%npart var_nd_part=>nlU%variable_node%parts(ipart) PetscCallA(SNESSetFunction(var_nd_part%petsc_snes,var_nd_part%petsc_PVector,FormFunctionU,nlU,ierr)) PetscCallA(SNESSetJacobian(var_nd_part%petsc_snes,var_nd_part%petsc_KMatrix,var_nd_part%petsc_KMatrix,FormJacobianU,nlU,ierr)) PetscCallA(SNESSetFromOptions(var_nd_part%petsc_snes,ierr)) end do var_nd_part=>null() end subroutine set_function_jacobian_to_solver_u subroutine FormFunctionU(petsc_snes,petsc_vector,petsc_PVector,nlU,ierr) type(tSNES)::petsc_snes type(tVec)::petsc_vector type(tVec)::petsc_PVector type(PetscScalar),pointer::vec_ptr_a(:),vec_ptr_b(:) PetscErrorCode::ierr type(variable_node_u),target::nlU type(PetscScalar),pointer::val_petsc_PVector(:) type(parameter),pointer::para_data type(mesh),pointer::mesh_data type(element),pointer::eles(:) integer::k,i,flag,pid PetscInt::its para_data=>nlU%para_data mesh_data=>nlU%mesh_data eles=>nlU%eles call MPI_Comm_rank(MPI_COMM_WORLD, pid, ierr) call update_data_by_petsc_vector(nlU%variable_node,petsc_vector) call set_petsc_PVector_to_zero(nlU%variable_node) call update_petsc_PVector(nlU%variable_node) flag = 0 do k=1,mesh_data%nel call get_element_value(nlU%data_e, eles(k)%nnd_e, mesh_data%elements(:,k), nlU%nvar_total, nlU%data) call cal_element_ResidualU(nlU%PVector_e, eles(k), nlU%data_e, para_data) call set_petsc_PVector_values_by_element(nlU%variable_node, nlU%PVector_e, eles(k)%nnd_e, mesh_data%elements(:,k)) end do call update_petsc_PVector(nlU%variable_node) para_data=>null() mesh_data=>null() eles=>null() return end subroutine FormFunctionU call SNESSolve(nlU%variable_node%petsc_snes,PETSC_NULL_VEC, nlU%variable_node%petsc_vector,error_petsc)

1 0