[petsc-users] 3D Poisson solver failed in KSPSolve when number of cores is larger than one
Matthew Knepley
knepley at gmail.com
Mon Jul 24 20:49:37 CDT 2023
On Mon, Jul 24, 2023 at 8:16 PM Thuc Bui <bui at calcreek.com> wrote:
> Dear PETSc Users/Developers.
>
>
>
> I have been successfully using PETsc on Windows without MPI for a while
> now. I have now attempted to implement PETSc with MPI on Windows 10. I have
> built a release version of PETSc 3.18.6 with MS MPI 10.1.2, Intel MKL 3.279
> (2020) and Visual Studio 2019 as a static library. I am testing a finite
> difference 3D Poisson solver (cylindrical coordinates) with this MPI PETSc.
>
>
>
> On a Windows 10 laptop with 4 cores (8 logical processors), the solver
> works fine with one core. However, it gives errors with 2 cores. The errors
> are listed below followed by a complete code. The errors are generated in
> KSPSolve. I perhaps miss some settings to cause this error “No method
> muldiagonalblock for Mat of type seqaij”. I don’t know why for a 2-core MPI
> program that a seqaij matrix is used in KSPSolve.
>
>
>
> I would be very grateful if someone can provide me with some direction to
> fix these issues.
>
That is a bug with PCEISENSTAT. It is not frequently used. Since your
matrix does not use inodes, a default implementation of
MatMultDIagonalBlock is not found. We just need to add this for MATSEQAIJ.
It should work with a different PC. We will fix this as soon as we can.
Thanks,
Matt
> Many thanks for your help,
>
> Thuc Bui
>
> Senior R&D Engineer
>
> Calabazas Creek Research, Inc
>
> (650) 948-5361
>
>
>
>
>
> ///////////////errors/////////////////////////
>
> rank=0 iStart=0 iEnd=2375
>
> rank=1 iStart=2375 iEnd=4750
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> --------------------- Error Message
> --------------------------------------------------------------
>
> --------------------- Error Message
> --------------------------------------------------------------
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> No support for this operation for this object type
>
> No support for this operation for this object type
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> No method multdiagonalblock for Mat of type seqaij
>
> No method multdiagonalblock for Mat of type seqaij
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> See https://petsc.org/release/faq/ for trouble shooting.
>
> See https://petsc.org/release/faq/ for trouble shooting.
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> Petsc Release Version 3.18.6, Mar 30, 2023
>
> Petsc Release Version 3.18.6, Mar 30, 2023
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> Poisson3D.exe on a v140-x64SReleaseSP named HERMES by bui Mon Jul 24
> 15:16:55 2023
>
> [1]PETSC ERROR: Configure options --with-cc="win32fe cl" --with-fc=0
> --with-cxx="win32fe cl" --with-debugging=0 --with-openmp
> --with-mpi-include=/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/include:/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/include/x64
> --with-mpi-lib="[/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/Lib/x64/msmpi.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/Lib/x64/msmpifec.lib]"
> --with-blaslapack-lib="[/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_intel_lp64.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_core.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_intel_thread.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/libiomp5md.lib]"
> --with-mpiexec=/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/bin/x64/mpiexec
> --with-shared-libraries=0
>
> Poisson3D.exe on a v140-x64SReleaseSP named HERMES by bui Mon Jul 24
> 15:16:55 2023
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #1 MatMultDiagonalBlock() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538
>
> Configure options --with-cc="win32fe cl" --with-fc=0 --with-cxx="win32fe
> cl" --with-debugging=0 --with-openmp
> --with-mpi-include=/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/include:/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/include/x64
> --with-mpi-lib="[/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/Lib/x64/msmpi.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/Lib/x64/msmpifec.lib]"
> --with-blaslapack-lib="[/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_intel_lp64.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_core.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/mkl_intel_thread.lib,/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/IntelMKL/latest/x64/lib/libiomp5md.lib]"
> --with-mpiexec=/cygdrive/c/Users/bui/Documents/Petsc/externalpackages/MSMPI/latest/bin/x64/mpiexec
> --with-shared-libraries=0
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #2 MatMultDiagonalBlock_MPIAIJ() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\impls\aij\mpi\mpiaij.c:978
>
> #1 MatMultDiagonalBlock() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #3 MatMultDiagonalBlock() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538
>
> #2 MatMultDiagonalBlock_MPIAIJ() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\impls\aij\mpi\mpiaij.c:978
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #4 PCApply_Eisenstat() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\impls\eisens\eisen.c:51
>
> #3 MatMultDiagonalBlock() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\mat\INTERF~1\matrix.c:2538
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #5 PCApply() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\INTERF~1\precon.c:441
>
> #4 PCApply_Eisenstat() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\impls\eisens\eisen.c:51
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #6 KSP_PCApply() at
> C:\Users\bui\DOCUME~1\Petsc\latest\include\petsc/private/kspimpl.h:380
>
> #5 PCApply() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\pc\INTERF~1\precon.c:441
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #7 KSPInitialResidual() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itres.c:64
>
> #6 KSP_PCApply() at
> C:\Users\bui\DOCUME~1\Petsc\latest\include\petsc/private/kspimpl.h:380
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #8 KSPSolve_GMRES() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\impls\gmres\gmres.c:227
>
> #7 KSPInitialResidual() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itres.c:64
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #9 KSPSolve_Private() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:899
>
> #8 KSPSolve_GMRES() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\impls\gmres\gmres.c:227
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #10 KSPSolve() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:1071
>
> #9 KSPSolve_Private() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:899
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> #11 main() at C:\Users\bui\Documents\Nemesis\Poisson3D\Main.cpp:253
>
> #10 KSPSolve() at
> C:\Users\bui\DOCUME~1\Petsc\latest\src\ksp\ksp\INTERF~1\itfunc.c:1071
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> No PETSc Option Table entries
>
> #11 main() at C:\Users\bui\Documents\Nemesis\Poisson3D\Main.cpp:253
>
> [1]PETSC ERROR:
>
> [0]PETSC ERROR:
>
> ----------------End of Error Message -------send entire error message to
> petsc-maint at mcs.anl.gov----------
>
> No PETSc Option Table entries
>
> [0]PETSC ERROR: ----------------End of Error Message -------send entire
> error message to petsc-maint at mcs.anl.gov----------
>
>
>
> //complete codes///////////////////////////////////////////////
>
> #include "petscsys.h"
>
> #include "petscksp.h"
>
> #ifndef PETSC_SUCCESS
>
> #define PETSC_SUCCESS 0
>
> #endif
>
> //rectangular matrix mxn
>
> int** allocateIMatrixR(int m, int n)
>
> {
>
> int** v = (int**)malloc(m * sizeof(int*));
>
> for (int i = 0; i < m; ++i) {
>
> v[i] = (int*)malloc(n * sizeof(int));
>
> }
>
> for (int i = 0; i < m; ++i) {
>
> for (int j = 0; j < n; ++n) {
>
> v[i][j] = -1;
>
> }
>
> }
>
> return v;
>
> }
>
> int mapIndex(int k, int i, int j, int Nz, int Nr)
>
> {
>
> //k for z, i for r, j for theta
>
> //base 0
>
> //3D table: theta r z
>
> // 0 0 0
>
> // 0 0 1
>
> // 0 0 2
>
> // ...
>
> // 0 1 0
>
> // 0 1 1
>
> // 0 1 2
>
> // ...
>
> int n = k + i * Nz + j * Nz * Nr;
>
> return n;
>
> }
>
> PetscErrorCode MPIAssembler(Mat* A, Vec* b, int Kz, int Ir, int Jt)
>
> {
>
> PetscFunctionBeginUser;
>
> int M;
>
> PetscCall(VecGetSize(*b, &M));
>
>
> //contruct an array of Mx3, row index is the equation number, each row is (i,j,k),
>
> //i is r-index, j theta-index, k z-index
>
> const int dim = 3;
>
> int** Nijk = allocateIMatrixR(M, dim);
>
> for (int k = 0; k < Kz; ++k) {
>
> int k1 = k + 1;
>
> for (int i = 0; i < Ir; ++i) {
>
> for (int j = 0; j < Jt; ++j) {
>
> //for each equation row m
>
> int m = mapIndex(k, i, j, Kz, Ir);
>
> Nijk[m][0] = i;
>
> Nijk[m][1] = j;
>
> Nijk[m][2] = k;
>
> }
>
> }
>
> }
>
>
> /////////////////////////////////////////////////////////////////////////////////
>
> //Cylinder dimensions
>
> double R = 11.0; //cm, cylinder radius
>
> double L = 11.0; //cm, cylinder length in z
>
> double Dr = R / (Ir + 1);
>
> double Dz = L / (Kz + 1);
>
> double Dr2 = Dr * Dr;
>
> double Dr_2 = 0.5 * Dr;
>
> double Dr_22 = Dr_2 * Dr_2;
>
> double Dz_2 = 1.0 / (Dz * Dz);
>
> double DT = 2.0 * M_PI / Jt;
>
> double DT2 = DT * DT;
>
> double Vr = 0.0; //Dirichlet on R
>
> double chargeDens = -4.0;
>
> int rank;
>
> PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
>
> //get the range of matrix rows owned by this processor
>
> int startingNumber, endingNumber;
>
>
> PetscCall(MatGetOwnershipRange(*A, &startingNumber, &endingNumber));
>
> fprintf(stderr, "rank=%i iStart=%i iEnd=%i\n"
> , rank, startingNumber, endingNumber);
>
> int Ir_1 = Ir - 1;
>
> int Jt_1 = Jt - 1;
>
> for (int m = startingNumber; m < endingNumber; ++m) {
>
> int i = Nijk[m][0];
>
> int j = Nijk[m][1];
>
> int k = Nijk[m][2];
>
> double ri = Dr * (i + 1);
>
> double ri_2 = ri - Dr_2;
>
> double ri2 = ri + Dr_2;
>
> double sqri = ri * ri;
>
> double Ai = (ri_2 / ri) / Dr2;
>
> double Bi = (ri2 / ri) / Dr2;
>
> double Ci = 1.0 / (sqri * DT2);
>
> double
> Ei = (ri2 + ri_2) / (ri * Dr2) + 2.0 / (sqri * DT2) + 2.0 * Dz_2;
>
> //contribution of charge density to RHS
>
> PetscCall(VecSetValue(*b, m, -chargeDens, ADD_VALUES));
>
> //contribution to RHS by Drichlet on R
>
> if (i == Ir_1) {
>
>
> PetscCall(VecSetValue(*b, i, -Bi * Vr, ADD_VALUES));
>
> }
>
>
>
> if (!i) { //from the axis
>
> //contributing to RHS
>
>
> PetscCall(VecSetValue(*b, m, -chargeDens * Ai * Dr_22, ADD_VALUES));
>
> //contributing to LHS
>
> double xTheta = (Ai / Jt);
>
> for (int jt = 0; jt < Jt; ++jt) {
>
> int nj = mapIndex(k, i, jt, Kz, Ir);
>
>
> PetscCall(MatSetValues(*A, 1, &m, 1, &nj, &xTheta, ADD_VALUES));
>
> }
>
> }
>
> else {
>
> //0 < i < Ir: contribution of Ai to LHS
>
> int i_jk = mapIndex(k, i - 1, j, Kz, Ir);
>
>
> PetscCall(MatSetValues(*A, 1, &m, 1, &i_jk, &Ai, ADD_VALUES));
>
> }
>
> if (i != Ir_1) {
>
> //0 <= i < It-1: contribution of Bi to LHS
>
> int ki1j = mapIndex(k, i + 1, j, Kz, Ir);
>
>
> PetscCall(MatSetValues(*A, 1, &m, 1, &ki1j, &Bi, ADD_VALUES));
>
> }
>
> //apply periodic BC in theta////////////////////////////
>
> int kij_;
>
> if (!j) {
>
> kij_ = mapIndex(k, i, Jt_1, Kz, Ir);
>
> }
>
> else {
>
> kij_ = mapIndex(k, i, j - 1, Kz, Ir);
>
> }
>
> //1st contribution of Ci to LHS
>
>
> PetscCall(MatSetValues(*A, 1, &m, 1, &kij_, &Ci, ADD_VALUES));
>
> int kij1;
>
> if (j == Jt_1) {
>
> kij1 = mapIndex(k, i, 0, Kz, Ir);
>
> }
>
> else {
>
> kij1 = mapIndex(k, i, j + 1, Kz, Ir);
>
> }
>
> //2nd contribution of Ci to LHS
>
>
> PetscCall(MatSetValues(*A, 1, &m, 1, &kij1, &Ci, ADD_VALUES));
>
> ////////////////////////////////////////////////
>
> //apply z Neumann BC ///////////////////////////////////
>
> int k_ij;
>
> if (!k) { //z=0
>
> k_ij = mapIndex(k, i, j, Kz, Ir);
>
> }
>
> else { //0 < z < L
>
> k_ij = mapIndex(k - 1, i, j, Kz, Ir);
>
> }
>
>
> PetscCall(MatSetValues(*A, 1, &m, 1, &k_ij, &Dz_2, ADD_VALUES));
>
> /////////////////////////////////////////////////////////
>
> //contribution of Ei to LHS
>
> int kij = mapIndex(k, i, j, Kz, Ir);
>
>
> PetscCall(MatSetValues(*A, 1, &m, 1, &kij, &Ei, ADD_VALUES));
>
> }
>
> if (Nijk) {
>
> for (int i = 0; i < M; ++i) {
>
> free(Nijk[i]);
>
> }
>
> free(Nijk);
>
> }
>
> PetscCall(VecAssemblyBegin(*b));
>
> PetscCall(VecAssemblyEnd(*b));
>
> PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
>
> PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
>
> PetscFunctionReturn(PETSC_SUCCESS);
>
> }
>
> int main(int argc, char** argv)
>
> {
>
> PetscFunctionBeginUser;
>
> const char help[] = "Poisson3D solver using PETSc.";
>
> PetscCall(PetscInitialize(&argc, &argv, PETSC_NULL, help));
>
> Mat mat_;
>
> Vec B_;
>
> Vec X_;
>
> int Nz = 26; //number of z-segments
>
> int Nr = 20; //number of r-segments
>
> int Nt = 10; //number of theta-segments
>
> PetscCall(PetscOptionsGetInt(NULL, NULL, "-z", &Nz, NULL));
>
> PetscCall(PetscOptionsGetInt(NULL, NULL, "-r", &Nr, NULL));
>
> PetscCall(PetscOptionsGetInt(NULL, NULL, "-t", &Nt, NULL));
>
> //Neumann at z=0 and z=L
>
> //Dirichlet at R, the cylinder radius
>
> //On axis, compute field
>
> //number of unknowns of cylinder interior nodes minus axis
>
> int Kz = Nz - 1;
>
> int Ir = Nr - 1;
>
> int M = Kz * Ir * Nt;
>
> int bw = 15; //for 7-point stencil
>
> PetscCall(MatCreate(PETSC_COMM_WORLD, &mat_));
>
> PetscCall(MatSetSizes(mat_, PETSC_DECIDE, PETSC_DECIDE, M, M));
>
> PetscCall(MatSetFromOptions(mat_));
>
> PetscCall(MatMPIAIJSetPreallocation(mat_, bw, 0, bw, 0));
>
> PetscCall(MatSeqAIJSetPreallocation(mat_, bw, 0));
>
> //source
>
> PetscCall(VecCreate(PETSC_COMM_WORLD, &B_));
>
> PetscCall(VecSetSizes(B_, PETSC_DECIDE, M));
>
> PetscCall(VecSetFromOptions(B_));
>
> //sol
>
> PetscCall(VecDuplicate(B_, &X_));
>
> MPIAssembler(&mat_, &B_, Kz, Ir, Nt);
>
> PCType pcType = PCEISENSTAT;
>
> KSPType solverType = KSPGMRES;
>
> int maxIters = 100;
>
> double relTol = 1.0e-6;
>
> KSP ksp;
>
> PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
>
> PetscCall(KSPSetOperators(ksp, mat_, mat_));
>
> PetscCall(KSPSetType(ksp, solverType));
>
> //preconditioner
>
> PC pc;
>
> PetscCall(KSPGetPC(ksp, &pc));
>
> PetscCall(PCSetType(pc, pcType));
>
>
> PetscCall(KSPSetTolerances(ksp, relTol, PETSC_DEFAULT, PETSC_DEFAULT, maxIters));
>
> PetscCall(KSPSetFromOptions(ksp));
>
> //solving
>
> PetscCall(KSPSolve(ksp, B_, X_));
>
> double xMin, xMax;
>
> int mnIndx, mxIndx;
>
> PetscCall(VecMin(X_, &mnIndx, &xMin));
>
> PetscCall(VecMax(X_, &mxIndx, &xMax));
>
> fprintf(stderr, "mnIndx=%i mxIndx=%i xMin=%g xMax=%g\n"
> , mnIndx, mxIndx, xMin, xMax);
>
>
>
> PetscCall(PetscFinalize());
>
> return 0;
>
> }
>
>
>
--
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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230724/4b3c04ff/attachment-0001.html>
More information about the petsc-users
mailing list