Re: [petsc-users] Cannot eagerly initialize cuda, as doing so results in cuda error 35 (cudaErrorInsufficientDriver) : CUDA driver version is insufficient for CUDA runtime version

2022-01-19 Thread Jed Brown
Fande Kong  writes:

> On Wed, Jan 19, 2022 at 11:39 AM Jacob Faibussowitsch 
> wrote:
>
>> Are you running on login nodes or compute nodes (I can’t seem to tell from
>> the configure.log)?
>>
>
> I was compiling codes on login nodes, and running codes on compute nodes.
> Login nodes do not have GPUs, but compute nodes do have GPUs.
>
> Just to be clear, the same thing (code, machine) with PETSc-3.16.1 worked
> perfectly. I have this trouble with PETSc-main.

I assume you can

export PETSC_OPTIONS='-device_enable lazy'

and it'll work.

I think this should be the default. The main complaint is that timing the first 
GPU-using event isn't accurate if it includes initialization, but I think this 
is mostly hypothetical because you can't trust any timing that doesn't preload 
in some form and the first GPU-using event will almost always be something 
uninteresting so I think it will rarely lead to confusion. Meanwhile, eager 
initialization is viscerally disruptive for lots of people.


[petsc-users] Does mpiaijkok intend to support 64-bit integers?

2022-01-19 Thread Fande Kong
Hi All,

It seems that mpiaijkok does not support 64-bit integers at this time. Do
we have any motivation for this? Or Is it just a bug?

Thanks,

Fande

petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx(306): error: a
value of type "MatColumnIndexType *" cannot be assigned to an entity of
type "int *"


petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx(308): error: a
value of type "PetscInt *" cannot be assigned to an entity of type "int *"


petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx(310): error: a
value of type "PetscInt *" cannot be assigned to an entity of type "int *"


petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx(316): error: a
value of type "MatColumnIndexType *" cannot be assigned to an entity of
type "int *"


petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx(329): error: a
value of type "PetscInt *" cannot be assigned to an entity of type "int *"


petsc/src/mat/impls/aij/mpi/kokkos/mpiaijkok.kokkos.cxx(331): error: a
value of type "PetscInt *" cannot be assigned to an entity of type "int *"


6 errors detected in the compilation of
"/tmp/tmpxft_00017e46_-6_mpiaijkok.kokkos.cpp1.ii".

gmake[3]: *** [arch-linux-c-opt/obj/mat/impls/aij/mpi/kokkos/mpiaijkok.o]
Error 1


Re: [petsc-users] PETSc MUMPS interface

2022-01-19 Thread Zhang, Hong via petsc-users
Varun,
This feature is merged to petsc main
https://gitlab.com/petsc/petsc/-/merge_requests/4727
Hong

From: petsc-users  on behalf of Zhang, Hong 
via petsc-users 
Sent: Wednesday, January 19, 2022 9:37 AM
To: Varun Hiremath 
Cc: Peder Jørgensgaard Olesen via petsc-users 
Subject: Re: [petsc-users] PETSc MUMPS interface

Varun,
Good to know it works. FactorSymbolic function is still being called twice, but 
the 2nd call is a no-op, thus it still appears in '-log_view'. I made changes 
in the low level of mumps routine, not within PCSetUp() because I feel your use 
case is limited to mumps, not other matrix package solvers.
Hong

From: Varun Hiremath 
Sent: Wednesday, January 19, 2022 2:44 AM
To: Zhang, Hong 
Cc: Peder Jørgensgaard Olesen via petsc-users 
Subject: Re: [petsc-users] PETSc MUMPS interface

Hi Hong,

Thanks, I tested your branch and I think it is working fine. I don't see any 
increase in runtime, however with -log_view I see that the MatLUFactorSymbolic 
function is still being called twice, so is this expected? Is the second call a 
no-op?

$ ./ex52.o -use_mumps_lu -print_mumps_memory -log_view | grep MatLUFactorSym
MatLUFactorSym 2 1.0 4.4411e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 
0.0e+00  2  0  0  0  0   2  0  0  0  0 0

Thanks,
Varun

On Mon, Jan 17, 2022 at 7:49 PM Zhang, Hong 
mailto:hzh...@mcs.anl.gov>> wrote:
Varun,
I created a branch hzhang/feature-mumps-mem-estimate,
see https://gitlab.com/petsc/petsc/-/merge_requests/4727

You may give it a try and let me know if this is what you want.
src/ksp/ksp/tutorials/ex52.c is an example.

Hong

From: Varun Hiremath mailto:varunhirem...@gmail.com>>
Sent: Monday, January 17, 2022 12:41 PM
To: Zhang, Hong mailto:hzh...@mcs.anl.gov>>
Cc: Jose E. Roman mailto:jro...@dsic.upv.es>>; Peder 
Jørgensgaard Olesen via petsc-users 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] PETSc MUMPS interface

Hi Hong,

Thanks for looking into this. Here is the workflow that I might use:

MatLUFactorSymbolic(F,A,perm,iperm,);

// get memory estimates from MUMPS e.g. INFO(3), INFOG(16), INFOG(17)
// find available memory on the system e.g. RAM size
if (estimated_memory > available_memory)
{
   // inform and stop; or
   // switch MUMPS to out-of-core factorization
   ICNTL(22) = 1;
}
else
{
   // set appropriate settings for in-core factorization
}

// Now we call the solve and inside if MatLUFactorSymbolic is already called 
then it should be skipped
EPSSolve(eps);

Thanks,
Varun

On Mon, Jan 17, 2022 at 9:18 AM Zhang, Hong 
mailto:hzh...@mcs.anl.gov>> wrote:
Varun,
I am trying to find a way to enable you to switch options after 
MatLUFactorSymbolic(). A hack is modifying the flag  'mumps->matstruc' inside  
MatLUFactorSymbolic_AIJMUMPS() and MatFactorNumeric_MUMPS().

My understanding of what you want is:
  // collect mumps memory info
  ...
  MatLUFactorSymbolic(F,A,perm,iperm,);
  printMumpsMemoryInfo(F);
  //-
if (memory is available) {
EPSSolve(eps); --> skip calling of MatLUFactorSymbolic()
} else {
   //out-of-core (OOC) option in MUMPS
}

Am I correct? I'll let you know once I work out a solution.
Hong


From: Varun Hiremath mailto:varunhirem...@gmail.com>>
Sent: Sunday, January 16, 2022 10:10 PM
To: Zhang, Hong mailto:hzh...@mcs.anl.gov>>
Cc: Jose E. Roman mailto:jro...@dsic.upv.es>>; Peder 
Jørgensgaard Olesen via petsc-users 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] PETSc MUMPS interface

Hi Jose, Hong,

Thanks for the explanation. I have verified using -log_view that 
MatLUFactorSymbolic is indeed getting called twice.

Hong, we use MUMPS solver for other things, and we typically run the symbolic 
analysis first and get memory estimates to ensure that we have enough memory 
available to run the case. If the available memory is not enough, we can stop 
or switch to the out-of-core (OOC) option in MUMPS. We wanted to do the same 
when using MUMPS via SLEPc/PETSc. Please let me know if there are other ways of 
getting these memory stats and switching options during runtime with PETSc.
Appreciate your help!

Thanks,
Varun

On Sun, Jan 16, 2022 at 4:01 PM Zhang, Hong 
mailto:hzh...@mcs.anl.gov>> wrote:
Varun,
I believe Jose is correct. You may verify it by running your code with option 
'-log_view', then check the number of calls to MatLUFactorSym.

I guess I can add a flag in PCSetUp() to check if user has already called 
MatLUFactorSymbolic() and wants to skip it. Normally, users simply allocate 
sufficient memory in the symbolic factorization. Why do you want to check it?
Hong


From: Jose E. Roman mailto:jro...@dsic.upv.es>>
Sent: Sunday, January 16, 2022 5:11 AM
To: Varun Hiremath mailto:varunhirem...@gmail.com>>
Cc: Zhang, Hong mailto:hzh...@mcs.anl.gov>>; Peder 
Jørgensgaard Olesen via petsc-users 
mailto:petsc-users@mcs.anl.gov>>

Re: [petsc-users] Cannot eagerly initialize cuda, as doing so results in cuda error 35 (cudaErrorInsufficientDriver) : CUDA driver version is insufficient for CUDA runtime version

2022-01-19 Thread Fande Kong
On Wed, Jan 19, 2022 at 11:39 AM Jacob Faibussowitsch 
wrote:

> Are you running on login nodes or compute nodes (I can’t seem to tell from
> the configure.log)?
>

I was compiling codes on login nodes, and running codes on compute nodes.
Login nodes do not have GPUs, but compute nodes do have GPUs.

Just to be clear, the same thing (code, machine) with PETSc-3.16.1 worked
perfectly. I have this trouble with PETSc-main.

I might do "git bisect" when I have time

Thanks,

Fande


If running from login nodes, do they support running with GPU’s? Some
> clusters will install stub versions of cuda runtime on login nodes (such
> that configuration can find them), but that won’t actually work in
> practice.
>
> If this is the case then CUDA will fail to initialize with this exact
> error. IIRC It wasn’t until CUDA 11.1 that they created a specific error
> code (cudaErrorStubLibrary) for it.
>
> Best regards,
>
> Jacob Faibussowitsch
> (Jacob Fai - booss - oh - vitch)
>
> On Jan 19, 2022, at 12:07, Fande Kong  wrote:
>
> Thanks, Jacob, and Junchao
>
> The log was attached.  I am using Sawtooth at INL
> https://hpc.inl.gov/SitePages/Home.aspx
>
>
> Thanks,
>
> Fande
>
> On Wed, Jan 19, 2022 at 10:32 AM Jacob Faibussowitsch 
> wrote:
>
>> Hi Fande,
>>
>> What machine are you running this on? Please attach configure.log so I
>> can troubleshoot this.
>>
>> Best regards,
>>
>> Jacob Faibussowitsch
>> (Jacob Fai - booss - oh - vitch)
>>
>> On Jan 19, 2022, at 10:04, Fande Kong  wrote:
>>
>> Hi All,
>>
>> Upgraded PETSc from 3.16.1 to the current main branch. I suddenly got the
>> following error message:
>>
>> 2d_diffusion]$ ../../../moose_test-dbg -i 2d_diffusion_test.i
>> -use_gpu_aware_mpi 0 -gpu_mat_type aijcusparse -gpu_vec_type cuda
>>  -log_view
>> [0]PETSC ERROR: - Error Message
>> --
>> [0]PETSC ERROR: Missing or incorrect user input
>> [0]PETSC ERROR: Cannot eagerly initialize cuda, as doing so results in
>> cuda error 35 (cudaErrorInsufficientDriver) : CUDA driver version is
>> insufficient for CUDA runtime version
>> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
>> [0]PETSC ERROR: Petsc Development GIT revision: v3.16.3-618-gad32f7e  GIT
>> Date: 2022-01-18 16:04:31 +
>> [0]PETSC ERROR: ../../../moose_test-dbg on a arch-linux-c-opt named
>> r8i3n0 by kongf Wed Jan 19 08:30:13 2022
>> [0]PETSC ERROR: Configure options --with-debugging=no
>> --with-shared-libraries=1 --download-fblaslapack=1 --download-metis=1
>> --download-ptscotch=1 --download-parmetis=1 --download-mumps=1
>> --download-strumpack=1 --download-scalapack=1 --download-slepc=1
>> --with-mpi=1 --with-cxx-dialect=C++14 --with-fortran-bindings=0
>> --with-sowing=0 --with-64-bit-indices --with-make-np=24 --with-cuda
>> --with-cudac=nvcc --with-cuda-arch=70 --download-kokkos=1
>> [0]PETSC ERROR: #1 initialize() at
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/impls/cupm/cupmdevice.cxx:298
>> [0]PETSC ERROR: #2 PetscDeviceInitializeTypeFromOptions_Private() at
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/interface/device.cxx:299
>> [0]PETSC ERROR: #3 PetscDeviceInitializeFromOptions_Internal() at
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/interface/device.cxx:425
>> [0]PETSC ERROR: #4 PetscInitialize_Common() at
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/pinit.c:963
>> [0]PETSC ERROR: #5 PetscInitialize() at
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/pinit.c:1238
>> [0]PETSC ERROR: #6 SlepcInitialize() at
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/arch-linux-c-opt/externalpackages/git.slepc/src/sys/slepcinit.c:275
>> [0]PETSC ERROR: #7 LibMeshInit() at ../src/base/libmesh.C:522
>> [r8i3n0:mpi_rank_0][MPIDI_CH3_Abort] application called
>> MPI_Abort(MPI_COMM_WORLD, 95) - process 0: No such file or directory (2)
>>
>> Thanks,
>>
>> Fande
>>
>>
>> 
>
>
>


Re: [petsc-users] Cannot eagerly initialize cuda, as doing so results in cuda error 35 (cudaErrorInsufficientDriver) : CUDA driver version is insufficient for CUDA runtime version

2022-01-19 Thread Jacob Faibussowitsch
Are you running on login nodes or compute nodes (I can’t seem to tell from the 
configure.log)? If running from login nodes, do they support running with 
GPU’s? Some clusters will install stub versions of cuda runtime on login nodes 
(such that configuration can find them), but that won’t actually work in 
practice. 

If this is the case then CUDA will fail to initialize with this exact error. 
IIRC It wasn’t until CUDA 11.1 that they created a specific error code 
(cudaErrorStubLibrary) for it.

Best regards,

Jacob Faibussowitsch
(Jacob Fai - booss - oh - vitch)

> On Jan 19, 2022, at 12:07, Fande Kong  wrote:
> 
> Thanks, Jacob, and Junchao
> 
> The log was attached.  I am using Sawtooth at INL 
> https://hpc.inl.gov/SitePages/Home.aspx 
> 
> 
> 
> Thanks,
> 
> Fande
> 
> On Wed, Jan 19, 2022 at 10:32 AM Jacob Faibussowitsch  > wrote:
> Hi Fande,
> 
> What machine are you running this on? Please attach configure.log so I can 
> troubleshoot this.
> 
> Best regards,
> 
> Jacob Faibussowitsch
> (Jacob Fai - booss - oh - vitch)
> 
>> On Jan 19, 2022, at 10:04, Fande Kong > > wrote:
>> 
>> Hi All,
>> 
>> Upgraded PETSc from 3.16.1 to the current main branch. I suddenly got the 
>> following error message:
>> 
>> 2d_diffusion]$ ../../../moose_test-dbg -i 2d_diffusion_test.i 
>> -use_gpu_aware_mpi 0 -gpu_mat_type aijcusparse -gpu_vec_type cuda  -log_view 
>> [0]PETSC ERROR: - Error Message 
>> --
>> [0]PETSC ERROR: Missing or incorrect user input 
>> [0]PETSC ERROR: Cannot eagerly initialize cuda, as doing so results in cuda 
>> error 35 (cudaErrorInsufficientDriver) : CUDA driver version is insufficient 
>> for CUDA runtime version
>> [0]PETSC ERROR: See https://petsc.org/release/faq/ 
>>  for trouble shooting.
>> [0]PETSC ERROR: Petsc Development GIT revision: v3.16.3-618-gad32f7e  GIT 
>> Date: 2022-01-18 16:04:31 +
>> [0]PETSC ERROR: ../../../moose_test-dbg on a arch-linux-c-opt named r8i3n0 
>> by kongf Wed Jan 19 08:30:13 2022
>> [0]PETSC ERROR: Configure options --with-debugging=no 
>> --with-shared-libraries=1 --download-fblaslapack=1 --download-metis=1 
>> --download-ptscotch=1 --download-parmetis=1 --download-mumps=1 
>> --download-strumpack=1 --download-scalapack=1 --download-slepc=1 
>> --with-mpi=1 --with-cxx-dialect=C++14 --with-fortran-bindings=0 
>> --with-sowing=0 --with-64-bit-indices --with-make-np=24 --with-cuda 
>> --with-cudac=nvcc --with-cuda-arch=70 --download-kokkos=1
>> [0]PETSC ERROR: #1 initialize() at 
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/impls/cupm/cupmdevice.cxx:298
>> [0]PETSC ERROR: #2 PetscDeviceInitializeTypeFromOptions_Private() at 
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/interface/device.cxx:299
>> [0]PETSC ERROR: #3 PetscDeviceInitializeFromOptions_Internal() at 
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/interface/device.cxx:425
>> [0]PETSC ERROR: #4 PetscInitialize_Common() at 
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/pinit.c:963
>> [0]PETSC ERROR: #5 PetscInitialize() at 
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/pinit.c:1238
>> [0]PETSC ERROR: #6 SlepcInitialize() at 
>> /home/kongf/workhome/sawtooth/moosegpu/petsc/arch-linux-c-opt/externalpackages/git.slepc/src/sys/slepcinit.c:275
>> [0]PETSC ERROR: #7 LibMeshInit() at ../src/base/libmesh.C:522
>> [r8i3n0:mpi_rank_0][MPIDI_CH3_Abort] application called 
>> MPI_Abort(MPI_COMM_WORLD, 95) - process 0: No such file or directory (2)
>> 
>> Thanks,
>> 
>> Fande
> 
> 



Re: [petsc-users] Nullspaces

2022-01-19 Thread Mark Adams
On Wed, Jan 19, 2022 at 12:54 PM Marco Cisternino <
marco.cistern...@optimad.it> wrote:

> Thank you, Matthew.
>
> I’m going to pay attention to our non-dimensionalization, avoiding
> division by cell volume helps a lot.
>
>
>
> Sorry, Mark, I cannot get your point: which 1D problem are you referring
> to? The case I’m talking about is based on a 3D octree mesh.
>
>
>
Woops, getting my threads mixed up.


Re: [petsc-users] Nullspaces

2022-01-19 Thread Marco Cisternino
Thank you, Matthew.
I’m going to pay attention to our non-dimensionalization, avoiding division by 
cell volume helps a lot.

Sorry, Mark, I cannot get your point: which 1D problem are you referring to? 
The case I’m talking about is based on a 3D octree mesh.

Thank you all for your support.

Bests,
Marco Cisternino

From: Mark Adams 
Sent: mercoledì 19 gennaio 2022 15:16
To: Matthew Knepley 
Cc: Marco Cisternino ; petsc-users 

Subject: Re: [petsc-users] Nullspaces

ILU is LU for this 1D problem and its singular.
ILU might have some logic to deal with a singular system. Not sure. LU should 
fail.
You might try -ksp_type preonly and -pc_type lu
And -pc_type jacobi should not have any numerical problems. Try that.

On Wed, Jan 19, 2022 at 7:19 AM Matthew Knepley 
mailto:knep...@gmail.com>> wrote:
On Wed, Jan 19, 2022 at 4:52 AM Marco Cisternino 
mailto:marco.cistern...@optimad.it>> wrote:
Thank you Matthew.
But I cannot get the point. I got the point about the test but to try to 
explain my doubt I’m going to prepare another toy code.

By words…
I usually have a finite volume discretization of the Laplace operator with 
homogeneous Neumann BC on an octree mesh and it reads
Aij * xj = bi,
being the discretization of
Int|Vi(nabla^2 pi dV) = Int|Vi(nabla dot ui)
(I omit constants), where Int|Vi(…dV) is a volume integral on the I cell, pi is 
cell pressure, ui is the cell velocity.
The computational domain contains 2 separated sub-domains.
Let’s consider 4 cells into the whole domain and 2 cells for each sub-domain.
I would expect a null space made of 2 vectors and from your patch they should 
look like [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)], i.e. norm2 = 
1 for both.
With MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr, ) 
I’m saying that [1/sqrt(4) 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)] is the null space, 
which is not but it is in the null space.
But this is not the case I sent you, not exactly.

The case I sent is 1/Vi * Aij * xj = 1/Vi bi, where Vi is the volume of the 
cell i.
Let’s admit that yj is in the null space of of Aij, it should be in the null 
space of 1/Vi * Aij, therefore Aij*yj = 0 and 1/Vi * Aij*yj = 0 too.
But in the framework of the test, this is true with infinite precision.
What happens if norm2(Aij*yj) = 10^-15 and Vi = 10^-5? Norm2(1/Vi * Aij * yj) = 
10^-10!!! Is yi still in the null space numerically?
Let’s say yi is the constant vector over the whole domain, i.e. [1/sqrt(4) 
1/sqrt(4) 1/sqrt(4) 1/sqrt(4)]. Should this be in the null space of 1/Vi * Aij, 
shouldn’t it?
An analogous argument should be for the compatibility condition that concerns 
bi.

My current problem is that properly creating the null space for Aij, i.e. 
[1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)], allows me to solve and 
find xi, but multiplying by 1/Vi, I cannot get any solution using both 
FGMRES+ILU and FGMRE+GAMG.

The tiny problem will load Aij, Vi and bi and show the problem by testing the 
proper null space and trying to solve. I will include the patch to my PETSc 
version. I hope to come back to you very soon.

This sounds like a dimensionalization problem to me. It is best to choose 
length (and other) units that make the matrix entries about 1. It seems
like you are far from this, and it is causing a loss of accuracy (your null 
vector has a residual of about 1e-7).

  Thanks,

  Matt

Thank you very much for your support!


Marco Cisternino

From: Matthew Knepley mailto:knep...@gmail.com>>
Sent: martedì 18 gennaio 2022 21:25
To: Marco Cisternino 
mailto:marco.cistern...@optimad.it>>
Cc: petsc-users mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Nullspaces

I made a fix for this:

  https://gitlab.com/petsc/petsc/-/merge_requests/4729

  Thanks,

 Matt

On Tue, Jan 18, 2022 at 3:20 PM Matthew Knepley 
mailto:knep...@gmail.com>> wrote:
On Thu, Dec 16, 2021 at 11:09 AM Marco Cisternino 
mailto:marco.cistern...@optimad.it>> wrote:
Hello Matthew,
as promised I prepared a minimal (112960 rows. I’m not able to produce anything 
smaller than this and triggering the issue) example of the behavior I was 
talking about some days ago.
What I did is to produce matrix, right hand side and initial solution of the 
linear system.

As I told you before, this linear system is the discretization of the pressure 
equation of a predictor-corrector method for NS equations in the framework of 
finite volume method.
This case has homogeneous Neumann boundary conditions. Computational domain has 
two independent and separated sub-domains.
I discretize the weak formulation and I divide every row of the linear system 
by the volume of the relative cell.
The underlying mesh is not uniform, therefore cells have different volumes.
The issue I’m going to explain does not show up if the mesh is uniform, same 
volume for all the cells.

I usually build the null space sub-domain by sub-domain with
MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, nConstants, constants, 
);

Re: [petsc-users] Cannot eagerly initialize cuda, as doing so results in cuda error 35 (cudaErrorInsufficientDriver) : CUDA driver version is insufficient for CUDA runtime version

2022-01-19 Thread Jacob Faibussowitsch
Hi Fande,

What machine are you running this on? Please attach configure.log so I can 
troubleshoot this.

Best regards,

Jacob Faibussowitsch
(Jacob Fai - booss - oh - vitch)

> On Jan 19, 2022, at 10:04, Fande Kong  wrote:
> 
> Hi All,
> 
> Upgraded PETSc from 3.16.1 to the current main branch. I suddenly got the 
> following error message:
> 
> 2d_diffusion]$ ../../../moose_test-dbg -i 2d_diffusion_test.i 
> -use_gpu_aware_mpi 0 -gpu_mat_type aijcusparse -gpu_vec_type cuda  -log_view 
> [0]PETSC ERROR: - Error Message 
> --
> [0]PETSC ERROR: Missing or incorrect user input 
> [0]PETSC ERROR: Cannot eagerly initialize cuda, as doing so results in cuda 
> error 35 (cudaErrorInsufficientDriver) : CUDA driver version is insufficient 
> for CUDA runtime version
> [0]PETSC ERROR: See https://petsc.org/release/faq/ 
>  for trouble shooting.
> [0]PETSC ERROR: Petsc Development GIT revision: v3.16.3-618-gad32f7e  GIT 
> Date: 2022-01-18 16:04:31 +
> [0]PETSC ERROR: ../../../moose_test-dbg on a arch-linux-c-opt named r8i3n0 by 
> kongf Wed Jan 19 08:30:13 2022
> [0]PETSC ERROR: Configure options --with-debugging=no 
> --with-shared-libraries=1 --download-fblaslapack=1 --download-metis=1 
> --download-ptscotch=1 --download-parmetis=1 --download-mumps=1 
> --download-strumpack=1 --download-scalapack=1 --download-slepc=1 --with-mpi=1 
> --with-cxx-dialect=C++14 --with-fortran-bindings=0 --with-sowing=0 
> --with-64-bit-indices --with-make-np=24 --with-cuda --with-cudac=nvcc 
> --with-cuda-arch=70 --download-kokkos=1
> [0]PETSC ERROR: #1 initialize() at 
> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/impls/cupm/cupmdevice.cxx:298
> [0]PETSC ERROR: #2 PetscDeviceInitializeTypeFromOptions_Private() at 
> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/interface/device.cxx:299
> [0]PETSC ERROR: #3 PetscDeviceInitializeFromOptions_Internal() at 
> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/interface/device.cxx:425
> [0]PETSC ERROR: #4 PetscInitialize_Common() at 
> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/pinit.c:963
> [0]PETSC ERROR: #5 PetscInitialize() at 
> /home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/pinit.c:1238
> [0]PETSC ERROR: #6 SlepcInitialize() at 
> /home/kongf/workhome/sawtooth/moosegpu/petsc/arch-linux-c-opt/externalpackages/git.slepc/src/sys/slepcinit.c:275
> [0]PETSC ERROR: #7 LibMeshInit() at ../src/base/libmesh.C:522
> [r8i3n0:mpi_rank_0][MPIDI_CH3_Abort] application called 
> MPI_Abort(MPI_COMM_WORLD, 95) - process 0: No such file or directory (2)
> 
> Thanks,
> 
> Fande



[petsc-users] Cannot eagerly initialize cuda, as doing so results in cuda error 35 (cudaErrorInsufficientDriver) : CUDA driver version is insufficient for CUDA runtime version

2022-01-19 Thread Fande Kong
Hi All,

Upgraded PETSc from 3.16.1 to the current main branch. I suddenly got the
following error message:

2d_diffusion]$ ../../../moose_test-dbg -i 2d_diffusion_test.i
-use_gpu_aware_mpi 0 -gpu_mat_type aijcusparse -gpu_vec_type cuda
 -log_view
[0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: Missing or incorrect user input
[0]PETSC ERROR: Cannot eagerly initialize cuda, as doing so results in cuda
error 35 (cudaErrorInsufficientDriver) : CUDA driver version is
insufficient for CUDA runtime version
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.16.3-618-gad32f7e  GIT
Date: 2022-01-18 16:04:31 +
[0]PETSC ERROR: ../../../moose_test-dbg on a arch-linux-c-opt named r8i3n0
by kongf Wed Jan 19 08:30:13 2022
[0]PETSC ERROR: Configure options --with-debugging=no
--with-shared-libraries=1 --download-fblaslapack=1 --download-metis=1
--download-ptscotch=1 --download-parmetis=1 --download-mumps=1
--download-strumpack=1 --download-scalapack=1 --download-slepc=1
--with-mpi=1 --with-cxx-dialect=C++14 --with-fortran-bindings=0
--with-sowing=0 --with-64-bit-indices --with-make-np=24 --with-cuda
--with-cudac=nvcc --with-cuda-arch=70 --download-kokkos=1
[0]PETSC ERROR: #1 initialize() at
/home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/impls/cupm/cupmdevice.cxx:298
[0]PETSC ERROR: #2 PetscDeviceInitializeTypeFromOptions_Private() at
/home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/interface/device.cxx:299
[0]PETSC ERROR: #3 PetscDeviceInitializeFromOptions_Internal() at
/home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/device/interface/device.cxx:425
[0]PETSC ERROR: #4 PetscInitialize_Common() at
/home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/pinit.c:963
[0]PETSC ERROR: #5 PetscInitialize() at
/home/kongf/workhome/sawtooth/moosegpu/petsc/src/sys/objects/pinit.c:1238
[0]PETSC ERROR: #6 SlepcInitialize() at
/home/kongf/workhome/sawtooth/moosegpu/petsc/arch-linux-c-opt/externalpackages/git.slepc/src/sys/slepcinit.c:275
[0]PETSC ERROR: #7 LibMeshInit() at ../src/base/libmesh.C:522
[r8i3n0:mpi_rank_0][MPIDI_CH3_Abort] application called
MPI_Abort(MPI_COMM_WORLD, 95) - process 0: No such file or directory (2)

Thanks,

Fande


Re: [petsc-users] PETSc MUMPS interface

2022-01-19 Thread Zhang, Hong via petsc-users
Varun,
Good to know it works. FactorSymbolic function is still being called twice, but 
the 2nd call is a no-op, thus it still appears in '-log_view'. I made changes 
in the low level of mumps routine, not within PCSetUp() because I feel your use 
case is limited to mumps, not other matrix package solvers.
Hong

From: Varun Hiremath 
Sent: Wednesday, January 19, 2022 2:44 AM
To: Zhang, Hong 
Cc: Peder Jørgensgaard Olesen via petsc-users 
Subject: Re: [petsc-users] PETSc MUMPS interface

Hi Hong,

Thanks, I tested your branch and I think it is working fine. I don't see any 
increase in runtime, however with -log_view I see that the MatLUFactorSymbolic 
function is still being called twice, so is this expected? Is the second call a 
no-op?

$ ./ex52.o -use_mumps_lu -print_mumps_memory -log_view | grep MatLUFactorSym
MatLUFactorSym 2 1.0 4.4411e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 
0.0e+00  2  0  0  0  0   2  0  0  0  0 0

Thanks,
Varun

On Mon, Jan 17, 2022 at 7:49 PM Zhang, Hong 
mailto:hzh...@mcs.anl.gov>> wrote:
Varun,
I created a branch hzhang/feature-mumps-mem-estimate,
see https://gitlab.com/petsc/petsc/-/merge_requests/4727

You may give it a try and let me know if this is what you want.
src/ksp/ksp/tutorials/ex52.c is an example.

Hong

From: Varun Hiremath mailto:varunhirem...@gmail.com>>
Sent: Monday, January 17, 2022 12:41 PM
To: Zhang, Hong mailto:hzh...@mcs.anl.gov>>
Cc: Jose E. Roman mailto:jro...@dsic.upv.es>>; Peder 
Jørgensgaard Olesen via petsc-users 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] PETSc MUMPS interface

Hi Hong,

Thanks for looking into this. Here is the workflow that I might use:

MatLUFactorSymbolic(F,A,perm,iperm,);

// get memory estimates from MUMPS e.g. INFO(3), INFOG(16), INFOG(17)
// find available memory on the system e.g. RAM size
if (estimated_memory > available_memory)
{
   // inform and stop; or
   // switch MUMPS to out-of-core factorization
   ICNTL(22) = 1;
}
else
{
   // set appropriate settings for in-core factorization
}

// Now we call the solve and inside if MatLUFactorSymbolic is already called 
then it should be skipped
EPSSolve(eps);

Thanks,
Varun

On Mon, Jan 17, 2022 at 9:18 AM Zhang, Hong 
mailto:hzh...@mcs.anl.gov>> wrote:
Varun,
I am trying to find a way to enable you to switch options after 
MatLUFactorSymbolic(). A hack is modifying the flag  'mumps->matstruc' inside  
MatLUFactorSymbolic_AIJMUMPS() and MatFactorNumeric_MUMPS().

My understanding of what you want is:
  // collect mumps memory info
  ...
  MatLUFactorSymbolic(F,A,perm,iperm,);
  printMumpsMemoryInfo(F);
  //-
if (memory is available) {
EPSSolve(eps); --> skip calling of MatLUFactorSymbolic()
} else {
   //out-of-core (OOC) option in MUMPS
}

Am I correct? I'll let you know once I work out a solution.
Hong


From: Varun Hiremath mailto:varunhirem...@gmail.com>>
Sent: Sunday, January 16, 2022 10:10 PM
To: Zhang, Hong mailto:hzh...@mcs.anl.gov>>
Cc: Jose E. Roman mailto:jro...@dsic.upv.es>>; Peder 
Jørgensgaard Olesen via petsc-users 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] PETSc MUMPS interface

Hi Jose, Hong,

Thanks for the explanation. I have verified using -log_view that 
MatLUFactorSymbolic is indeed getting called twice.

Hong, we use MUMPS solver for other things, and we typically run the symbolic 
analysis first and get memory estimates to ensure that we have enough memory 
available to run the case. If the available memory is not enough, we can stop 
or switch to the out-of-core (OOC) option in MUMPS. We wanted to do the same 
when using MUMPS via SLEPc/PETSc. Please let me know if there are other ways of 
getting these memory stats and switching options during runtime with PETSc.
Appreciate your help!

Thanks,
Varun

On Sun, Jan 16, 2022 at 4:01 PM Zhang, Hong 
mailto:hzh...@mcs.anl.gov>> wrote:
Varun,
I believe Jose is correct. You may verify it by running your code with option 
'-log_view', then check the number of calls to MatLUFactorSym.

I guess I can add a flag in PCSetUp() to check if user has already called 
MatLUFactorSymbolic() and wants to skip it. Normally, users simply allocate 
sufficient memory in the symbolic factorization. Why do you want to check it?
Hong


From: Jose E. Roman mailto:jro...@dsic.upv.es>>
Sent: Sunday, January 16, 2022 5:11 AM
To: Varun Hiremath mailto:varunhirem...@gmail.com>>
Cc: Zhang, Hong mailto:hzh...@mcs.anl.gov>>; Peder 
Jørgensgaard Olesen via petsc-users 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] PETSc MUMPS interface

Hong may give a better answer, but if you look at PCSetUp_LU() 
https://petsc.org/main/src/ksp/pc/impls/factor/lu/lu.c.html#PCSetUp_LU you will 
see that MatLUFactorSymbolic() is called unconditionally during the first 
PCSetUp(). Currently there is no way to check if the user has already called 

Re: [petsc-users] Downloaded superlu_dist could not be used. Please check install in $PREFIX

2022-01-19 Thread Fande Kong
Thanks, Sherry, and Satish,

I will try your suggestion, and report back to you as soon as possible.

Thanks,

Fande

On Tue, Jan 18, 2022 at 10:48 PM Satish Balay  wrote:

> Sherry,
>
> This is with superlu-dist-7.1.1 [not master branch]
>
>
> Fande,
>
> >>
> Executing: mpifort  -o /tmp/petsc-UYa6A8/config.compilers/conftest
> -fopenmp -fopenmp  -I$PREFIX/include -fPIC -O3  -fopenmp
> /tmp/petsc-UYa6A8/config.compilers/conftest.o
> /tmp/petsc-UYa6A8/config.compilers/confc.o  -Wl,-rpath,$PREFIX/lib
> -L$PREFIX/lib -lsuperlu_dist -lpthread -Wl,-rpath,$PREFIX/lib -L$PREFIX/lib
> -lparmetis -Wl,-rpath,$PREFIX/lib -L$PREFIX/lib -lmetis
> -Wl,-rpath,$PREFIX/lib -L$PREFIX/lib -lflapack -Wl,-rpath,$PREFIX/lib
> -L$PREFIX/lib -lfblas -lm -Wl,-rpath,$BUILD_PREFIX/lib -L$BUILD_PREFIX/lib
> -lstdc++ -ldl -lmpifort -lmpi -lgfortran -lm
> -Wl,-rpath,$BUILD_PREFIX/lib/gcc/x86_64-conda-linux-gnu/9.3.0
> -L$BUILD_PREFIX/lib/gcc/x86_64-conda-linux-gnu/9.3.0
> -Wl,-rpath,$BUILD_PREFIX/lib/gcc -L$BUILD_PREFIX/lib/gcc
> -Wl,-rpath,$BUILD_PREFIX/x86_64-conda-linux-gnu/lib
> -L$BUILD_PREFIX/x86_64-conda-linux-gnu/lib -Wl,-rpath,$BUILD_PREFIX/lib
> -lgfortran -lm -lgcc_s -lquadmath -lrt -lquadmath -lstdc++ -ldl
> Possible ERROR while running linker:
> stderr:
> $BUILD_PREFIX/bin/../lib/gcc/x86_64-conda-linux-gnu/9.3.0/../../../../x86_64-conda-linux-gnu/bin/ld:
> warning: libmpicxx.so.12, needed by $PREFIX/lib/libsuperlu_dist.so, not
> found (try using -rpath or -rpath-link)
> <<<
>
> I don't really understand why this error comes up [as with shared
> libraries we should be able to link with -lsuperlu_dist - without having to
> link with libmpicxx.so.12
>
> What do you get for:
>
> ldd $PREFIX/lib/libstdc++.so
>
>
> BTW: is configure.log modified to replace realpaths with $PREFIX
> $BUILD_PREFIX etc?
>
> Can you try additional configure option LIBS=-lmpicxx and see if that
> works around this problem?
>
> Satish
>
> On Tue, 18 Jan 2022, Xiaoye S. Li wrote:
>
> > There was a merge error in the master branch. I fixed it today. Not sure
> > whether that's causing your problem.   Can you try now?
> >
> > Sherry
> >
> > On Mon, Jan 17, 2022 at 11:55 AM Fande Kong  wrote:
> >
> > > I am trying to port PETSc-3.16.3 to the MOOSE ecosystem. I got an error
> > > that PETSc could not build  superlu_dist.  The log file was attached.
> > >
> > > PETSc-3.15.x worked correctly in the same environment.
> > >
> > > Thanks,
> > > Fande
> > >
> >
>
>


Re: [petsc-users] Nullspaces

2022-01-19 Thread Mark Adams
ILU is LU for this 1D problem and its singular.
ILU might have some logic to deal with a singular system. Not sure. LU
should fail.
You might try -ksp_type preonly and -pc_type lu
And -pc_type jacobi should not have any numerical problems. Try that.

On Wed, Jan 19, 2022 at 7:19 AM Matthew Knepley  wrote:

> On Wed, Jan 19, 2022 at 4:52 AM Marco Cisternino <
> marco.cistern...@optimad.it> wrote:
>
>> Thank you Matthew.
>>
>> But I cannot get the point. I got the point about the test but to try to
>> explain my doubt I’m going to prepare another toy code.
>>
>>
>>
>> By words…
>>
>> I usually have a finite volume discretization of the Laplace operator
>> with homogeneous Neumann BC on an octree mesh and it reads
>>
>> Aij * xj = bi,
>>
>> being the discretization of
>>
>> Int|Vi(nabla^2 pi dV) = Int|Vi(nabla dot ui)
>>
>> (I omit constants), where Int|Vi(…dV) is a volume integral on the I cell,
>> pi is cell pressure, ui is the cell velocity.
>>
>> The computational domain contains 2 separated sub-domains.
>>
>> Let’s consider 4 cells into the whole domain and 2 cells for each
>> sub-domain.
>>
>> I would expect a null space made of 2 vectors and from your patch they
>> should look like [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)],
>> i.e. norm2 = 1 for both.
>>
>> With MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr,
>> ) I’m saying that [1/sqrt(4) 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)] is
>> the null space, which is not but it is in the null space.
>>
>> But this is not the case I sent you, not exactly.
>>
>>
>>
>> The case I sent is 1/Vi * Aij * xj = 1/Vi bi, where Vi is the volume of
>> the cell i.
>>
>> Let’s admit that yj is in the null space of of Aij, it should be in the
>> null space of 1/Vi * Aij, therefore Aij*yj = 0 and 1/Vi * Aij*yj = 0 too.
>>
>> But in the framework of the test, this is true with infinite precision.
>>
>> What happens if norm2(Aij*yj) = 10^-15 and Vi = 10^-5? Norm2(1/Vi * Aij
>> * yj) = 10^-10!!! Is yi still in the null space numerically?
>>
>> Let’s say yi is the constant vector over the whole domain, i.e.
>> [1/sqrt(4) 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)]. Should this be in the null space
>> of 1/Vi * Aij, shouldn’t it?
>>
>> An analogous argument should be for the compatibility condition that
>> concerns bi.
>>
>>
>>
>> My current problem is that properly creating the null space for Aij, i.e.
>> [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)], allows me to solve
>> and find xi, but multiplying by 1/Vi, I cannot get any solution using both
>> FGMRES+ILU and FGMRE+GAMG.
>>
>>
>>
>> The tiny problem will load Aij, Vi and bi and show the problem by testing
>> the proper null space and trying to solve. I will include the patch to my
>> PETSc version. I hope to come back to you very soon.
>>
>
> This sounds like a dimensionalization problem to me. It is best to choose
> length (and other) units that make the matrix entries about 1. It seems
> like you are far from this, and it is causing a loss of accuracy (your
> null vector has a residual of about 1e-7).
>
>   Thanks,
>
>   Matt
>
>
>> Thank you very much for your support!
>>
>>
>>
>>
>>
>> Marco Cisternino
>>
>>
>>
>> *From:* Matthew Knepley 
>> *Sent:* martedì 18 gennaio 2022 21:25
>> *To:* Marco Cisternino 
>> *Cc:* petsc-users 
>> *Subject:* Re: [petsc-users] Nullspaces
>>
>>
>>
>> I made a fix for this:
>>
>>
>>
>>   https://gitlab.com/petsc/petsc/-/merge_requests/4729
>>
>>
>>
>>   Thanks,
>>
>>
>>
>>  Matt
>>
>>
>>
>> On Tue, Jan 18, 2022 at 3:20 PM Matthew Knepley 
>> wrote:
>>
>> On Thu, Dec 16, 2021 at 11:09 AM Marco Cisternino <
>> marco.cistern...@optimad.it> wrote:
>>
>> Hello Matthew,
>>
>> as promised I prepared a minimal (112960 rows. I’m not able to produce
>> anything smaller than this and triggering the issue) example of the
>> behavior I was talking about some days ago.
>>
>> What I did is to produce matrix, right hand side and initial solution of
>> the linear system.
>>
>>
>>
>> As I told you before, this linear system is the discretization of the
>> pressure equation of a predictor-corrector method for NS equations in the
>> framework of finite volume method.
>>
>> This case has homogeneous Neumann boundary conditions. Computational
>> domain has two independent and separated sub-domains.
>>
>> I discretize the weak formulation and I divide every row of the linear
>> system by the volume of the relative cell.
>>
>> The underlying mesh is not uniform, therefore cells have different
>> volumes.
>>
>> The issue I’m going to explain does not show up if the mesh is uniform,
>> same volume for all the cells.
>>
>>
>>
>> I usually build the null space sub-domain by sub-domain with
>>
>> MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, nConstants, constants,
>> );
>>
>> Where nConstants = 2 and constants contains two normalized arrays with
>> constant values on degrees of freedom relative to the associated sub-domain
>> and zeros elsewhere.
>>
>>
>>
>> However, as a test I tried the 

Re: [petsc-users] Nullspaces

2022-01-19 Thread Matthew Knepley
On Wed, Jan 19, 2022 at 4:52 AM Marco Cisternino <
marco.cistern...@optimad.it> wrote:

> Thank you Matthew.
>
> But I cannot get the point. I got the point about the test but to try to
> explain my doubt I’m going to prepare another toy code.
>
>
>
> By words…
>
> I usually have a finite volume discretization of the Laplace operator with
> homogeneous Neumann BC on an octree mesh and it reads
>
> Aij * xj = bi,
>
> being the discretization of
>
> Int|Vi(nabla^2 pi dV) = Int|Vi(nabla dot ui)
>
> (I omit constants), where Int|Vi(…dV) is a volume integral on the I cell,
> pi is cell pressure, ui is the cell velocity.
>
> The computational domain contains 2 separated sub-domains.
>
> Let’s consider 4 cells into the whole domain and 2 cells for each
> sub-domain.
>
> I would expect a null space made of 2 vectors and from your patch they
> should look like [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)],
> i.e. norm2 = 1 for both.
>
> With MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr,
> ) I’m saying that [1/sqrt(4) 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)] is
> the null space, which is not but it is in the null space.
>
> But this is not the case I sent you, not exactly.
>
>
>
> The case I sent is 1/Vi * Aij * xj = 1/Vi bi, where Vi is the volume of
> the cell i.
>
> Let’s admit that yj is in the null space of of Aij, it should be in the
> null space of 1/Vi * Aij, therefore Aij*yj = 0 and 1/Vi * Aij*yj = 0 too.
>
> But in the framework of the test, this is true with infinite precision.
>
> What happens if norm2(Aij*yj) = 10^-15 and Vi = 10^-5? Norm2(1/Vi * Aij *
> yj) = 10^-10!!! Is yi still in the null space numerically?
>
> Let’s say yi is the constant vector over the whole domain, i.e. [1/sqrt(4)
> 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)]. Should this be in the null space of 1/Vi *
> Aij, shouldn’t it?
>
> An analogous argument should be for the compatibility condition that
> concerns bi.
>
>
>
> My current problem is that properly creating the null space for Aij, i.e.
> [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)], allows me to solve
> and find xi, but multiplying by 1/Vi, I cannot get any solution using both
> FGMRES+ILU and FGMRE+GAMG.
>
>
>
> The tiny problem will load Aij, Vi and bi and show the problem by testing
> the proper null space and trying to solve. I will include the patch to my
> PETSc version. I hope to come back to you very soon.
>

This sounds like a dimensionalization problem to me. It is best to choose
length (and other) units that make the matrix entries about 1. It seems
like you are far from this, and it is causing a loss of accuracy (your null
vector has a residual of about 1e-7).

  Thanks,

  Matt


> Thank you very much for your support!
>
>
>
>
>
> Marco Cisternino
>
>
>
> *From:* Matthew Knepley 
> *Sent:* martedì 18 gennaio 2022 21:25
> *To:* Marco Cisternino 
> *Cc:* petsc-users 
> *Subject:* Re: [petsc-users] Nullspaces
>
>
>
> I made a fix for this:
>
>
>
>   https://gitlab.com/petsc/petsc/-/merge_requests/4729
>
>
>
>   Thanks,
>
>
>
>  Matt
>
>
>
> On Tue, Jan 18, 2022 at 3:20 PM Matthew Knepley  wrote:
>
> On Thu, Dec 16, 2021 at 11:09 AM Marco Cisternino <
> marco.cistern...@optimad.it> wrote:
>
> Hello Matthew,
>
> as promised I prepared a minimal (112960 rows. I’m not able to produce
> anything smaller than this and triggering the issue) example of the
> behavior I was talking about some days ago.
>
> What I did is to produce matrix, right hand side and initial solution of
> the linear system.
>
>
>
> As I told you before, this linear system is the discretization of the
> pressure equation of a predictor-corrector method for NS equations in the
> framework of finite volume method.
>
> This case has homogeneous Neumann boundary conditions. Computational
> domain has two independent and separated sub-domains.
>
> I discretize the weak formulation and I divide every row of the linear
> system by the volume of the relative cell.
>
> The underlying mesh is not uniform, therefore cells have different
> volumes.
>
> The issue I’m going to explain does not show up if the mesh is uniform,
> same volume for all the cells.
>
>
>
> I usually build the null space sub-domain by sub-domain with
>
> MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, nConstants, constants,
> );
>
> Where nConstants = 2 and constants contains two normalized arrays with
> constant values on degrees of freedom relative to the associated sub-domain
> and zeros elsewhere.
>
>
>
> However, as a test I tried the constant over the whole domain using 2
> alternatives that should produce the same null space:
>
>1. MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr,
>);
>2. Vec* nsp;
>
> VecDuplicateVecs(solution, 1, );
>
> VecSet(nsp[0],1.0);
>
> VecNormalize(nsp[0], nullptr);
>
> MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, 1, nsp, );
>
>
>
> Once I created the null space I test it using:
>
> MatNullSpaceTest(nullspace, m_A, );
>
>
>
> The case 1 pass the test 

Re: [petsc-users] Nullspaces

2022-01-19 Thread Marco Cisternino
Thank you Matthew.
But I cannot get the point. I got the point about the test but to try to 
explain my doubt I’m going to prepare another toy code.

By words…
I usually have a finite volume discretization of the Laplace operator with 
homogeneous Neumann BC on an octree mesh and it reads
Aij * xj = bi,
being the discretization of
Int|Vi(nabla^2 pi dV) = Int|Vi(nabla dot ui)
(I omit constants), where Int|Vi(…dV) is a volume integral on the I cell, pi is 
cell pressure, ui is the cell velocity.
The computational domain contains 2 separated sub-domains.
Let’s consider 4 cells into the whole domain and 2 cells for each sub-domain.
I would expect a null space made of 2 vectors and from your patch they should 
look like [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)], i.e. norm2 = 
1 for both.
With MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr, ) 
I’m saying that [1/sqrt(4) 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)] is the null space, 
which is not but it is in the null space.
But this is not the case I sent you, not exactly.

The case I sent is 1/Vi * Aij * xj = 1/Vi bi, where Vi is the volume of the 
cell i.
Let’s admit that yj is in the null space of of Aij, it should be in the null 
space of 1/Vi * Aij, therefore Aij*yj = 0 and 1/Vi * Aij*yj = 0 too.
But in the framework of the test, this is true with infinite precision.
What happens if norm2(Aij*yj) = 10^-15 and Vi = 10^-5? Norm2(1/Vi * Aij * yj) = 
10^-10!!! Is yi still in the null space numerically?
Let’s say yi is the constant vector over the whole domain, i.e. [1/sqrt(4) 
1/sqrt(4) 1/sqrt(4) 1/sqrt(4)]. Should this be in the null space of 1/Vi * Aij, 
shouldn’t it?
An analogous argument should be for the compatibility condition that concerns 
bi.

My current problem is that properly creating the null space for Aij, i.e. 
[1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)], allows me to solve and 
find xi, but multiplying by 1/Vi, I cannot get any solution using both 
FGMRES+ILU and FGMRE+GAMG.

The tiny problem will load Aij, Vi and bi and show the problem by testing the 
proper null space and trying to solve. I will include the patch to my PETSc 
version. I hope to come back to you very soon.

Thank you very much for your support!


Marco Cisternino

From: Matthew Knepley 
Sent: martedì 18 gennaio 2022 21:25
To: Marco Cisternino 
Cc: petsc-users 
Subject: Re: [petsc-users] Nullspaces

I made a fix for this:

  https://gitlab.com/petsc/petsc/-/merge_requests/4729

  Thanks,

 Matt

On Tue, Jan 18, 2022 at 3:20 PM Matthew Knepley 
mailto:knep...@gmail.com>> wrote:
On Thu, Dec 16, 2021 at 11:09 AM Marco Cisternino 
mailto:marco.cistern...@optimad.it>> wrote:
Hello Matthew,
as promised I prepared a minimal (112960 rows. I’m not able to produce anything 
smaller than this and triggering the issue) example of the behavior I was 
talking about some days ago.
What I did is to produce matrix, right hand side and initial solution of the 
linear system.

As I told you before, this linear system is the discretization of the pressure 
equation of a predictor-corrector method for NS equations in the framework of 
finite volume method.
This case has homogeneous Neumann boundary conditions. Computational domain has 
two independent and separated sub-domains.
I discretize the weak formulation and I divide every row of the linear system 
by the volume of the relative cell.
The underlying mesh is not uniform, therefore cells have different volumes.
The issue I’m going to explain does not show up if the mesh is uniform, same 
volume for all the cells.

I usually build the null space sub-domain by sub-domain with
MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, nConstants, constants, 
);
Where nConstants = 2 and constants contains two normalized arrays with constant 
values on degrees of freedom relative to the associated sub-domain and zeros 
elsewhere.

However, as a test I tried the constant over the whole domain using 2 
alternatives that should produce the same null space:

  1.  MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr, );
  2.  Vec* nsp;

VecDuplicateVecs(solution, 1, );

VecSet(nsp[0],1.0);

VecNormalize(nsp[0], nullptr);

MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, 1, nsp, );


Once I created the null space I test it using:
MatNullSpaceTest(nullspace, m_A, );

The case 1 pass the test while case 2 don’t.

I have a small code for matrix loading, null spaces creation and testing.
Unfortunately I cannot implement a small code able to produce that linear 
system.

As attachment you can find an archive containing the matrix, the initial 
solution (used to manually build the null space) and the rhs (not used in the 
test code) in binary format.
You can also find the testing code in the same archive.
I used petsc 3.12(gcc+openMPI) and petsc 3.15.2(intelOneAPI) same results.
If the attachment is not delivered, I can share a link to it.

Marco, please forgive me for taking so long to get to your issue. It has been 
crazy.


Re: [petsc-users] PETSc MUMPS interface

2022-01-19 Thread Varun Hiremath
Hi Hong,

Thanks, I tested your branch and I think it is working fine. I don't see
any increase in runtime, however with -log_view I see that the
MatLUFactorSymbolic function is still being called twice, so is this
expected? Is the second call a no-op?

$ ./ex52.o -use_mumps_lu -print_mumps_memory -log_view | grep MatLUFactorSym
MatLUFactorSym 2 1.0 4.4411e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00  2  0  0  0  0   2  0  0  0  0 0

Thanks,
Varun

On Mon, Jan 17, 2022 at 7:49 PM Zhang, Hong  wrote:

> Varun,
> I created a branch hzhang/feature-mumps-mem-estimate,
> see https://gitlab.com/petsc/petsc/-/merge_requests/4727
>
> You may give it a try and let me know if this is what you want.
> src/ksp/ksp/tutorials/ex52.c is an example.
>
> Hong
> --
> *From:* Varun Hiremath 
> *Sent:* Monday, January 17, 2022 12:41 PM
> *To:* Zhang, Hong 
> *Cc:* Jose E. Roman ; Peder Jørgensgaard Olesen via
> petsc-users 
> *Subject:* Re: [petsc-users] PETSc MUMPS interface
>
> Hi Hong,
>
> Thanks for looking into this. Here is the workflow that I might use:
>
> MatLUFactorSymbolic(F,A,perm,iperm,);
>
> // get memory estimates from MUMPS e.g. INFO(3), INFOG(16), INFOG(17)
> // find available memory on the system e.g. RAM size
> if (estimated_memory > available_memory)
> {
>// inform and stop; or
>// switch MUMPS to out-of-core factorization
>ICNTL(22) = 1;
> }
> else
> {
>// set appropriate settings for in-core factorization
> }
>
> // Now we call the solve and inside if MatLUFactorSymbolic is already
> called then it should be skipped
> EPSSolve(eps);
>
> Thanks,
> Varun
>
> On Mon, Jan 17, 2022 at 9:18 AM Zhang, Hong  wrote:
>
> Varun,
> I am trying to find a way to enable you to switch options after 
> MatLUFactorSymbolic().
> A hack is modifying the flag  'mumps->matstruc'
> inside  MatLUFactorSymbolic_AIJMUMPS() and MatFactorNumeric_MUMPS().
>
> My understanding of what you want is:
>   // collect mumps memory info
>   ...
>   MatLUFactorSymbolic(F,A,perm,iperm,);
>   printMumpsMemoryInfo(F);
>   //-
> if (memory is available) {
> EPSSolve(eps); --> skip calling of MatLUFactorSymbolic()
> } else {
>//out-of-core (OOC) option in MUMPS
> }
>
> Am I correct? I'll let you know once I work out a solution.
> Hong
>
> --
> *From:* Varun Hiremath 
> *Sent:* Sunday, January 16, 2022 10:10 PM
> *To:* Zhang, Hong 
> *Cc:* Jose E. Roman ; Peder Jørgensgaard Olesen via
> petsc-users 
> *Subject:* Re: [petsc-users] PETSc MUMPS interface
>
> Hi Jose, Hong,
>
> Thanks for the explanation. I have verified using -log_view that 
> MatLUFactorSymbolic
> is indeed getting called twice.
>
> Hong, we use MUMPS solver for other things, and we typically run the
> symbolic analysis first and get memory estimates to ensure that we have
> enough memory available to run the case. If the available memory is not
> enough, we can stop or switch to the out-of-core (OOC) option in MUMPS. We
> wanted to do the same when using MUMPS via SLEPc/PETSc. Please let me know
> if there are other ways of getting these memory stats and switching options
> during runtime with PETSc.
> Appreciate your help!
>
> Thanks,
> Varun
>
> On Sun, Jan 16, 2022 at 4:01 PM Zhang, Hong  wrote:
>
> Varun,
> I believe Jose is correct. You may verify it by running your code with
> option '-log_view', then check the number of calls to MatLUFactorSym.
>
> I guess I can add a flag in PCSetUp() to check if user has already called
> MatLUFactorSymbolic() and wants to skip it. Normally, users simply allocate
> sufficient memory in the symbolic factorization. Why do you want to check
> it?
> Hong
>
> --
> *From:* Jose E. Roman 
> *Sent:* Sunday, January 16, 2022 5:11 AM
> *To:* Varun Hiremath 
> *Cc:* Zhang, Hong ; Peder Jørgensgaard Olesen via
> petsc-users 
> *Subject:* Re: [petsc-users] PETSc MUMPS interface
>
> Hong may give a better answer, but if you look at PCSetUp_LU()
> https://petsc.org/main/src/ksp/pc/impls/factor/lu/lu.c.html#PCSetUp_LU
> you will see that MatLUFactorSymbolic() is called unconditionally during
> the first PCSetUp(). Currently there is no way to check if the user has
> already called MatLUFactorSymbolic().
>
> Jose
>
>
> > El 16 ene 2022, a las 10:40, Varun Hiremath 
> escribió:
> >
> > Hi Hong,
> >
> > Thank you, this is very helpful!
> >
> > Using this method I am able to get the memory estimates before the
> actual solve, however, I think my code may be causing the symbolic
> factorization to be run twice. Attached is my code where I am using SLEPc
> to compute eigenvalues, and I use MUMPS for factorization. I have commented
> above the code that computes the memory estimates, could you please check
> and tell me if this would cause the symbolic factor to be computed twice (a
> second time inside EPSSolve?), as I am seeing a slight increase in the
> overall computation time?
> >
> > Regards,
> > Varun
> >
> > On Wed, Jan 12, 2022