from:"Danyang Su"

Hi Junchao,

 

I cannot configure your branch with same options due to the error in sowing. I 
had similar error before on other clusters with very old openmpi version. 
Problem was solved when openmpi was updated to a newer one. 

 

At this moment, I configured a PETSc version with Openmpi 2.1.6 version and it 
seems working properly. 

 

Thanks and have a good rest of the weekend,

 

Danyang

 

From: Danyang Su 
Date: Saturday, April 10, 2021 at 4:08 PM
To: Junchao Zhang 
Cc: Barry Smith , "petsc-users@mcs.anl.gov" 

Subject: Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI 
version

 

Hi Junchao,

 

The configuration is successful. The error comes from the last step when I run 

 

make PETSC_DIR=/home/danyangs/soft/petsc/petsc-3.13.6 
PETSC_ARCH=linux-intel-openmpi check

 

***Error detected during compile or link!***

See http://www.mcs.anl.gov/petsc/documentation/faq.html

/home/danyangs/soft/petsc/petsc-3.13.6/src/snes/tutorials ex5f

*

mpif90 -fPIC -O3 -march=native -mtune=nativels   
-I/home/danyangs/soft/petsc/petsc-3.13.6/include 
-I/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/include  
ex5f.F90  
-Wl,-rpath,/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib 
-L/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib 
-Wl,-rpath,/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib 
-L/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib 
-Wl,-rpath,/global/software/intel/composer_xe_2013_sp1.2.144/mkl/lib/intel64 
-L/global/software/intel/composer_xe_2013_sp1.2.144/mkl/lib/intel64 
-Wl,-rpath,/global/software/intel/composer_xe_2013_sp1.2.144/compiler/lib/intel64
 -L/global/software/intel/composer_xe_2013_sp1.2.144/compiler/lib/intel64 
-Wl,-rpath,/global/software/openmpi-1.6.5/intel/lib64 
-L/global/software/openmpi-1.6.5/intel/lib64 
-Wl,-rpath,/global/software/intel/composerxe/mkl/lib/intel64 
-L/global/software/intel/composerxe/mkl/lib/intel64 
-Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/4.4.7 
-L/usr/lib/gcc/x86_64-redhat-linux/4.4.7 
-Wl,-rpath,/global/software/intel/composerxe/lib/intel64 -lpetsc -lHYPRE 
-lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lscalapack -lsuperlu 
-lflapack -lfblas -lX11 -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5 
-lparmetis -lmetis -lstdc++ -ldl -lmpi_f90 -lmpi_f77 -lmpi -lm -lnuma -lrt 
-lnsl -lutil -limf -lifport -lifcore -lsvml -lipgo -lintlc -lpthread -lgcc_s 
-lirc_s -lstdc++ -ldl -o ex5f

ifort: command line warning #10159: invalid argument for option '-m'

/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_allocate'

/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_attach'

/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_create_dynamic'

gmake[4]: *** [ex5f] Error 1

 

Thanks,

 

Danyang

 

From: Junchao Zhang 
Date: Saturday, April 10, 2021 at 3:57 PM
To: Danyang Su 
Cc: Barry Smith , "petsc-users@mcs.anl.gov" 

Subject: Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI 
version

 

You sent a wrong one. This configure.log was from a successful configuration. 
Note FOPTFLAGS="-O3 -march=native -mtune=nativels" looks suspicious.

 

--Junchao Zhang

 

 

On Sat, Apr 10, 2021 at 5:32 PM Danyang Su  wrote:

 

Hi Junchao,

 

Thanks for looking into this problem. The configuration log is attached.

 

All the best,

 

Danyang

From: Junchao Zhang 
Date: Saturday, April 10, 2021 at 2:36 PM
To: Danyang Su 
Cc: Barry Smith , "petsc-users@mcs.anl.gov" 

Subject: Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI 
version

 

Hi, Danyang, 

 

Send the configure.log.  Also, PETSc does not need MPI_Win_allocate etc to 
work. I will have a look.


--Junchao Zhang

 

 

On Sat, Apr 10, 2021 at 2:47 PM Danyang Su  wrote:

Hi Barry,

 

I tried this option before but get ‘Error running configure on OpenMPI’

 

***

 UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for 
details):

---

Error running configure on OPENMPI

***

  File "/global/home/danyangs/soft/petsc/petsc-3.14.6/config/configure.py", 
line 456, in petsc_configure

framework.configure(out = sys.stdout)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1253, in configure

self.processChildren()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1242, in processChildren

self.

Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI version

Hi Junchao,

 

The configuration is successful. The error comes from the last step when I run 

 

make PETSC_DIR=/home/danyangs/soft/petsc/petsc-3.13.6 
PETSC_ARCH=linux-intel-openmpi check

 

***Error detected during compile or link!***

See http://www.mcs.anl.gov/petsc/documentation/faq.html

/home/danyangs/soft/petsc/petsc-3.13.6/src/snes/tutorials ex5f

*

mpif90 -fPIC -O3 -march=native -mtune=nativels   
-I/home/danyangs/soft/petsc/petsc-3.13.6/include 
-I/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/include  
ex5f.F90  
-Wl,-rpath,/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib 
-L/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib 
-Wl,-rpath,/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib 
-L/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib 
-Wl,-rpath,/global/software/intel/composer_xe_2013_sp1.2.144/mkl/lib/intel64 
-L/global/software/intel/composer_xe_2013_sp1.2.144/mkl/lib/intel64 
-Wl,-rpath,/global/software/intel/composer_xe_2013_sp1.2.144/compiler/lib/intel64
 -L/global/software/intel/composer_xe_2013_sp1.2.144/compiler/lib/intel64 
-Wl,-rpath,/global/software/openmpi-1.6.5/intel/lib64 
-L/global/software/openmpi-1.6.5/intel/lib64 
-Wl,-rpath,/global/software/intel/composerxe/mkl/lib/intel64 
-L/global/software/intel/composerxe/mkl/lib/intel64 
-Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/4.4.7 
-L/usr/lib/gcc/x86_64-redhat-linux/4.4.7 
-Wl,-rpath,/global/software/intel/composerxe/lib/intel64 -lpetsc -lHYPRE 
-lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lscalapack -lsuperlu 
-lflapack -lfblas -lX11 -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5 
-lparmetis -lmetis -lstdc++ -ldl -lmpi_f90 -lmpi_f77 -lmpi -lm -lnuma -lrt 
-lnsl -lutil -limf -lifport -lifcore -lsvml -lipgo -lintlc -lpthread -lgcc_s 
-lirc_s -lstdc++ -ldl -o ex5f

ifort: command line warning #10159: invalid argument for option '-m'

/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_allocate'

/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_attach'

/home/danyangs/soft/petsc/petsc-3.13.6/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_create_dynamic'

gmake[4]: *** [ex5f] Error 1

 

Thanks,

 

Danyang

 

From: Junchao Zhang 
Date: Saturday, April 10, 2021 at 3:57 PM
To: Danyang Su 
Cc: Barry Smith , "petsc-users@mcs.anl.gov" 

Subject: Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI 
version

 

You sent a wrong one. This configure.log was from a successful configuration. 
Note FOPTFLAGS="-O3 -march=native -mtune=nativels" looks suspicious.

 

--Junchao Zhang

 

 

On Sat, Apr 10, 2021 at 5:32 PM Danyang Su  wrote:

 

Hi Junchao,

 

Thanks for looking into this problem. The configuration log is attached.

 

All the best,

 

Danyang

From: Junchao Zhang 
Date: Saturday, April 10, 2021 at 2:36 PM
To: Danyang Su 
Cc: Barry Smith , "petsc-users@mcs.anl.gov" 

Subject: Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI 
version

 

Hi, Danyang, 

 

Send the configure.log.  Also, PETSc does not need MPI_Win_allocate etc to 
work. I will have a look.


--Junchao Zhang

 

 

On Sat, Apr 10, 2021 at 2:47 PM Danyang Su  wrote:

Hi Barry,

 

I tried this option before but get ‘Error running configure on OpenMPI’

 

***

 UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for 
details):

---

Error running configure on OPENMPI

***

  File "/global/home/danyangs/soft/petsc/petsc-3.14.6/config/configure.py", 
line 456, in petsc_configure

framework.configure(out = sys.stdout)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1253, in configure

self.processChildren()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1242, in processChildren

self.serialEvaluation(self.childGraph)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1217, in serialEvaluation

child.configure()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 1144, in configure

self.executeTest(self.configureLibrary)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/base.py",
 line 140, in executeTest

ret = test(*args,**kargs)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py&q

Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI version

Hi Junchao,

 

Thanks. I will test this branch and get back to you later.

 

All the best,

 

Danyang

 

From: Junchao Zhang 
Date: Saturday, April 10, 2021 at 3:32 PM
To: Danyang Su 
Cc: Barry Smith , "petsc-users@mcs.anl.gov" 

Subject: Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI 
version

 

Danyang,

  Could you try branch jczhang/fix-mpi3-win with your old configuration (i.e., 
use system mpicc)? Note the MR 3849 is based off latest petsc-3.15 release

  Thanks.

--Junchao Zhang

 

 

On Sat, Apr 10, 2021 at 4:36 PM Junchao Zhang  wrote:

Hi, Danyang, 

 

Send the configure.log.  Also, PETSc does not need MPI_Win_allocate etc to 
work. I will have a look.


--Junchao Zhang

 

 

On Sat, Apr 10, 2021 at 2:47 PM Danyang Su  wrote:

Hi Barry,

 

I tried this option before but get ‘Error running configure on OpenMPI’

 

***

 UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for 
details):

---

Error running configure on OPENMPI

***

  File "/global/home/danyangs/soft/petsc/petsc-3.14.6/config/configure.py", 
line 456, in petsc_configure

framework.configure(out = sys.stdout)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1253, in configure

self.processChildren()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1242, in processChildren

self.serialEvaluation(self.childGraph)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1217, in serialEvaluation

child.configure()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 1144, in configure

self.executeTest(self.configureLibrary)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/base.py",
 line 140, in executeTest

ret = test(*args,**kargs)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 902, in configureLibrary

for location, directory, lib, incl in self.generateGuesses():

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 476, in generateGuesses

d = self.checkDownload()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/packages/OpenMPI.py",
 line 56, in checkDownload

return self.getInstallDir()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 365, in getInstallDir

installDir = self.Install()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/packages/OpenMPI.py",
 line 63, in Install

installDir = config.package.GNUPackage.Install(self)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 1667, in Install

raise RuntimeError('Error running configure on ' + self.PACKAGE)



Finishing configure run at Sat, 10 Apr 2021 11:57:20 -0700

========

 

Thanks,

 

Danyang

 

From: Barry Smith 
Date: Saturday, April 10, 2021 at 10:31 AM
To: Danyang Su 
Cc: "petsc-users@mcs.anl.gov" 
Subject: Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI 
version

 

 

  Depending on the network you can remove the ./configure options 
--with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90  and use instead 
--with-cc=icc --with-cxx=icpc and--with-fc=ifort --download-openmpi 

 

  Barry

 

 

On Apr 10, 2021, at 12:18 PM, Danyang Su  wrote:

 

Dear PETSc developers and users,

 

I am trying to install the latest PETSc version on an ancient cluster. The 
OpenMPI version is 1.6.5 and Compiler is Intel 14.0, which are the newest on 
that cluster. I have no problem to install PETSc up to version 3.12.5. However, 
if I try to use PETSc 3.13+, there are three undefined reference errors in 
MPI_Win_allocate, MPI_Win_attach and MPI_Win_create_dynamic. I know these three 
functions are available from OpenMPI 2.0+. Because the cluster is not in 
technical support anymore, there is no way I can install new OpenMPI version or 
do some update. Is it possible to disable these three functions in PETSc 3.13+ 
version?

 

The errors occur in ‘make check’ step:

/home/dsu/soft/petsc/petsc-3.13.0/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_allocate'

/home/dsu/soft/petsc/petsc-3.13.0/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference

Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI version

Hi Barry,

 

I tried this option before but get ‘Error running configure on OpenMPI’

 

***

 UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for 
details):

---

Error running configure on OPENMPI

***

  File "/global/home/danyangs/soft/petsc/petsc-3.14.6/config/configure.py", 
line 456, in petsc_configure

framework.configure(out = sys.stdout)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1253, in configure

self.processChildren()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1242, in processChildren

self.serialEvaluation(self.childGraph)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/framework.py",
 line 1217, in serialEvaluation

child.configure()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 1144, in configure

self.executeTest(self.configureLibrary)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/base.py",
 line 140, in executeTest

ret = test(*args,**kargs)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 902, in configureLibrary

for location, directory, lib, incl in self.generateGuesses():

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 476, in generateGuesses

d = self.checkDownload()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/packages/OpenMPI.py",
 line 56, in checkDownload

return self.getInstallDir()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 365, in getInstallDir

installDir = self.Install()

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/packages/OpenMPI.py",
 line 63, in Install

installDir = config.package.GNUPackage.Install(self)

  File 
"/global/home/danyangs/soft/petsc/petsc-3.14.6/config/BuildSystem/config/package.py",
 line 1667, in Install

raise RuntimeError('Error running configure on ' + self.PACKAGE)



Finishing configure run at Sat, 10 Apr 2021 11:57:20 -0700

============

 

Thanks,

 

Danyang

 

From: Barry Smith 
Date: Saturday, April 10, 2021 at 10:31 AM
To: Danyang Su 
Cc: "petsc-users@mcs.anl.gov" 
Subject: Re: [petsc-users] Undefined reference in PETSc 3.13+ with old MPI 
version

 

 

  Depending on the network you can remove the ./configure options 
--with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90  and use instead 
--with-cc=icc --with-cxx=icpc and--with-fc=ifort --download-openmpi 

 

  Barry

 



On Apr 10, 2021, at 12:18 PM, Danyang Su  wrote:

 

Dear PETSc developers and users,

 

I am trying to install the latest PETSc version on an ancient cluster. The 
OpenMPI version is 1.6.5 and Compiler is Intel 14.0, which are the newest on 
that cluster. I have no problem to install PETSc up to version 3.12.5. However, 
if I try to use PETSc 3.13+, there are three undefined reference errors in 
MPI_Win_allocate, MPI_Win_attach and MPI_Win_create_dynamic. I know these three 
functions are available from OpenMPI 2.0+. Because the cluster is not in 
technical support anymore, there is no way I can install new OpenMPI version or 
do some update. Is it possible to disable these three functions in PETSc 3.13+ 
version?

 

The errors occur in ‘make check’ step:

/home/dsu/soft/petsc/petsc-3.13.0/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_allocate'

/home/dsu/soft/petsc/petsc-3.13.0/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_attach'

/home/dsu/soft/petsc/petsc-3.13.0/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_create_dynamic'

 

The configuration used is shown below:

./configure --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90 --download-mumps 
--download-scalapack --download-parmetis --download-metis 
--download-fblaslapack --download-hypre --download-superlu --download-hdf5=yes 
--with-debugging=0 COPTFLAGS="-O3 -march=native -mtune=native" CXXOPTFLAGS="-O3 
-march=native -mtune=native" FOPTFLAGS="-O3 -march=native -mtune=nativels"

 

Thanks,

 

Danyang

[petsc-users] Undefined reference in PETSc 3.13+ with old MPI version

Dear PETSc developers and users,

 

I am trying to install the latest PETSc version on an ancient cluster. The 
OpenMPI version is 1.6.5 and Compiler is Intel 14.0, which are the newest on 
that cluster. I have no problem to install PETSc up to version 3.12.5. However, 
if I try to use PETSc 3.13+, there are three undefined reference errors in 
MPI_Win_allocate, MPI_Win_attach and MPI_Win_create_dynamic. I know these three 
functions are available from OpenMPI 2.0+. Because the cluster is not in 
technical support anymore, there is no way I can install new OpenMPI version or 
do some update. Is it possible to disable these three functions in PETSc 3.13+ 
version?

 

The errors occur in ‘make check’ step:

/home/dsu/soft/petsc/petsc-3.13.0/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_allocate'

/home/dsu/soft/petsc/petsc-3.13.0/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_attach'

/home/dsu/soft/petsc/petsc-3.13.0/linux-intel-openmpi/lib/libpetsc.so: 
undefined reference to `MPI_Win_create_dynamic'

 

The configuration used is shown below:

 ./configure --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90 
--download-mumps --download-scalapack --download-parmetis --download-metis 
--download-fblaslapack --download-hypre --download-superlu --download-hdf5=yes 
--with-debugging=0 COPTFLAGS="-O3 -march=native -mtune=native" CXXOPTFLAGS="-O3 
-march=native -mtune=native" FOPTFLAGS="-O3 -march=native -mtune=nativels"

 

Thanks,

 

Danyang

Re: [petsc-users] Quite different behaviours of PETSc solver on different clusters

2020-10-29 Thread Danyang Su

Hi Matt,

No, interations from both linear and nonlinear solvers are similar. The system 
administrator doubt that the latency in mpich makes the difference. We will 
test a petsc version with OpenMPI on that cluster to check if it makes 
difference.

Thanks, 

Danyang

On October 29, 2020 6:05:53 p.m. PDT, Matthew Knepley  wrote:
>On Thu, Oct 29, 2020 at 3:04 PM Su,D.S. Danyang 
>wrote:
>
>> Dear PETSc users,
>>
>>
>>
>> This is a question bother me for some time. I have the same code
>running
>> on different clusters and both clusters have good speedup. However, I
>> noticed some thing quite strange. On one cluster, the solver is quite
>> stable in computing time while on another cluster, the solver is
>unstable
>> in computing time. As shown in the figure below, the local
>calculation
>> almost has no communication and the computing time in this part is
>quite
>> stable. However, PETSc solver on Cluster B jumps quite a lot and the
>> performance is not as good as Cluster A, even though the local
>calculation
>> is a little better on Cluster B. There are some difference on
>hardware and
>> PETSc configuration and optimization. Cluster A uses OpenMPI + GCC
>compiler
>> and Cluster B uses MPICH + GCC compiler. The number of processors
>used is
>> 128 on Cluster A and 120 on Cluster B. I also tested different number
>of
>> processors but the problem is the same. Does anyone have any idea
>which
>> part might cause this problem?
>>
>
>First question: Does the solver take more iterates when the time bumps
>up?
>
>  Thanks,
>
>Matt
>
>
>>
>>
>>
>>
>> Thanks and regards,
>>
>>
>>
>> Danyang
>>
>>
>>
>
>
>-- 
>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/
>

-- 
Sent from my Android device with K-9 Mail. Please excuse my brevity.

Re: [petsc-users] Parallel writing in HDF5-1.12.0 when some processors have no data to write

2020-06-14 Thread Danyang Su

Hi Jed,

I cannot reproduce the same problem in C (attached example). The C code works 
fine on HDF 1.12.0 without problem. I am still confused why this happens.  
Anyway, will use HDF5-1.10.6 with PETSc-3.13 at this moment.

Thanks,

Danyang

On 2020-06-13, 1:39 PM, "Jed Brown"  wrote:

Can you reproduce in C?  You're missing three extra arguments that exist in 
the Fortran interface.

https://support.hdfgroup.org/HDF5/doc/RM/RM_H5D.html#Dataset-Create

    Danyang Su  writes:

> Hi Jed,
>
> Attached is the example for your test.  
>
> This example uses H5Sset_none to tell H5Dwrite call that there will be no 
data. 4-th process HAS to participate we are in a collective mode.
> The code is ported and modified based on the C example from 
https://support.hdfgroup.org/ftp/HDF5/examples/parallel/coll_test.c
>
> The compiling flags in the makefile are same as those used in my own code.
>
> To compile the code, please run 'make all'
> To test the code, please run 'mpiexec -n 4 ./hdf5_zero_data'. Any number 
of processors larger than 4 should help to detect the problem.
>
> The code may crash on HDF5 1.12.0 but works fine on HDF5 1.10.x.
>
> The following platforms have been tested:
>   Macos-Mojave + GNU-8.2 + HDF5-1.12.0 -> Works fine
>   Ubuntu-16.04 + GNU-5.4 + HDF5-1.12.0 -> Crashes
>   Ubuntu-16.04 + GNU-7.5 + HDF5-1.12.0 -> Crashes
>   Ubuntu-16.04 + GNU-5.4 + HDF5-1.10.x -> Works fine
>   Centos-7 + Intel2018 + HDF5-12.0  -> Works fine
>
> Possible error when code crashes
>  At line 6686 of file H5_gen.F90
>  Fortran runtime error: Index '1' of dimension 1 of array 'buf' above 
upper bound of 0
>
> Thanks,
>
> Danyang
>
> On 2020-06-12, 6:05 AM, "Jed Brown"  wrote:
>
> Danyang Su  writes:
>
> > Hi Jed,
> >
> > Thanks for your double check. 
> >
> > The HDF 1.10.6 version also works. But versions from 1.12.x stop 
working.
>
> I'd suggest making a reduced test case in order to submit a bug 
report.
>
> This was the relevant change in PETSc for hdf5-1.12.
>
> 
https://gitlab.com/petsc/petsc/commit/806daeb7de397195b5132278177f4d5553f9f612
>
> > Attached is the code section where I have problem.
> >
> > !c write the dataset collectively
> > 
!!!
> >  CODE CRASHES HERE IF SOME PROCESSORS HAVE NO DATA TO 
WRITE
> > 
!!!
> > call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, dataset, 
hdf5_dsize,   &
> > hdf5_ierr, file_space_id=filespace, 
   &
> > mem_space_id=memspace, xfer_prp = xlist_id)
> >
> > Please let me know if there is something wrong in the code that 
causes the problem.
>
> !c
> !c This example uses H5Sset_none to tell H5Dwrite call that 
> !c there will be no data. 4-th process HAS to participate since 
> !c we are in a collective mode.
> !c 
> !c The code is ported and modified based on the C example 
> !c from 
https://support.hdfgroup.org/ftp/HDF5/examples/parallel/coll_test.c
> !c by Danyang Su on June 12, 2020.
> !c
> !c To compile the code, please run 'make all'
> !c To test the code, please run 'mpiexec -n 4 ./hdf5_zero_data'
> !c
> !c IMPORTNANT NOTE
> !c The code may crash on HDF5 1.12.0 but works fine on HDF5 1.10.x.
> !c 
> !c The following platforms have been tested:
> !c Macos-Mojave + GNU-8.2 + HDF5-1.12.0 -> Works fine
> !c Ubuntu-16.04 + GNU-5.4 + HDF5-1.12.0 -> Crashes
> !c Ubuntu-16.04 + GNU-7.5 + HDF5-1.12.0 -> Crashes
> !c Ubuntu-16.04 + GNU-5.4 + HDF5-1.10.x -> Works fine
> !c Centos-7 + Intel2018 + HDF5-12.0 -> Works fine
> !c
> !c Possible error when code crashes
> !c At line 6686 of file H5_gen.F90
> !c Fortran runtime error: Index '1' of dimension 1 of array 'buf' 
above upper bound of 0
> !c 
>
> program hdf5_zero_data
>
>
> #include 
>
>   use petscsys
>   use hdf5
>
>   implicit none 
>
>   character(len=10), parameter :: h5File_Name = "SDS_row.h5"
>   character(len=8), parameter :: DatasetName = "IntArray"
>   integer, parameter :

Re: [petsc-users] Parallel writing in HDF5-1.12.0 when some processors have no data to write

2020-06-13 Thread Danyang Su

Hi Jed,

These three extra arguments in h5dcreate_f are optional so I am not sure if 
these three arguments cause the problem. The error is from h5dwrite_f where all 
the arguments are provided.

I will try to see if I can reproduce the problem in C.

Thanks,

Danyang

On 2020-06-13, 1:39 PM, "Jed Brown"  wrote:

Can you reproduce in C?  You're missing three extra arguments that exist in 
the Fortran interface.

https://support.hdfgroup.org/HDF5/doc/RM/RM_H5D.html#Dataset-Create

    Danyang Su  writes:

> Hi Jed,
>
> Attached is the example for your test.  
>
> This example uses H5Sset_none to tell H5Dwrite call that there will be no 
data. 4-th process HAS to participate we are in a collective mode.
> The code is ported and modified based on the C example from 
https://support.hdfgroup.org/ftp/HDF5/examples/parallel/coll_test.c
>
> The compiling flags in the makefile are same as those used in my own code.
>
> To compile the code, please run 'make all'
> To test the code, please run 'mpiexec -n 4 ./hdf5_zero_data'. Any number 
of processors larger than 4 should help to detect the problem.
>
> The code may crash on HDF5 1.12.0 but works fine on HDF5 1.10.x.
>
> The following platforms have been tested:
>   Macos-Mojave + GNU-8.2 + HDF5-1.12.0 -> Works fine
>   Ubuntu-16.04 + GNU-5.4 + HDF5-1.12.0 -> Crashes
>   Ubuntu-16.04 + GNU-7.5 + HDF5-1.12.0 -> Crashes
>   Ubuntu-16.04 + GNU-5.4 + HDF5-1.10.x -> Works fine
>   Centos-7 + Intel2018 + HDF5-12.0  -> Works fine
>
> Possible error when code crashes
>  At line 6686 of file H5_gen.F90
>  Fortran runtime error: Index '1' of dimension 1 of array 'buf' above 
upper bound of 0
>
> Thanks,
>
> Danyang
>
> On 2020-06-12, 6:05 AM, "Jed Brown"  wrote:
>
> Danyang Su  writes:
>
> > Hi Jed,
> >
> > Thanks for your double check. 
> >
> > The HDF 1.10.6 version also works. But versions from 1.12.x stop 
working.
>
> I'd suggest making a reduced test case in order to submit a bug 
report.
>
> This was the relevant change in PETSc for hdf5-1.12.
>
> 
https://gitlab.com/petsc/petsc/commit/806daeb7de397195b5132278177f4d5553f9f612
>
> > Attached is the code section where I have problem.
> >
> > !c write the dataset collectively
> > 
!!!
> >  CODE CRASHES HERE IF SOME PROCESSORS HAVE NO DATA TO 
WRITE
> > 
!!!
> > call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, dataset, 
hdf5_dsize,   &
> > hdf5_ierr, file_space_id=filespace, 
   &
> > mem_space_id=memspace, xfer_prp = xlist_id)
> >
> > Please let me know if there is something wrong in the code that 
causes the problem.
>
> !c
> !c This example uses H5Sset_none to tell H5Dwrite call that 
> !c there will be no data. 4-th process HAS to participate since 
> !c we are in a collective mode.
> !c 
> !c The code is ported and modified based on the C example 
> !c from 
https://support.hdfgroup.org/ftp/HDF5/examples/parallel/coll_test.c
> !c by Danyang Su on June 12, 2020.
> !c
> !c To compile the code, please run 'make all'
> !c To test the code, please run 'mpiexec -n 4 ./hdf5_zero_data'
> !c
> !c IMPORTNANT NOTE
> !c The code may crash on HDF5 1.12.0 but works fine on HDF5 1.10.x.
> !c 
> !c The following platforms have been tested:
> !c Macos-Mojave + GNU-8.2 + HDF5-1.12.0 -> Works fine
> !c Ubuntu-16.04 + GNU-5.4 + HDF5-1.12.0 -> Crashes
> !c Ubuntu-16.04 + GNU-7.5 + HDF5-1.12.0 -> Crashes
> !c Ubuntu-16.04 + GNU-5.4 + HDF5-1.10.x -> Works fine
> !c Centos-7 + Intel2018 + HDF5-12.0 -> Works fine
> !c
> !c Possible error when code crashes
> !c At line 6686 of file H5_gen.F90
> !c Fortran runtime error: Index '1' of dimension 1 of array 'buf' 
above upper bound of 0
> !c 
>
> program hdf5_zero_data
>
>
> #include 
>
>   use petscsys
>   use hdf5
>
>   implicit none 
>
>   character(len=10), parameter :: h5File_Name = "SDS_row.h5"
>   character(len=8), parameter :: DatasetName = "IntArray"
>   integer

Re: [petsc-users] Parallel writing in HDF5-1.12.0 when some processors have no data to write

2020-06-12 Thread Danyang Su

Hi Jed,

Attached is the example for your test.  

This example uses H5Sset_none to tell H5Dwrite call that there will be no data. 
4-th process HAS to participate we are in a collective mode.
The code is ported and modified based on the C example from 
https://support.hdfgroup.org/ftp/HDF5/examples/parallel/coll_test.c

The compiling flags in the makefile are same as those used in my own code.

To compile the code, please run 'make all'
To test the code, please run 'mpiexec -n 4 ./hdf5_zero_data'. Any number of 
processors larger than 4 should help to detect the problem.

The code may crash on HDF5 1.12.0 but works fine on HDF5 1.10.x.

The following platforms have been tested:
  Macos-Mojave + GNU-8.2 + HDF5-1.12.0 -> Works fine
  Ubuntu-16.04 + GNU-5.4 + HDF5-1.12.0 -> Crashes
  Ubuntu-16.04 + GNU-7.5 + HDF5-1.12.0 -> Crashes
  Ubuntu-16.04 + GNU-5.4 + HDF5-1.10.x -> Works fine
  Centos-7 + Intel2018 + HDF5-12.0  -> Works fine

Possible error when code crashes
 At line 6686 of file H5_gen.F90
 Fortran runtime error: Index '1' of dimension 1 of array 'buf' above upper 
bound of 0

Thanks,

Danyang

On 2020-06-12, 6:05 AM, "Jed Brown"  wrote:

Danyang Su  writes:

> Hi Jed,
>
> Thanks for your double check. 
>
> The HDF 1.10.6 version also works. But versions from 1.12.x stop working.

I'd suggest making a reduced test case in order to submit a bug report.

This was the relevant change in PETSc for hdf5-1.12.

https://gitlab.com/petsc/petsc/commit/806daeb7de397195b5132278177f4d5553f9f612

> Attached is the code section where I have problem.
>
> !c write the dataset collectively
> !!!
>  CODE CRASHES HERE IF SOME PROCESSORS HAVE NO DATA TO WRITE
> !!!
> call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, dataset, hdf5_dsize,   &
> hdf5_ierr, file_space_id=filespace,&
> mem_space_id=memspace, xfer_prp = xlist_id)
>
> Please let me know if there is something wrong in the code that causes 
the problem.

hdf5_zero_data.F90
Description: Binary data

makefile
Description: Binary data

Re: [petsc-users] Parallel writing in HDF5-1.12.0 when some processors have no data to write

2020-06-12 Thread Danyang Su

Hi Jed,

Thanks for your double check. 

The HDF 1.10.6 version also works. But versions from 1.12.x stop working.

Attached is the code section where I have problem.

!c write the dataset collectively
!!!
 CODE CRASHES HERE IF SOME PROCESSORS HAVE NO DATA TO WRITE
!!!
call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, dataset, hdf5_dsize,   &
hdf5_ierr, file_space_id=filespace,&
mem_space_id=memspace, xfer_prp = xlist_id)

Please let me know if there is something wrong in the code that causes the 
problem.

Thanks,

Danyang

On 2020-06-11, 8:32 PM, "Jed Brown"  wrote:

Danyang Su  writes:

> Hi Barry,
>
> The HDF5 calls fail. I reconfigure PETSc with HDF 1.10.5 version and it 
works fine on different platforms. So, it is more likely there is a bug in the 
latest HDF version.

I would double-check that you have not subtly violated a collective 
requirement in the interface, then report to upstream.

example.F90
Description: Binary data

Re: [petsc-users] Parallel writing in HDF5-1.12.0 when some processors have no data to write

2020-06-11 Thread Danyang Su

Hi Barry,

The HDF5 calls fail. I reconfigure PETSc with HDF 1.10.5 version and it works 
fine on different platforms. So, it is more likely there is a bug in the latest 
HDF version.

Thanks.

All the best,

Danyang



On June 11, 2020 5:58:28 a.m. PDT, Barry Smith  wrote:
>
>Are you making HDF5 calls that fail or is it PETSc routines calling
>HDF5 that fail? 
>
>Regardless it sounds like the easiest fix is to switch back to the
>previous HDF5 and wait for HDF5 to fix what sounds to be a bug.
>
>   Barry
>
>
>> On Jun 11, 2020, at 1:05 AM, Danyang Su  wrote:
>> 
>> Hi All,
>>  
>> Sorry to send the previous incomplete email accidentally. 
>>  
>> After updating to HDF5-1.12.0, I got some problem if some processors
>have no data to write or not necessary to write. Since parallel writing
>is collective, I cannot disable those processors from writing. For the
>old version, there seems no such problem. So far, the problem only
>occurs on Linux using GNU compiler. The same code has no problem using
>intel compiler or latest gnu compiler on MacOS. 
>>  
>> I have already included h5sselect_none in the code for those
>processors without data. But it does not take effect. The problem is
>documented in the following link (How do you write data when one
>process doesn't have or need to write data ?).
>>  https://support.hdfgroup.org/HDF5/hdf5-quest.html#par-nodata
><https://support.hdfgroup.org/HDF5/hdf5-quest.html#par-nodata>
>>  
>> Similar problem has also been reported on HDF Forum by others.
>> https://forum.hdfgroup.org/t/bug-on-hdf5-1-12-0-fortran-parallel/6864
><https://forum.hdfgroup.org/t/bug-on-hdf5-1-12-0-fortran-parallel/6864>
>>  
>> Any suggestion for that?
>>  
>> Thanks,
>>  
>> Danyang

-- 
Sent from my Android device with K-9 Mail. Please excuse my brevity.

[petsc-users] FW: Parallel writing in HDF5-1.12.0 when some processors have no data to write

2020-06-11 Thread Danyang Su

Hi All,

 

Sorry to send the previous incomplete email accidentally. 

 

After updating to HDF5-1.12.0, I got some problem if some processors have no 
data to write or not necessary to write. Since parallel writing is collective, 
I cannot disable those processors from writing. For the old version, there 
seems no such problem. So far, the problem only occurs on Linux using GNU 
compiler. The same code has no problem using intel compiler or latest gnu 
compiler on MacOS. 

 

I have already included h5sselect_none in the code for those processors without 
data. But it does not take effect. The problem is documented in the following 
link (How do you write data when one process doesn't have or need to write data 
?).

 https://support.hdfgroup.org/HDF5/hdf5-quest.html#par-nodata

 

Similar problem has also been reported on HDF Forum by others.

https://forum.hdfgroup.org/t/bug-on-hdf5-1-12-0-fortran-parallel/6864

 

Any suggestion for that?

 

Thanks,

 

Danyang

[petsc-users] Parallel writing in HDF5-1.12.0 when some processors have no data to write

2020-06-11 Thread Danyang Su

Hi All,

 

After updating to HDF5-1.12.0, I got some problem if some processors have no 
data to write or not necessary to write. Since parallel writing is collective, 
I cannot disable those processors from writing. For the old version, there 
seems no such problem. So far, the problem only occurs on Linux using GNU 
compiler. The same code has no problem using intel compiler or latest gnu 
compiler on MacOS.

 

Looks like it is caused by zero memory space. However, as documented in

Re: [petsc-users] Bug in ex14f.F90 when debug flags are used?

2020-06-05 Thread Danyang Su

Thanks, Satish.

Danyang

On 2020-06-05, 2:11 PM, "Satish Balay"  wrote:

VecGetArray() is for F77 - it relies on  out-of-bound access.

The safer call is  VecGetArrayF90()

Now that PETSc requires F90  - perhaps VecGetArray() should be deprecated 
[and all examples fixed to use VecGetArrayF90]..

Satish

On Fri, 5 Jun 2020, Danyang Su wrote:

> Hi All,
> 
>  
> 
> I have a question regarding the following example. 
> 
> 
https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/tutorials/ex14f.F90.html
> 
>  
> 
>  
> 
> When debug flags are used in the make file, the code crashed with 
following error.
> 
>  
> 
> At line 335 of file ex14f.F90
> 
> Fortran runtime error: Index '-14450582413' of dimension 1 of array 'xx' 
below lower bound of 1
> 
>  
> 
> FFLAGS   = -g -fcheck=all -fbacktrace -Wall
> 
> CPPFLAGS = -g -fcheck=all -fbacktrace -Wall
> 
>  
> 
> Does this make sense?
> 
>  
> 
> Danyang
> 
>

[petsc-users] Bug in ex14f.F90 when debug flags are used?

2020-06-05 Thread Danyang Su

Hi All,

 

I have a question regarding the following example. 

https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/tutorials/ex14f.F90.html

 

 

When debug flags are used in the make file, the code crashed with following 
error.

 

At line 335 of file ex14f.F90

Fortran runtime error: Index '-14450582413' of dimension 1 of array 'xx' below 
lower bound of 1

 

FFLAGS   = -g -fcheck=all -fbacktrace -Wall

CPPFLAGS = -g -fcheck=all -fbacktrace -Wall

 

Does this make sense?

 

Danyang

Re: [petsc-users] Agglomeration for Multigrid on Unstructured Meshes

2020-06-01 Thread Danyang Su

Thanks Jed for the quick response. Yes I am asking about the repartitioning of 
coarse grids in geometric multigrid for unstructured mesh. I am happy with AMG. 
Thanks for letting me know.

Danyang

On 2020-06-01, 1:47 PM, "Jed Brown"  wrote:

I assume you're talking about repartitioning of coarse grids in
geometric multigrid -- that hasn't been implemented.

https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCTELESCOPE.html

But you can use an algebraic multigrid that does similar communicator
reduction, and can be applied to the original global problem or just on
the "coarse" problem of an initial geometric hierarchy.

Danyang Su  writes:

> Dear All,
>
>  
>
> I recalled there was a presentation ‘Extreme-scale multigrid components 
with PETSc’  taling about agglomeration in parallel multigrid, with future plan 
to extend to support unstructured meshes. Is this under development or to be 
added? 
>
>  
>
> Thanks and regards,
>
>  
>
> Danyang

[petsc-users] Agglomeration for Multigrid on Unstructured Meshes

2020-06-01 Thread Danyang Su

Dear All,

 

I recalled there was a presentation ‘Extreme-scale multigrid components with 
PETSc’  taling about agglomeration in parallel multigrid, with future plan to 
extend to support unstructured meshes. Is this under development or to be 
added? 

 

Thanks and regards,

 

Danyang

Re: [petsc-users] Domain decomposition using DMPLEX

2019-11-26 Thread Danyang Su


On 2019-11-26 10:18 a.m., Matthew Knepley wrote:
On Tue, Nov 26, 2019 at 11:43 AM Danyang Su <mailto:danyang...@gmail.com>> wrote:


On 2019-11-25 7:54 p.m., Matthew Knepley wrote:

On Mon, Nov 25, 2019 at 6:25 PM Swarnava Ghosh
mailto:swarnav...@gmail.com>> wrote:

Dear PETSc users and developers,

I am working with dmplex to distribute a 3D unstructured mesh
made of tetrahedrons in a cuboidal domain. I had a few queries:
1) Is there any way of ensuring load balancing based on the
number of vertices per MPI process.


You can now call DMPlexRebalanceSharedPoints() to try and get
better balance of vertices.


Hi Matt,

I just want to follow up if this new function can help to solve
the "Strange Partition in PETSc 3.11" problem I mentioned before.
Would you please let me know when shall I call this function?
Right before DMPlexDistribute?

This is not the problem. I believe the problem is that you are 
partitioning hybrid cells, and the way we handle
them internally changed, which I think screwed up the dual mesh for 
partitioning in your example. I have been

sick, so I have not gotten to your example yet, but I will.


Hope you are getting well soon. The mesh is not hybrid, only prism cells 
layer by layer. But the height of the prism varies significantly.


Thanks,

Danyang



  Sorry about that,

    Matt

call DMPlexCreateFromCellList

call DMPlexGetPartitioner

call PetscPartitionerSetFromOptions

call DMPlexDistribute

Thanks,

Danyang


2) As the global domain is cuboidal, is the resulting domain
decomposition also cuboidal on every MPI process? If not, is
there a way to ensure this? For example in DMDA, the default
domain decomposition for a cuboidal domain is cuboidal.


It sounds like you do not want something that is actually
unstructured. Rather, it seems like you want to
take a DMDA type thing and split it into tets. You can get a
cuboidal decomposition of a hex mesh easily.
Call DMPlexCreateBoxMesh() with one cell for every process,
distribute, and then uniformly refine. This
will not quite work for tets since the mesh partitioner will tend
to violate that constraint. You could:

  a) Prescribe the distribution yourself using the Shell
partitioner type

or

  b) Write a refiner that turns hexes into tets

We already have a refiner that turns tets into hexes, but we
never wrote the other direction because it was not clear
that it was useful.

  Thanks,

     Matt

Sincerely,
SG



-- 
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/>

Re: [petsc-users] Domain decomposition using DMPLEX

2019-11-26 Thread Danyang Su


On 2019-11-25 7:54 p.m., Matthew Knepley wrote:
On Mon, Nov 25, 2019 at 6:25 PM Swarnava Ghosh > wrote:


Dear PETSc users and developers,

I am working with dmplex to distribute a 3D unstructured mesh made
of tetrahedrons in a cuboidal domain. I had a few queries:
1) Is there any way of ensuring load balancing based on the number
of vertices per MPI process.


You can now call DMPlexRebalanceSharedPoints() to try and get better 
balance of vertices.


Hi Matt,

I just want to follow up if this new function can help to solve the 
"Strange Partition in PETSc 3.11" problem I mentioned before. Would you 
please let me know when shall I call this function? Right before 
DMPlexDistribute?


call DMPlexCreateFromCellList

call DMPlexGetPartitioner

call PetscPartitionerSetFromOptions

call DMPlexDistribute

Thanks,

Danyang


2) As the global domain is cuboidal, is the resulting domain
decomposition also cuboidal on every MPI process? If not, is there
a way to ensure this? For example in DMDA, the default domain
decomposition for a cuboidal domain is cuboidal.


It sounds like you do not want something that is actually 
unstructured. Rather, it seems like you want to
take a DMDA type thing and split it into tets. You can get a cuboidal 
decomposition of a hex mesh easily.
Call DMPlexCreateBoxMesh() with one cell for every process, 
distribute, and then uniformly refine. This
will not quite work for tets since the mesh partitioner will tend to 
violate that constraint. You could:


  a) Prescribe the distribution yourself using the Shell partitioner type

or

  b) Write a refiner that turns hexes into tets

We already have a refiner that turns tets into hexes, but we never 
wrote the other direction because it was not clear

that it was useful.

  Thanks,

     Matt

Sincerely,
SG



--
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/

Re: [petsc-users] DMPlex memory problem in scaling test

2019-10-10 Thread Danyang Su via petsc-users

Hi Matt,

My previous test is terminated after calling subroutine A as shown below.

>> In Subroutine A

  call DMPlexDistribute(dmda_flow%da,stencil_width,    &
    PETSC_NULL_SF,distributedMesh,ierr)
  CHKERRQ(ierr)

  if (distributedMesh /= PETSC_NULL_DM) then

    call DMDestroy(dmda_flow%da,ierr)
    CHKERRQ(ierr)
    !c set the global mesh as distributed mesh
    dmda_flow%da = distributedMesh

    call DMDestroy(distributedMesh,ierr)

   !If DMDestroy(distributedMesh,ierr) called, then everything is 
destroyed and there is nothing output with -malloc_test. However, I got 
error in the next subroutine [0]PETSC ERROR: DMGetCoordinatesLocal() 
line 5545 in /home/dsu/Soft/PETSc/petsc-dev/src/dm/interface/dm.c Object 
already free: Parameter # 1

    CHKERRQ(ierr)

 end if

>> In Subroutine B

  !c get local mesh DM and set coordinates

  call DMGetCoordinatesLocal(dmda_flow%da,gc,ierr)
  CHKERRQ(ierr)

  call DMGetCoordinateDM(dmda_flow%da,cda,ierr)
  CHKERRQ(ierr)

Thanks,

Danyang

On 2019-10-10 6:15 p.m., Matthew Knepley wrote:

On Thu, Oct 10, 2019 at 9:00 PM Danyang Su <mailto:danyang...@gmail.com>> wrote:

Labels should be destroyed with the DM. Just make a small code
that does nothing but distribute the mesh and end. If you
run with -malloc_test you should see if everythign is destroyed
properly.

  Thanks,

    Matt

Attached is the output run with -malloc_test using 2 processor.
It's a big file. How can I quick check if something is not
properly destroyed?

Everything output has not been destroyed. It looks like you did not 
destroy the distributed DM.

  Thanks,

    Matt

Thanks,

Danyang

--
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/>

Re: [petsc-users] Makefile change for PETSc3.12.0???

2019-10-02 Thread Danyang Su via petsc-users


On 2019-10-02 11:00 a.m., Balay, Satish wrote:

Can you retry with this fix:

https://gitlab.com/petsc/petsc/commit/3ae65d51d08dba2e118033664acfd64a46c9bf1d

[You can use maint branch for it]

Satish


This works. Thanks.

Danyang



On Wed, 2 Oct 2019, Danyang Su via petsc-users wrote:


Dear All,

I installed PETSc3.12.0 version and got problem in compiling my code (Fortran
and C++). The code and makefile are the same as I used for previous PETSc
version.

The error information seems the make command does not know the compiler
information. I tested this on two linux workstations and both return the same
error.

make: *** No rule to make target '../../usg/math_common.o', needed by 'exe'.
Stop.

The makefile I use is shown below:

#PETSc variables for development version, version V3.6.0 and later
include ${PETSC_DIR}/lib/petsc/conf/variables
include ${PETSC_DIR}/lib/petsc/conf/rules

CFLAGS =
CXXFLAGS = -std=c++11 -O3
CPPFLAGS = -DUSECGAL_NO
FFLAGS = -frounding-math -O3
FPPFLAGS = -DLINUX -DRELEASE -DRELEASE_X64 -DPETSC -DPETSC_HAVE_MUMPS
-DPETSC_HAVE_SUPERLU
CLEANFILES = executable-linux

SRC =./../../

OBJS = $(SRC)usg/math_common.o\
     $(SRC)usg/geometry_definition.o\
     ...
     $(SRC)updtrootdensity.o

exe: $(OBJS) chkopts
     -${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS) -o executable-linux $(OBJS)
${PETSC_LIB}

Any idea on this?

Thanks,

Danyang

Re: [petsc-users] Error running configure on SOWING

2019-09-10 Thread Danyang Su via petsc-users


Hi Barry,

With both --download-sowing-cc= and --download-sowing-cxx= specified, it 
can be configured now.


Thanks as always for all your help,

Danyang

On 2019-09-10 11:19 a.m., Smith, Barry F. wrote:

   Ahh, sorry it also needs the C++ compiler provided with 
--download-sowing-cxx= something



On Sep 10, 2019, at 1:11 PM, Danyang Su  wrote:

Sorry I forgot to attached the log file.

Attached are the log files using the following configuration:

./configure COPTFLAGS="-march=native -O2" CXXOPTFLAGS="-march=native -O2" 
FOPTFLAGS="-march=native -O2" --with-scalapack=1 
--with-scalapack-lib="[/scinet/niagara/intel/2019.1/compilers_and_libraries_2019.1.144/linux/mkl/lib/intel64/libmkl_scalapack_lp64.so,/scinet/niagara/intel/2019.1/compilers_and_libraries_2019.1.144/linux/mkl/lib/intel64/libmkl_blacs_intelmpi_lp64.so]"
 --download-parmetis=1 --download-metis=1 --download-ptscotch=1 --download-fblaslapack=1 --download-hypre=1 
--download-superlu_dist=1 --with-hdf5=1 
--with-hdf5-dir=/scinet/niagara/software/2019a/opt/intel-2019.1-intelmpi-2019.1/hdf5-mpi/1.10.4 --download-zlib=1 
--download-szlib=1 --download-ctetgen=1 --with-debugging=0 --with-cxx-dialect=C++11 
--with-mpi-dir=/scinet/niagara/intel/2019.1/compilers_and_libraries_2019.1.144/linux/mpi/intel64 
--download-sowing-cc=/scinet/niagara/intel/2019.1/compilers_and_libraries_2019.1.144/linux/bin/intel64/icc

Thanks,

Danyang

On 2019-09-10 11:03 a.m., Smith, Barry F. wrote:

   Please send the configure.log file when run with 
--download-sowing-cc=yourCcompiler  and also 
$PETSC_ARCH/externalpackages/git.sowing/config.log this will tell us why it is 
rejecting the C compiler.

    Barry



On Sep 10, 2019, at 12:43 PM, Danyang Su via petsc-users 
 wrote:

Dear All,

I am trying to install petsc-dev on a cluster with intel compiler. However, the 
configuration get stuck on SOWING.

Error running configure on SOWING: Could not execute "['./configure 
--prefix=/home/m/min3p/danyangs/soft/petsc/petsc-dev/linux-intel-opt']":
checking for ranlib... ranlib
checking for a BSD-compatible install... /usr/bin/install -c
checking whether install works... yes
checking for ar... ar
checking for gcc... no
checking for cc... no
checking for cl.exe... noconfigure: error: in 
`/gpfs/fs1/home/m/min3p/danyangs/soft/petsc/petsc-dev/linux-intel-opt/externalpackages/git.sowing':
configure: error: no acceptable C compiler found in $PATH
See `config.log' for more details

Actually the C compiler is there.

If I use GNU compiler, there is no problem. I also tried to use different 
sowing configuration as discussed on 
https://lists.mcs.anl.gov/pipermail/petsc-dev/2018-June/023070.html, but 
without success.

The configuration is

./configure COPTFLAGS="-march=native -O2" CXXOPTFLAGS="-march=native -O2" 
FOPTFLAGS="-march=native -O2" --with-scalapack=1 
--with-scalapack-lib="[/scinet/niagara/intel/2019.1/compilers_and_libraries_2019.1.144/linux/mkl/lib/intel64/libmkl_scalapack_lp64.so,/scinet/niagara/intel/2019.1/compilers_and_libraries_2019.1.144/linux/mkl/lib/intel64/libmkl_blacs_intelmpi_lp64.so]"
 --download-parmetis=1 --download-metis=1 --download-ptscotch=1 --download-fblaslapack=1 --download-hypre=1 
--download-superlu_dist=1 --with-hdf5=1 
--with-hdf5-dir=/scinet/niagara/software/2019a/opt/intel-2019.1-intelmpi-2019.1/hdf5-mpi/1.10.4 --download-zlib=1 
--download-szlib=1 --download-ctetgen=1 --with-debugging=0 --with-cxx-dialect=C++11 
--with-mpi-dir=/scinet/niagara/intel/2019.1/compilers_and_libraries_2019.1.144/linux/mpi/intel64 -download-sowing

Any suggestion on this?

Thanks and regards,

danyang

[petsc-users] Error in creating compressed data using HDF5

2019-09-02 Thread Danyang Su via petsc-users


Dear All,

Not sure if this is the right place to ask hdf5 question. I installed 
hdf5 through PETSc configuration --download-hdf5=yes. The code runs 
without problem except the function to create compressed data (red part 
shown below).


    !c create local memory space and hyperslab
    call h5screate_simple_f(hdf5_ndim, hdf5_dsize, memspace,   &
    hdf5_ierr)
    call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, &
   hdf5_offset, hdf5_count, hdf5_ierr, &
   hdf5_stride, hdf5_block)

    !c create the global file space and hyperslab
    call h5screate_simple_f(hdf5_ndim,hdf5_gdsize,filespace, &
    hdf5_ierr)
    call h5sselect_hyperslab_f(filespace, H5S_SELECT_SET_F,    &
   hdf5_goffset, hdf5_count, hdf5_ierr,    &
   hdf5_stride, hdf5_block)

    !c create a data chunking property
    call h5pcreate_f(H5P_DATASET_CREATE_F, chunk_id, hdf5_ierr)
    call h5pset_chunk_f(chunk_id, hdf5_ndim, hdf5_csize, hdf5_ierr)

    !c create compressed data, dataset must be chunked for compression
    !c the following cause crash in hdf5 library, check when new
    !c hdf5 version is available

    ! Set ZLIB / DEFLATE Compression using compression level 6.
    ! To use SZIP Compression comment out these lines.
    !call h5pset_deflate_f(chunk_id, 6, hdf5_ierr)

    ! Uncomment these lines to set SZIP Compression
    !szip_options_mask = H5_SZIP_NN_OM_F
    !szip_pixels_per_block = 16
    !call H5Pset_szip_f(chunk_id, szip_options_mask,    &
    !   szip_pixels_per_block, hdf5_ierr)

    !c create the dataset id
    call h5dcreate_f(group_id, dataname, H5T_NATIVE_INTEGER,   &
 filespace, dset_id, hdf5_ierr,    &
 dcpl_id=chunk_id)

    !c create a data transfer property
    call h5pcreate_f(H5P_DATASET_XFER_F, xlist_id, hdf5_ierr)
    call h5pset_dxpl_mpio_f(xlist_id, H5FD_MPIO_COLLECTIVE_F,  &
    hdf5_ierr)

    !c write the dataset collectively
    call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, dataset, hdf5_dsize,  &
    hdf5_ierr, file_space_id=filespace,    &
    mem_space_id=memspace, xfer_prp = xlist_id)

    call h5dclose_f(dset_id, hdf5_ierr)

    !c close resources
    call h5sclose_f(filespace, hdf5_ierr)
    call h5sclose_f(memspace, hdf5_ierr)
    call h5pclose_f(chunk_id, hdf5_ierr)
    call h5pclose_f(xlist_id, hdf5_ierr)


Both h5pset_deflate_f and H5Pset_szip_f crashes the code with error 
information as shown below. If I comment out h5pset_deflate_f and 
H5Pset_szip_f, then everything works fine.


HDF5-DIAG: Error detected in HDF5 (1.8.18) MPI-process 0:
  #000: H5D.c line 194 in H5Dcreate2(): unable to create dataset
    major: Dataset
    minor: Unable to initialize object
  #001: H5Dint.c line 455 in H5D__create_named(): unable to create and 
link to dataset

    major: Dataset
    minor: Unable to initialize object
  #002: H5L.c line 1638 in H5L_link_object(): unable to create new link 
to object

    major: Links
    minor: Unable to initialize object
  #003: H5L.c line 1882 in H5L_create_real(): can't insert link
    major: Symbol table
    minor: Unable to insert object
  #004: H5Gtraverse.c line 861 in H5G_traverse(): internal path 
traversal failed

    major: Symbol table
    minor: Object not found

Does anyone encounter this kind of error before?

Kind regards,

Danyang

Re: [petsc-users] different dof in DMDA creating

2019-08-15 Thread Danyang Su via petsc-users


Hi Barry and Matt,

Would you please give me some advice on the functions I need to use to 
set different dof to the specified nodes. For now, I use 
DMPlexCreateSection that dof is uniform throughout the domain. I am a 
bit lost in choosing the right DMPlex functions, unfortunately.


Thanks and regards,

Danyang

On 2019-02-07 1:53 p.m., Danyang Su wrote:

Thanks, Barry. DMPlex also works for my code.

Danyang

On 2019-02-07 1:14 p.m., Smith, Barry F. wrote:

   No, you would need to use the more flexible DMPlex


On Feb 7, 2019, at 3:04 PM, Danyang Su via petsc-users 
 wrote:


Dear PETSc Users,

Does DMDA support different number of degrees of freedom for 
different node? For example I have a 2D subsurface flow problem with 
the default dof = 1 throughout the domain. Now I want to add some 
sparse fractures in the domain. For the nodes connected to the 
sparse fractures, I want to set dof to 2. Is it possible to set dof 
to 2 for those nodes only?


Thanks,

Danyang

[petsc-users] Strange compiling error in DMPlexDistribute after updating PETSc to V3.11.0

2019-04-05 Thread Danyang Su via petsc-users


Hi All,

I got a strange error in calling DMPlexDistribute after updating PETSc 
to V3.11.0. There sounds no change in the interface of DMPlexDistribute 
as documented in


https://www.mcs.anl.gov/petsc/petsc-3.10/docs/manualpages/DMPLEX/DMPlexDistribute.html#DMPlexDistribute

https://www.mcs.anl.gov/petsc/petsc-3.11/docs/manualpages/DMPLEX/DMPlexDistribute.html#DMPlexDistribute

The code section is shown below.

  !c distribute mesh over processes
  call 
DMPlexDistribute(dm,stencil_width,PETSC_NULL_SF,distributed_dm,ierr)


When I use PETSc V3.10 and earlier versions, it works fine. After 
updating to latest PETSc V3.11.0, I got the following error during compiling


  call 
DMPlexDistribute(dm,stencil_width,PETSC_NULL_SF,distributed_dm,ierr)

 1
Error: Non-variable expression in variable definition context (actual 
argument to INTENT = OUT/INOUT) at (1)
/home/dsu/Soft/PETSc/petsc-3.11.0/linux-gnu-opt/lib/petsc/conf/petscrules:31: 
recipe for target '../../solver/solver_ddmethod.o' failed


The fortran example 
/home/dsu/Soft/PETSc/petsc-3.11.0/src/dm/label/examples/tutorials/ex1f90.F90 
which also uses DMPlexDistribute can be compiled without problem. Is 
there any updates in the compiler flags I need to change?


Thanks,

Danyang

Re: [petsc-users] Installation error on macOS Mojave using GNU compiler

2019-01-04 Thread Danyang Su via petsc-users




On 2019-01-03, 4:59 PM, "Balay, Satish"  wrote:

On Thu, 3 Jan 2019, Matthew Knepley via petsc-users wrote:

> On Thu, Jan 3, 2019 at 7:02 PM Danyang Su via petsc-users <
> petsc-users@mcs.anl.gov> wrote:
> 
> > Hi All,
> >
> > I am trying to install PETSc on macOS Mojave using GNU compiler.
> > First, I tried the debug version using the following configuration and 
it
> > works fine.
> >
> > ./configure --with-cc=gcc --with-cxx=g++ --with-fc=gfortran
> > --download-mpich --download-scalapack --download-parmetis 
--download-metis
> > --download-ptscotch --download-fblaslapack --download-hypre
> > --download-superlu_dist --download-hdf5=yes --download-ctetgen
> >
> > After testing debug version, I reconfigured PETSc with optimization 
turned
> > on using the following configuration. However, I got error during this 
step.
> >
> 
> Your optimization flags are not right because the compiler is producing
> AVX-512, but your linker cannot handle it. However, it looks like
> it might be that your Fortran compiler can't handle it. Do you need
> Fortran? If not, turn it off (configure with --with-fc=0) and try again.

Or use 'FOPTFLAGS=-O3' [if its indeed fortran sources causing grief]

If '-march=native -mtune=native' gives compiler/linker errors - don't use 
them.

Satish

Hi Satish,

After removing "-march=native -mtune=native", it works now.

Thanks,

Danyang

> 
>   Thanks,
> 
> Matt
> 
> 
> > ./configure --with-cc=gcc --with-cxx=g++ --with-fc=gfortran
> > --with-cxx-dialect=C++11 --download-mpich --download-scalapack
> > --download-parmetis --download-metis --download-ptscotch
> > --download-fblaslapack --download-hypre --download-superlu_dist
> > --download-hdf5=yes --download-ctetgen --with-debugging=0 COPTFLAGS="-O3
> > -march=native -mtune=native" CXXOPTFLAGS="-O3 -march=native 
-mtune=native"
> > FOPTFLAGS="-O3 -march=native -mtune=native"
> >
> > The error information is
> >
> > ctoolchain/usr/bin/ranlib: file: .libs/libmpl.a(mpl_dbg.o) has no 
symbols
> > 
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/ranlib:
> > file: .libs/libmpl.a(mpl_dbg.o) has no symbols
> > 
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/ranlib:
> > file: .libs/libmpl.a(mpl_dbg.o) has no symbols
> > /var/folders/jm/wcm4mv8s3v1gqz383tcf_4c0gp/T//ccrlfFuo.s:14:2: 
error:
> > instruction requires: AVX-512 ISA AVX-512 VL ISA
> > vmovdqu64   (%rdi), %xmm0
> > ^
> > make[2]: *** [src/binding/fortran/use_mpi/mpi_constants.mod-stamp] 
Error 1
> > make[2]: *** Waiting for unfinished jobs
> > make[1]: *** [all-recursive] Error 1
> > make: *** [all] Error 2
> >
> >
> > Thanks,
> >
> > Danyang
> >
> 
> 
>

Re: [petsc-users] [petsc-maint] Fwd: DMPlex global to natural problem using DmPlexGetVertexNumbering or DMPlexGlobalToNatural

2018-11-29 Thread Danyang Su via petsc-users



On 18-11-29 06:13 PM, Matthew Knepley wrote:
On Thu, Nov 29, 2018 at 7:40 PM Danyang Su via petsc-maint 
mailto:petsc-ma...@mcs.anl.gov>> wrote:


Dear PETSc developers & users,

Sorry to bother you again. I just encounter some difficulties in
DMPlex
global to natural order. This is not the very necessary function
in my
code. But it's best to have it in case someone wants to feed code
with
initial conditions or parameters from external file using natural
ordering.

First I tried to use DmPlexGetVertexNumbering which seems pretty
straightforward, but I always get error saying "You need a ISO C
conforming compiler to use the glibc headers". I use Gfortran on
Linux.

Then I switched to use Label to save global-natural order before
distributing the mesh, everything works fine. The only problem is
that
it takes a very long time during mesh distribution when the size
of mesh
is very big.

Finally I tried to use GlobalToNatural ordering. Set Natural order to
TRUE before calling DMPlexDistribute and then call
DMPlexSetMigrationSF.
However, the function DMPlexGlobalToNaturalEnd seems doing nothing as
the returned vector is always unchanged. I got some help from Josh
who
has done similar work before. But I still cannot figure out what is
wrong in the code.


Can you run the example that Blaise recommended? Sometimes its easiest 
to start

with an example and change it into your code.

I cannot compile ex26 due to missing exodusII.h file. Error information is
make ex26
/home/dsu/Soft/PETSc/petsc-3.10.2/linux-gnu-dbg/bin/mpicc -o ex26.o -c 
-Wall -Wwrite-strings -Wno-strict-aliasing -Wno-unknown-pragmas 
-fstack-protector -fvisibility=hidden -g3 
-I/home/dsu/Soft/PETSc/petsc-3.10.2/include 
-I/home/dsu/Soft/PETSc/petsc-3.10.2/linux-gnu-dbg/include `pwd`/ex26.c
/home/dsu/Soft/PETSc/petsc-3.10.2/src/dm/impls/plex/examples/tests/ex26.c:10:22: 
fatal error: exodusII.h: No such file or directory

compilation terminated.
/home/dsu/Soft/PETSc/petsc-3.10.2/lib/petsc/conf/rules:359: recipe for 
target 'ex26.o' failed


Anyway, I will look at this example to get familiar with GlobalToNatural 
related functions.


Thanks,

Danyang



  Thanks,

Matt

Attached is the code section with related functions included,
together
with a mesh and output of global vector and natural vector. If
something
is wrong in the code, it should be in these four functions. The
related
functions are called in the following order.

call solver_dd_create_dmplex

call solver_dd_mapping_set_dmplex

call solver_dd_DMDACreate_flow

!call solver_dd_DMDACreate_reactNot used in current testing

call solver_dd_mapping_global_natural

I really appreciate your help.

Regards,

Danyang



--
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/%7Eknepley/>

Re: [petsc-users] Makefile for mixed C++ and Fortran code

2018-06-03 Thread Danyang Su


Hi Satish,

Your makefile works.

Please ignore the CGAL dependency. I just modified from CGAL example and 
it needs to be further optimized with the dependency that is required 
for the algorithm I use.


Thanks,

Danyang


On 18-06-02 02:51 PM, Satish Balay wrote:

Try the attached makefile. [with correct PETSC_DIR and PETSC_ARCH values]

If you have issues - send the complete makefiles - and complete error log..


On Sat, 2 Jun 2018, Danyang Su wrote:


Hi Barry,

For the code without PETSc, the rules used to compile the code with CGAL is

note: DLIB  can probably be simlified - and a fewof the options eliminated.



DLIB = -lstdc++ -lmetis -lm -L/usr/local/lib -rdynamic

-lstdc++ is setup by petsc

-lmetis can be a dependency of petsc - so care should be taken to have only one 
copy of metis


/usr/local/lib/libmpfr.so /usr/local/lib/libgmp.so

Why are these here? Normally these are dependencies of gcc/gfortran [and PETSc 
configure picks up the correct ones]


/usr/local/lib/libCGAL_ImageIO.so.11.0.1 /usr/local/lib/libCGAL.so.11.0.1

I have no idea what these are..

Satish


/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so -lpthread
-lGLU -lGL -lX11 -lXext -lz /usr/local/lib/libCGAL_ImageIO.so.11.0.1
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext -lz
/usr/local/lib/libmpfr.so /usr/local/lib/libgmp.so
/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so -lpthread
-Wl,-rpath,/usr/local/lib

FFLAGS = -O3 -I$(LIS_INC)
CXXFLAGS = -std=c++11 -O3 -I$(LIS_INC)

However, after adding these to the makefile using PETSc, I got error telling
me that all the petsc include files cannot be found.

../../solver/solver_snes_common.F90:27:0: fatal error: petscversion.h: No such
\
file or directory
  #include 

../../solver/solver_snes_common.F90:27:0: fatal error: petscversion.h: No such
\
file or directory
  #include 

Similar for other head files. However, If I change the file path to the the
full path, the code cannot compiled. Does the rule I use destroy the PETSc
relative path/

The make commands I use is

executable: $(SOURCES) chkopts
 -${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out
$(SOURCES) ${PETSC_LIB} ${LIS_LIB} ${DLIB}
%.o:%.F90
 $(FLINKER) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
 $(CLINKER) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@

Thanks,

Danyang

On 18-06-01 10:41 AM, Smith, Barry F. wrote:

 You need to determine exactly what flags are passed to the C++ compiler
 for your compile that works and make sure those same flags are used in
 "PETSc version" of the makefile. You could add the flags directly to the
 rule


%.o:%.cpp
 $(CLINKER) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@

Barry



On Jun 1, 2018, at 12:37 PM, Danyang Su  wrote:

Follow up:

With following command

executable: $(SOURCES) chkopts
 -${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out
 $(SOURCES) ${PETSC_LIB}

%.o:%.F90
 $(FLINKER) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
 $(CLINKER) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@

The compiler return error: no match function.

../../usg/cgal_triangulation_2d.cpp: In function ‘void
outputTriangulation2d(in\
t, const char*, int, const char*)’:
../../usg/cgal_triangulation_2d.cpp:485:20: error: no matching function for
cal\
l to ‘std::basic_ofstream::open(std::string&)’
 out.open(strfile);

Thanks,


Danyang

On 18-06-01 10:07 AM, Danyang Su wrote:

Hi All,

My code needs to link to an external C++ library (CGAL). The code is
written in Fortran and I have already written interface to let Fortran
call C++ function. For the sequential version without PETSc, it can be
compiled without problem using the following makefile. The parallel
version without CGAL can also be compiled successfully. However, when I
tried to use PETSc together with CGAL library, I cannot compile the code.
My questions is: How can I modify the makefile? Do I need to reconfigure
PETSc with special flags? All the makefile samples are shown below.

#makefile for sequential version

FC = gfortran
#FC = ifort
CXX = g++ -std=c++11

DLIB = -lstdc++ -lm -L/usr/local/lib -rdynamic /usr/local/lib/libmpfr.so
/usr/local/lib/libgmp.so /usr/local/lib/libCGAL_ImageIO.so.11.0.1
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext -lz
/usr/local/lib/libCGAL_ImageIO.so.11.0.1 /usr/local/lib/libCGAL.so.11.0.1
/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so
-lpthread -lGLU -lGL -lX11 -lXext -lz /usr/local/lib/libmpfr.so
/usr/local/lib/libgmp.so /usr/local/lib/libboost_thread.so
/usr/local/lib/libboost_system.so -lpthread -Wl,-rpath,/usr/local/lib

FFLAGS = -O3
CXXFLAGS = -O3

FPPFLAGS =  -DUSECGAL

SRC =./../../

SOURCES = $(SRC)usg/

Re: [petsc-users] Makefile for mixed C++ and Fortran code

2018-06-02 Thread Danyang Su


Hi Barry,

For the code without PETSc, the rules used to compile the code with CGAL is

DLIB = -lstdc++ -lmetis -lm -L/usr/local/lib -rdynamic 
/usr/local/lib/libmpfr.so /usr/local/lib/libgmp.so 
/usr/local/lib/libCGAL_ImageIO.so.11.0.1 
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so 
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext -lz 
/usr/local/lib/libCGAL_ImageIO.so.11.0.1 
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so 
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext -lz 
/usr/local/lib/libmpfr.so /usr/local/lib/libgmp.so 
/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so 
-lpthread -Wl,-rpath,/usr/local/lib


FFLAGS = -O3 -I$(LIS_INC)
CXXFLAGS = -std=c++11 -O3 -I$(LIS_INC)

However, after adding these to the makefile using PETSc, I got error 
telling me that all the petsc include files cannot be found.


../../solver/solver_snes_common.F90:27:0: fatal error: petscversion.h: 
No such \

file or directory
 #include 

../../solver/solver_snes_common.F90:27:0: fatal error: petscversion.h: 
No such \

file or directory
 #include 

Similar for other head files. However, If I change the file path to the 
the full path, the code cannot compiled. Does the rule I use destroy the 
PETSc relative path/


The make commands I use is

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out 
$(SOURCES) ${PETSC_LIB} ${LIS_LIB} ${DLIB}

%.o:%.F90
$(FLINKER) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CLINKER) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@

Thanks,

Danyang

On 18-06-01 10:41 AM, Smith, Barry F. wrote:

You need to determine exactly what flags are passed to the C++ compiler for your 
compile that works and make sure those same flags are used in "PETSc version" 
of the makefile. You could add the flags directly to the rule


%.o:%.cpp
$(CLINKER) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@

   Barry



On Jun 1, 2018, at 12:37 PM, Danyang Su  wrote:

Follow up:

With following command

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out $(SOURCES) 
${PETSC_LIB}

%.o:%.F90
$(FLINKER) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CLINKER) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@

The compiler return error: no match function.

../../usg/cgal_triangulation_2d.cpp: In function ‘void outputTriangulation2d(in\
t, const char*, int, const char*)’:
../../usg/cgal_triangulation_2d.cpp:485:20: error: no matching function for cal\
l to ‘std::basic_ofstream::open(std::string&)’
out.open(strfile);

Thanks,


Danyang

On 18-06-01 10:07 AM, Danyang Su wrote:

Hi All,

My code needs to link to an external C++ library (CGAL). The code is written in 
Fortran and I have already written interface to let Fortran call C++ function. 
For the sequential version without PETSc, it can be compiled without problem 
using the following makefile. The parallel version without CGAL can also be 
compiled successfully. However, when I tried to use PETSc together with CGAL 
library, I cannot compile the code. My questions is: How can I modify the 
makefile? Do I need to reconfigure PETSc with special flags? All the makefile 
samples are shown below.

#makefile for sequential version

FC = gfortran
#FC = ifort
CXX = g++ -std=c++11

DLIB = -lstdc++ -lm -L/usr/local/lib -rdynamic /usr/local/lib/libmpfr.so 
/usr/local/lib/libgmp.so /usr/local/lib/libCGAL_ImageIO.so.11.0.1 
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so 
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext -lz 
/usr/local/lib/libCGAL_ImageIO.so.11.0.1 /usr/local/lib/libCGAL.so.11.0.1 
/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so -lpthread 
-lGLU -lGL -lX11 -lXext -lz /usr/local/lib/libmpfr.so /usr/local/lib/libgmp.so 
/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so -lpthread 
-Wl,-rpath,/usr/local/lib

FFLAGS = -O3
CXXFLAGS = -O3

FPPFLAGS =  -DUSECGAL

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES)
$(FC) $(FFLAGS) $(FPPFLAGS) -o executable.out $(SOURCES) ${LIS_LIB} $(DLIB)
%.o:%.F90
$(FC) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CXX) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@


#makefile for parallel version with PETSc, without CGAL

#FC = ifort
#FC = gfortran

DLIB = -lm

FFLAGS = -O3

FPPFLAGS =  -DUSEPETSC

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out $(SOURCES) 
${PETSC_LIB}


#makefile for parallel version with PETSc, with CGAL, CANNOT work

#FC = ifort
#FC = gfortran

Re: [petsc-users] Makefile for mixed C++ and Fortran code

2018-06-01 Thread Danyang Su


On 18-06-01 10:29 AM, Smith, Barry F. wrote:

   What happens if you add ${DLIB} to the end of the line


executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out $(SOURCES) 
${PETSC_LIB}

   Send all the output from trying this.

Barry

Thanks, Barry. Still not work after including ${DLIB}






On Jun 1, 2018, at 12:07 PM, Danyang Su  wrote:

Hi All,

My code needs to link to an external C++ library (CGAL). The code is written in 
Fortran and I have already written interface to let Fortran call C++ function. 
For the sequential version without PETSc, it can be compiled without problem 
using the following makefile. The parallel version without CGAL can also be 
compiled successfully. However, when I tried to use PETSc together with CGAL 
library, I cannot compile the code. My questions is: How can I modify the 
makefile? Do I need to reconfigure PETSc with special flags? All the makefile 
samples are shown below.

#makefile for sequential version

FC = gfortran
#FC = ifort
CXX = g++ -std=c++11

DLIB = -lstdc++ -lm -L/usr/local/lib -rdynamic /usr/local/lib/libmpfr.so 
/usr/local/lib/libgmp.so /usr/local/lib/libCGAL_ImageIO.so.11.0.1 
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so 
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext -lz 
/usr/local/lib/libCGAL_ImageIO.so.11.0.1 /usr/local/lib/libCGAL.so.11.0.1 
/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so -lpthread 
-lGLU -lGL -lX11 -lXext -lz /usr/local/lib/libmpfr.so /usr/local/lib/libgmp.so 
/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so -lpthread 
-Wl,-rpath,/usr/local/lib

FFLAGS = -O3
CXXFLAGS = -O3

FPPFLAGS =  -DUSECGAL

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES)
$(FC) $(FFLAGS) $(FPPFLAGS) -o executable.out $(SOURCES) ${LIS_LIB} $(DLIB)
%.o:%.F90
$(FC) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CXX) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@


#makefile for parallel version with PETSc, without CGAL

#FC = ifort
#FC = gfortran

DLIB = -lm

FFLAGS = -O3

FPPFLAGS =  -DUSEPETSC

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out $(SOURCES) 
${PETSC_LIB}


#makefile for parallel version with PETSc, with CGAL, CANNOT work

#FC = ifort
#FC = gfortran

DLIB = -lm

FFLAGS = -O3

FPPFLAGS =  -DUSEPETSC -DUSECGAL

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out $(SOURCES) 
${PETSC_LIB}

%.o:%.F90
$(FC) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CXX) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@


Thanks,

Danyang

Re: [petsc-users] Makefile for mixed C++ and Fortran code

2018-06-01 Thread Danyang Su


Follow up:

With following command

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out 
$(SOURCES) ${PETSC_LIB}


%.o:%.F90
$(FLINKER) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CLINKER) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@

The compiler return error: no match function.

../../usg/cgal_triangulation_2d.cpp: In function ‘void 
outputTriangulation2d(in\

t, const char*, int, const char*)’:
../../usg/cgal_triangulation_2d.cpp:485:20: error: no matching function 
for cal\

l to ‘std::basic_ofstream::open(std::string&)’
out.open(strfile);

Thanks,


Danyang

On 18-06-01 10:07 AM, Danyang Su wrote:

Hi All,

My code needs to link to an external C++ library (CGAL). The code is 
written in Fortran and I have already written interface to let Fortran 
call C++ function. For the sequential version without PETSc, it can be 
compiled without problem using the following makefile. The parallel 
version without CGAL can also be compiled successfully. However, when 
I tried to use PETSc together with CGAL library, I cannot compile the 
code. My questions is: How can I modify the makefile? Do I need to 
reconfigure PETSc with special flags? All the makefile samples are 
shown below.


#makefile for sequential version

FC = gfortran
#FC = ifort
CXX = g++ -std=c++11

DLIB = -lstdc++ -lm -L/usr/local/lib -rdynamic 
/usr/local/lib/libmpfr.so /usr/local/lib/libgmp.so 
/usr/local/lib/libCGAL_ImageIO.so.11.0.1 
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so 
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext 
-lz /usr/local/lib/libCGAL_ImageIO.so.11.0.1 
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so 
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext 
-lz /usr/local/lib/libmpfr.so /usr/local/lib/libgmp.so 
/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so 
-lpthread -Wl,-rpath,/usr/local/lib


FFLAGS = -O3
CXXFLAGS = -O3

FPPFLAGS =  -DUSECGAL

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES)
$(FC) $(FFLAGS) $(FPPFLAGS) -o executable.out $(SOURCES) 
${LIS_LIB} $(DLIB)

%.o:%.F90
$(FC) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CXX) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@


#makefile for parallel version with PETSc, without CGAL

#FC = ifort
#FC = gfortran

DLIB = -lm

FFLAGS = -O3

FPPFLAGS =  -DUSEPETSC

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out 
$(SOURCES) ${PETSC_LIB}



#makefile for parallel version with PETSc, with CGAL, CANNOT work

#FC = ifort
#FC = gfortran

DLIB = -lm

FFLAGS = -O3

FPPFLAGS =  -DUSEPETSC -DUSECGAL

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out 
$(SOURCES) ${PETSC_LIB}


%.o:%.F90
$(FC) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CXX) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@


Thanks,

Danyang

[petsc-users] Makefile for mixed C++ and Fortran code

2018-06-01 Thread Danyang Su


Hi All,

My code needs to link to an external C++ library (CGAL). The code is 
written in Fortran and I have already written interface to let Fortran 
call C++ function. For the sequential version without PETSc, it can be 
compiled without problem using the following makefile. The parallel 
version without CGAL can also be compiled successfully. However, when I 
tried to use PETSc together with CGAL library, I cannot compile the 
code. My questions is: How can I modify the makefile? Do I need to 
reconfigure PETSc with special flags? All the makefile samples are shown 
below.


#makefile for sequential version

FC = gfortran
#FC = ifort
CXX = g++ -std=c++11

DLIB = -lstdc++ -lm -L/usr/local/lib -rdynamic /usr/local/lib/libmpfr.so 
/usr/local/lib/libgmp.so /usr/local/lib/libCGAL_ImageIO.so.11.0.1 
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so 
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext -lz 
/usr/local/lib/libCGAL_ImageIO.so.11.0.1 
/usr/local/lib/libCGAL.so.11.0.1 /usr/local/lib/libboost_thread.so 
/usr/local/lib/libboost_system.so -lpthread -lGLU -lGL -lX11 -lXext -lz 
/usr/local/lib/libmpfr.so /usr/local/lib/libgmp.so 
/usr/local/lib/libboost_thread.so /usr/local/lib/libboost_system.so 
-lpthread -Wl,-rpath,/usr/local/lib


FFLAGS = -O3
CXXFLAGS = -O3

FPPFLAGS =  -DUSECGAL

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES)
$(FC) $(FFLAGS) $(FPPFLAGS) -o executable.out $(SOURCES) ${LIS_LIB} 
$(DLIB)

%.o:%.F90
$(FC) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CXX) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@


#makefile for parallel version with PETSc, without CGAL

#FC = ifort
#FC = gfortran

DLIB = -lm

FFLAGS = -O3

FPPFLAGS =  -DUSEPETSC

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out 
$(SOURCES) ${PETSC_LIB}



#makefile for parallel version with PETSc, with CGAL, CANNOT work

#FC = ifort
#FC = gfortran

DLIB = -lm

FFLAGS = -O3

FPPFLAGS =  -DUSEPETSC -DUSECGAL

SRC =./../../

SOURCES = $(SRC)usg/math_common.o\
$(SRC)usg/geometry_definition.o\
$(SRC)usg/cgal_common.o\

...

executable: $(SOURCES) chkopts
-${FLINKER} $(FFLAGS) $(FPPFLAGS) $(CPPFLAGS)  -o executable.out 
$(SOURCES) ${PETSC_LIB}


%.o:%.F90
$(FC) $(FFLAGS) $(FPPFLAGS) -c -frounding-math $< -o $@
%.o:%.cpp
$(CXX) $(CXXFLAGS) $(CPPFLAGS) -c -frounding-math $< -o $@


Thanks,

Danyang

Re: [petsc-users] Segmentation Violation in getting DMPlex coordinates

2018-04-28 Thread Danyang Su


Hi Matt and Barry,

Thanks for your quick response. After changing DMDAVecGetArrayF90 to 
VecGetArrayF90, everything works now.


Thanks,

Danyang


On 18-04-28 01:38 PM, Smith, Barry F. wrote:

   Added runtime error checking for such incorrect calls in 
barry/dmda-calls-type-check



On Apr 28, 2018, at 9:19 AM, Matthew Knepley <knep...@gmail.com> wrote:

On Sat, Apr 28, 2018 at 2:08 AM, Danyang Su <danyang...@gmail.com> wrote:
Hi All,

I use DMPlex and need to get coordinates back after distribution. However, I 
always get segmentation violation in getting coords values in the following 
codes if using multiple processors. If only one processor is used, it works 
fine.

For each processors, the off value starts from 0 which looks good. I also tried 
0-based index, which gives the same error. Would any one help to check what is 
wrong here?

  idof   1 off   0
  idof   2 off   0
  idof   1 off   2
  idof   2 off   2
  idof   1 off   4
  idof   2 off   4
  idof   1 off   6
  idof   2 off   6
  idof   1 off   8
  idof   2 off   8


   DM :: distributedMesh, cda
   Vec :: gc
   PetscScalar, pointer :: coords(:)
   PetscSection ::  cs

   ...

   call DMGetCoordinatesLocal(dmda_flow%da,gc,ierr)
   CHKERRQ(ierr)

   call DMGetCoordinateDM(dmda_flow%da,cda,ierr)
   CHKERRQ(ierr)

   call DMGetDefaultSection(cda,cs,ierr)
   CHKERRQ(ierr)

   call PetscSectionGetChart(cs,istart,iend,ierr)
   CHKERRQ(ierr)

   !c get coordinates array
   call DMDAVecGetArrayF90(cda,gc,coords,ierr)

You cannot call DMDA function if you have a DMPlex. You jsut call 
VecGetArrayF90()

Matt
  
   CHKERRQ(ierr)


   do ipoint = istart, iend-1

 call PetscSectionGetDof(cs,ipoint,dof,ierr)
 CHKERRQ(ierr)

 call PetscSectionGetOffset(cs,ipoint,off,ierr)
 CHKERRQ(ierr)

 inode = ipoint-istart+1

 if (cell_coords == coords_xyz) then
   nodes(inode)%x = coords(off+1)
   nodes(inode)%y = coords(off+2)
   nodes(inode)%z = coords(off+3)
 else if (cell_coords == coords_xy) then
   nodes(inode)%x = coords(off+1)
   nodes(inode)%y = coords(off+2)
   nodes(inode)%z = 0.0d0
 else if (cell_coords == coords_yz) then
   nodes(inode)%x = 0.0d0
   nodes(inode)%y = coords(off+1)
   nodes(inode)%z = coords(off+2)
 else if (cell_coords ==coords_xz) then
   nodes(inode)%x = coords(off+1)
   nodes(inode)%y = 0.0d0
   nodes(inode)%z = coords(off+2)
 end if
   end do

   call DMDAVecRestoreArrayF90(cda,gc,coords,ierr)
   CHKERRQ(ierr)

Thanks,

Danyang





--
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/

[petsc-users] Segmentation Violation in getting DMPlex coordinates

2018-04-28 Thread Danyang Su


Hi All,

I use DMPlex and need to get coordinates back after distribution. 
However, I always get segmentation violation in getting coords values in 
the following codes if using multiple processors. If only one processor 
is used, it works fine.


For each processors, the off value starts from 0 which looks good. I 
also tried 0-based index, which gives the same error. Would any one help 
to check what is wrong here?


 idof   1 off   0
 idof   2 off   0
 idof   1 off   2
 idof   2 off   2
 idof   1 off   4
 idof   2 off   4
 idof   1 off   6
 idof   2 off   6
 idof   1 off   8
 idof   2 off   8


  DM :: distributedMesh, cda
  Vec :: gc
  PetscScalar, pointer :: coords(:)
  PetscSection ::  cs

  ...

  call DMGetCoordinatesLocal(dmda_flow%da,gc,ierr)
  CHKERRQ(ierr)

  call DMGetCoordinateDM(dmda_flow%da,cda,ierr)
  CHKERRQ(ierr)

  call DMGetDefaultSection(cda,cs,ierr)
  CHKERRQ(ierr)

  call PetscSectionGetChart(cs,istart,iend,ierr)
  CHKERRQ(ierr)

  !c get coordinates array
  call DMDAVecGetArrayF90(cda,gc,coords,ierr)
  CHKERRQ(ierr)

  do ipoint = istart, iend-1

call PetscSectionGetDof(cs,ipoint,dof,ierr)
CHKERRQ(ierr)

call PetscSectionGetOffset(cs,ipoint,off,ierr)
CHKERRQ(ierr)

inode = ipoint-istart+1

if (cell_coords == coords_xyz) then
  nodes(inode)%x = coords(off+1)
  nodes(inode)%y = coords(off+2)
  nodes(inode)%z = coords(off+3)
else if (cell_coords == coords_xy) then
  nodes(inode)%x = coords(off+1)
  nodes(inode)%y = coords(off+2)
  nodes(inode)%z = 0.0d0
else if (cell_coords == coords_yz) then
  nodes(inode)%x = 0.0d0
  nodes(inode)%y = coords(off+1)
  nodes(inode)%z = coords(off+2)
else if (cell_coords ==coords_xz) then
  nodes(inode)%x = coords(off+1)
  nodes(inode)%y = 0.0d0
  nodes(inode)%z = coords(off+2)
end if
  end do

  call DMDAVecRestoreArrayF90(cda,gc,coords,ierr)
  CHKERRQ(ierr)

Thanks,

Danyang

Re: [petsc-users] Get vertex index of each cell in DMPlex after distribution

2018-04-27 Thread Danyang Su

On 2018-04-27 04:11 AM, Matthew Knepley wrote:
On Fri, Apr 27, 2018 at 2:09 AM, Danyang Su <danyang...@gmail.com
<mailto:danyang...@gmail.com>> wrote:

Hi Matt,

Sorry if this is a stupid question.

In the previous code for unstructured grid, I create labels to
mark the original node/cell index from VTK file and then
distribute it so that each subdomain has a copy of its original
node and cell index, as well as the PETSc numbering. Now I am
trying to get avoid of using large number of keys in
DMSetLabelValue since this costs lot of time for large problem.

I can get the coordinates of subdomain after distribution by using
DMGetCoordinatesLocal and DMGetCoordinateDM.

How can I get the vertex index of each cell after distribution?
Would you please give me a hint or functions that I can use.

You can permute the vectors back to the natural ordering using

http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMPLEX/DMPlexNaturalToGlobalBegin.html

which says you have to call DMPlexSetUseNaturalSF() before
distributing the mesh. It is tested in

src/dm/impls/plex/examples/tests/ex15.c

so you can see how its intended to work. It is very new and has not
been tested by many people.

I can see how you might want this for small tests. Why would you want
it for production models?

Hi Matt,

This is indeed what I need. As some of years old cases import initial
conditions from external files, which are in natural ordering as the
original mesh. Just want to make the code compatible to the old input files.

Thanks,

Danyang

Thanks,

Matt

Thanks,

Danyang

On 18-04-25 02:12 PM, Danyang Su wrote:

On 2018-04-25 09:47 AM, Matthew Knepley wrote:

On Wed, Apr 25, 2018 at 12:40 PM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:

Hi Matthew,

In the worst case, every node/cell may have different label.

Do not use Label for this. Its not an appropriate thing. If
every cell is different, just use the cell number.
Labels are for mapping a relatively small number of keys (like
material IDs) to sets of points (cells, vertices, etc.)
Its not a great data structure for a permutation.

Yes. If there is small number of keys, it runs very fast, even
for more than one million DMSetLabelValue calls. The performance
just deteriorates as the number of keys increases.

I cannot get avoid of DMSetLabelValue as node/cell index of
original mesh is needed for the previous input file that uses
some of global node/cell index to set value. But if I can get the
natural order of nodes/cells from DMPlex, I can discard the use
of DMSetLabelValue. Is there any function can do this job?

Thanks,

Danyang

However, I still do not believe these numbers. The old code does
a string comparison every time. I will setup a test.

Matt

Below is one of the worst scenario with 102299 nodes and
102299 different labels for test. I found the time cost
increase during the loop. The first 9300 loop takes least
time (<0.5) while the last 9300 loops takes much more time
(>7.7), as shown below. If I use larger mesh with >1 million
nodes, it runs very very slowly in this part. The PETSc is
configured with optimization on.

Configure options --with-cc=gcc --with-cxx=g++
--with-fc=gfortran --download-mpich --download-scalapack
--download-parmetis --download-metis --download-ptscotch
--download-fblaslapack --download-hypre
--download-superlu_dist --download-hdf5=yes
--download-ctetgen --with-debugging=0 COPTFLAGS="-O3
-march=native -mtune=native" CXXOPTFLAGS="-O3 -march=native
-mtune=native" FOPTFLAGS="-O3 -march=native -mtune=native"

istart iendprogressCPU_Timetime cost - old (sec)
time cost - new (sec)
0 92990 1524670045.51166

930018599 0.100010753 1524670045.996050.4843890667
0.497246027
18600 27899 0.200010747 1524670047.326351.330302
1.3820912838
27900 37199 0.300010741 1524670049.3066 1.9802515507
2.2439446449
37200 46499 0.400010765 1524670052.1594 2.852804184
3.0739262104
46500 55799 0.500010729 1524670055.909613.7502081394
3.9270553589
55800 65099 0.600010753 1524670060.476544.5669286251
4.7571902275
65100 74399 0.700010777 1524670066.0941 5.6175630093
5.7428796291
74400 83699 0.800010741 1524670072.538866.44475317
6.5761549473
83700 92998 0.900010765 1524670079.990727.4518604279
7.4606924057

[petsc-users] Get vertex index of each cell in DMPlex after distribution

2018-04-27 Thread Danyang Su


Hi Matt,

Sorry if this is a stupid question.

In the previous code for unstructured grid, I create labels to mark the 
original node/cell index from VTK file and then distribute it so that 
each subdomain has a copy of its original node and cell index, as well 
as the PETSc numbering. Now I am trying to get avoid of using large 
number of keys in DMSetLabelValue since this costs lot of time for large 
problem.


I can get the coordinates of subdomain after distribution by using 
DMGetCoordinatesLocal and DMGetCoordinateDM.


How can I get the vertex index of each cell after distribution? Would 
you please give me a hint or functions that I can use.


Thanks,

Danyang


On 18-04-25 02:12 PM, Danyang Su wrote:

On 2018-04-25 09:47 AM, Matthew Knepley wrote:
On Wed, Apr 25, 2018 at 12:40 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi Matthew,

In the worst case, every node/cell may have different label.

Do not use Label for this. Its not an appropriate thing. If every 
cell is different, just use the cell number.
Labels are for mapping a relatively small number of keys (like 
material IDs) to sets of points (cells, vertices, etc.)

Its not a great data structure for a permutation.
Yes. If there is small number of keys, it runs very fast, even for 
more than one million DMSetLabelValue calls. The performance just 
deteriorates as the number of keys increases.


I cannot get avoid of DMSetLabelValue as node/cell index of original 
mesh is needed for the previous input file that uses some of global 
node/cell index to set value. But if I can get the natural order of 
nodes/cells from DMPlex, I can discard the use of DMSetLabelValue. Is 
there any function can do this job?


Thanks,

Danyang


However, I still do not believe these numbers. The old code does a 
string comparison every time. I will setup a test.


   Matt

Below is one of the worst scenario with 102299 nodes and 102299
different labels for test. I found the time cost increase during
the loop. The first 9300 loop takes least time (<0.5) while the
last 9300 loops takes much more time (>7.7), as shown below. If I
use larger mesh with >1 million nodes, it runs very very slowly
in this part. The PETSc is configured with optimization on.

Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran
--download-mpich --download-scalapack --download-parmetis
--download-metis --download-ptscotch --download-fblaslapack
--download-hypre --download-superlu_dist --download-hdf5=yes
--download-ctetgen --with-debugging=0 COPTFLAGS="-O3
-march=native -mtune=native" CXXOPTFLAGS="-O3 -march=native
-mtune=native" FOPTFLAGS="-O3 -march=native -mtune=native"

istart  iendprogressCPU_Timetime cost - old (sec)   
time
cost - new (sec)
0   92990   1524670045.51166

930018599   0.100010753 1524670045.996050.4843890667
0.497246027
18600   27899   0.200010747 1524670047.326351.330302
1.3820912838
27900   37199   0.300010741 1524670049.3066 1.9802515507
2.2439446449
37200   46499   0.400010765 1524670052.1594 2.852804184
3.0739262104
46500   55799   0.500010729 1524670055.909613.7502081394
3.9270553589
55800   65099   0.600010753 1524670060.476544.5669286251
4.7571902275
65100   74399   0.700010777 1524670066.0941 5.6175630093
5.7428796291
74400   83699   0.800010741 1524670072.538866.44475317
6.5761549473
83700   92998   0.900010765 1524670079.990727.4518604279
7.4606924057
92999   102298  1   1524670087.710667.7199423313
8.2424075603



old code

do ipoint = 0, istart-1
  !c output time cost, use 1 processor to test
  if (b_enable_output .and. rank == 0) then
if (mod(ipoint,iprogress) == 0 .or. ipoint ==
istart-1) then
  !write(*,'(f3.1,1x)',advance="no") (ipoint+1.0)/istart
  write(*,*) ipoint,
(ipoint+1.0)/istart,"time",MPI_Wtime()
end if
  end if

  call DMSetLabelValue(dmda_flow%da,"cid_lg2g",ipoint, &
   ipoint+1,ierr)
  CHKERRQ(ierr)
end do


new code

call DMCreateLabel(dmda_flow%da,'cid_lg2g',ierr)
CHKERRQ(ierr)

call DMGetLabel(dmda_flow%da,'cid_lg2g',label, ierr)
CHKERRQ(ierr)

do ipoint = 0, istart-1
  !c output time cost, use 1 processor to test
  if (b_enable_output .and. rank == 0) then
if (mod(ipoint,iprogress) == 0 .or. ipoint ==
istart-1) then
  !writ

Re: [petsc-users] DMSetLabelValue takes a lot of time for large domain

2018-04-25 Thread Danyang Su


On 2018-04-25 09:47 AM, Matthew Knepley wrote:
On Wed, Apr 25, 2018 at 12:40 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi Matthew,

In the worst case, every node/cell may have different label.

Do not use Label for this. Its not an appropriate thing. If every cell 
is different, just use the cell number.
Labels are for mapping a relatively small number of keys (like 
material IDs) to sets of points (cells, vertices, etc.)

Its not a great data structure for a permutation.
Yes. If there is small number of keys, it runs very fast, even for more 
than one million DMSetLabelValue calls. The performance just 
deteriorates as the number of keys increases.


I cannot get avoid of DMSetLabelValue as node/cell index of original 
mesh is needed for the previous input file that uses some of global 
node/cell index to set value. But if I can get the natural order of 
nodes/cells from DMPlex, I can discard the use of DMSetLabelValue. Is 
there any function can do this job?


Thanks,

Danyang


However, I still do not believe these numbers. The old code does a 
string comparison every time. I will setup a test.


   Matt

Below is one of the worst scenario with 102299 nodes and 102299
different labels for test. I found the time cost increase during
the loop. The first 9300 loop takes least time (<0.5) while the
last 9300 loops takes much more time (>7.7), as shown below. If I
use larger mesh with >1 million nodes, it runs very very slowly in
this part. The PETSc is configured with optimization on.

Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran
--download-mpich --download-scalapack --download-parmetis
--download-metis --download-ptscotch --download-fblaslapack
--download-hypre --download-superlu_dist --download-hdf5=yes
--download-ctetgen --with-debugging=0 COPTFLAGS="-O3 -march=native
-mtune=native" CXXOPTFLAGS="-O3 -march=native -mtune=native"
FOPTFLAGS="-O3 -march=native -mtune=native"

istart  iendprogressCPU_Timetime cost - old (sec)   
time cost
- new (sec)
0   92990   1524670045.51166

930018599   0.100010753 1524670045.996050.4843890667
0.497246027
18600   27899   0.200010747 1524670047.326351.330302
1.3820912838
27900   37199   0.300010741 1524670049.3066 1.9802515507
2.2439446449
37200   46499   0.400010765 1524670052.1594 2.852804184 
3.0739262104
46500   55799   0.500010729 1524670055.909613.7502081394
3.9270553589
55800   65099   0.600010753 1524670060.476544.5669286251
4.7571902275
65100   74399   0.700010777 1524670066.0941 5.6175630093
5.7428796291
74400   83699   0.800010741 1524670072.538866.44475317  
6.5761549473
83700   92998   0.900010765 1524670079.990727.4518604279
7.4606924057
92999   102298  1   1524670087.710667.7199423313
8.2424075603



old code

    do ipoint = 0, istart-1
  !c output time cost, use 1 processor to test
  if (b_enable_output .and. rank == 0) then
    if (mod(ipoint,iprogress) == 0 .or. ipoint ==
istart-1) then
  !write(*,'(f3.1,1x)',advance="no") (ipoint+1.0)/istart
  write(*,*) ipoint,
(ipoint+1.0)/istart,"time",MPI_Wtime()
    end if
  end if

  call DMSetLabelValue(dmda_flow%da,"cid_lg2g",ipoint, &
   ipoint+1,ierr)
  CHKERRQ(ierr)
    end do


new code

    call DMCreateLabel(dmda_flow%da,'cid_lg2g',ierr)
    CHKERRQ(ierr)

    call DMGetLabel(dmda_flow%da,'cid_lg2g',label, ierr)
    CHKERRQ(ierr)

    do ipoint = 0, istart-1
  !c output time cost, use 1 processor to test
  if (b_enable_output .and. rank == 0) then
    if (mod(ipoint,iprogress) == 0 .or. ipoint ==
istart-1) then
  !write(*,'(f3.1,1x)',advance="no") (ipoint+1.0)/istart
  write(*,*) ipoint,
(ipoint+1.0)/istart,"time",MPI_Wtime()
    end if
  end if

  call DMLabelSetValue(label,ipoint,ipoint+1,ierr)
  CHKERRQ(ierr)
    end do

Thanks,

Danyang

On 2018-04-25 03:16 AM, Matthew Knepley wrote:

On Tue, Apr 24, 2018 at 11:57 PM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:

Hi All,

I use DMPlex in unstructured grid code and recently found
DMSetLabelValue takes a lot of time for large problem, e.g.,
num. of cells > 1 million. In my c

Re: [petsc-users] DMSetLabelValue takes a lot of time for large domain

2018-04-25 Thread Danyang Su


Hi Matthew,

In the worst case, every node/cell may have different label.

Below is one of the worst scenario with 102299 nodes and 102299 
different labels for test. I found the time cost increase during the 
loop. The first 9300 loop takes least time (<0.5) while the last 9300 
loops takes much more time (>7.7), as shown below. If I use larger mesh 
with >1 million nodes, it runs very very slowly in this part. The PETSc 
is configured with optimization on.


Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran 
--download-mpich --download-scalapack --download-parmetis 
--download-metis --download-ptscotch --download-fblaslapack 
--download-hypre --download-superlu_dist --download-hdf5=yes 
--download-ctetgen --with-debugging=0 COPTFLAGS="-O3 -march=native 
-mtune=native" CXXOPTFLAGS="-O3 -march=native -mtune=native" 
FOPTFLAGS="-O3 -march=native -mtune=native"


istart 	iend 	progress 	CPU_Time 	time cost - old (sec) 	time cost - new 
(sec)

0   92990   1524670045.51166

930018599   0.100010753 1524670045.996050.4843890667
0.497246027
18600   27899   0.200010747 1524670047.326351.330302
1.3820912838
27900   37199   0.300010741 1524670049.3066 1.9802515507
2.2439446449
37200   46499   0.400010765 1524670052.1594 2.852804184 
3.0739262104
46500   55799   0.500010729 1524670055.909613.7502081394
3.9270553589
55800   65099   0.600010753 1524670060.476544.5669286251
4.7571902275
65100   74399   0.700010777 1524670066.0941 5.6175630093
5.7428796291
74400   83699   0.800010741 1524670072.538866.44475317  
6.5761549473
83700   92998   0.900010765 1524670079.990727.4518604279
7.4606924057
92999   102298  1   1524670087.710667.71994233138.2424075603



old code

    do ipoint = 0, istart-1
  !c output time cost, use 1 processor to test
  if (b_enable_output .and. rank == 0) then
    if (mod(ipoint,iprogress) == 0 .or. ipoint == istart-1) then
  !write(*,'(f3.1,1x)',advance="no") (ipoint+1.0)/istart
  write(*,*) ipoint, (ipoint+1.0)/istart,"time",MPI_Wtime()
    end if
  end if

  call DMSetLabelValue(dmda_flow%da,"cid_lg2g",ipoint, &
   ipoint+1,ierr)
  CHKERRQ(ierr)
    end do


new code

    call DMCreateLabel(dmda_flow%da,'cid_lg2g',ierr)
    CHKERRQ(ierr)

    call DMGetLabel(dmda_flow%da,'cid_lg2g',label, ierr)
    CHKERRQ(ierr)

    do ipoint = 0, istart-1
  !c output time cost, use 1 processor to test
  if (b_enable_output .and. rank == 0) then
    if (mod(ipoint,iprogress) == 0 .or. ipoint == istart-1) then
  !write(*,'(f3.1,1x)',advance="no") (ipoint+1.0)/istart
  write(*,*) ipoint, (ipoint+1.0)/istart,"time",MPI_Wtime()
    end if
  end if

  call DMLabelSetValue(label,ipoint,ipoint+1,ierr)
  CHKERRQ(ierr)
    end do

Thanks,

Danyang

On 2018-04-25 03:16 AM, Matthew Knepley wrote:
On Tue, Apr 24, 2018 at 11:57 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi All,

I use DMPlex in unstructured grid code and recently found
DMSetLabelValue takes a lot of time for large problem, e.g., num.
of cells > 1 million. In my code, I use


I read your code wrong. For large loop, you should not use the 
convenience function. You should use


DMPlexCreateFromCellList ()


DMGetLabel(dm, name, )


Loop over all cells/nodes{

DMSetLabelValue


Replace this by DMLabelSetValue(label, point, val)

}

DMPlexDistribute

The code works fine except DMSetLabelValue takes a lot of time for
large problem. I use DMSetLabelValue to set material id for all
the nodes or cells so that each subdomain has a copy of material
id. Is there any other functions that can be used more efficient,
e.g. set labels by array, not 1 by 1?


That should take much less time.

  Thanks,

     Matt

Thanks,

Danyang




--
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.caam.rice.edu/%7Emk51/>

[petsc-users] DMSetLabelValue takes a lot of time for large domain

2018-04-24 Thread Danyang Su


Hi All,

I use DMPlex in unstructured grid code and recently found 
DMSetLabelValue takes a lot of time for large problem, e.g., num. of 
cells > 1 million. In my code, I use


DMPlexCreateFromCellList ()

Loop over all cells/nodes{

DMSetLabelValue

}

DMPlexDistribute

The code works fine except DMSetLabelValue takes a lot of time for large 
problem. I use DMSetLabelValue to set material id for all the nodes or 
cells so that each subdomain has a copy of material id. Is there any 
other functions that can be used more efficient, e.g. set labels by 
array, not 1 by 1?


Thanks,

Danyang

Re: [petsc-users] [petsc-maint] how to check if cell is local owned in DMPlex

2018-03-08 Thread Danyang Su


On 18-03-07 03:54 PM, Jed Brown wrote:

Danyang Su <danyang...@gmail.com> writes:


Based on my test, this function works fine using current PETSc-dev
version, but I cannot get it compiled correctly using other versions for
Fortran code, as mentioned in the previous emails. I asked this question
because some of the clusters we use do not have the PETSc-dev version
and it takes time get staff reinstall another version.

You can install PETSc yourself as a normal user.

Thanks, Jed. This is the way I will follow.

Danyang

Re: [petsc-users] [petsc-maint] how to check if cell is local owned in DMPlex

2018-03-07 Thread Danyang Su


Hi All,

Thanks again for all your help during my code development and it turns 
out the DMPlex for unstructured grid works pretty well. I have another 
question.


 Is there any alternative method to get number of leaves without using 
PetscSFGetGraph?


Based on my test, this function works fine using current PETSc-dev 
version, but I cannot get it compiled correctly using other versions for 
Fortran code, as mentioned in the previous emails. I asked this question 
because some of the clusters we use do not have the PETSc-dev version 
and it takes time get staff reinstall another version.


Thanks,

Danyang

On 18-03-05 11:50 AM, Smith, Barry F. wrote:

MatSolverPackage

became MatSolverType




On Mar 5, 2018, at 1:35 PM, Danyang Su <danyang...@gmail.com> wrote:

Hi Barry and Matt,

The compiling problem should be caused by the PETSc version installed on my 
computer. When updated to PETSc-Dev version, the ex1f example works fine. 
However, I cannot compile this example under PETSc-3.8.3 version.

After updating to PETSc-Dev verison, I encounter another compiling problem in 
my code.

 MatSolverPackage :: solver_pkg_flow
 1
Error: Unclassifiable statement at (1)

Including petscmat.h or petscpc.h does not help to solve this problem. I can 
rewrite this part to get rid of this. But I would rather to keep this if there 
is alternative way to go. What is the head file should I include in order to 
use MatSolverPackage?

Thanks,

Danyang


On 18-03-04 11:15 AM, Smith, Barry F. wrote:

   See src/vec/is/sf/examples/tutorials/ex1f.F90 in the master branch of the 
PETSc git repository

   BTW:

git grep -i petscsfgetgraph

will show every use of the function in the source code. Very useful tool

Barry



On Mar 4, 2018, at 1:05 PM, Danyang Su <danyang...@gmail.com> wrote:



On 18-03-04 08:08 AM, Matthew Knepley wrote:

On Fri, Mar 2, 2018 at 3:22 PM, Danyang Su <danyang...@gmail.com> wrote:
Hi Matt,
I use the latest Fortran style in PETSc 3.8. Enclosed are the PETSc 
configuration, code compiling log and the function that causes compiling error. 
The compiling error happens after I include petscsf.h in the following section. 
I didn't find petscsf.h in petsc/finclude/ folder so I use the head file in the 
'include' folder and this seems not allowed.

I apologize for taking so long. The PetscSF definitions are in

#include 

Hi Matt,

After including
#include 
   use petscis

I still get error saying undefined reference to `petscsfgetgraph_'

Did I miss any other head file?

Thanks,

Danyang

You are correct that they should be moved out.

   Thanks,

  Matt

#ifdef PETSC_V3_8_X

#include 
#include 
#include 
#include 
   use petscsys
   use petscdmplex
   use petscsf

#endif

Thanks,

Danyang


On 18-03-02 12:08 PM, Matthew Knepley wrote:

On Fri, Mar 2, 2018 at 3:00 PM, Danyang Su <danyang...@gmail.com> wrote:
On 18-03-02 10:58 AM, Matthew Knepley wrote:

On Fri, Mar 2, 2018 at 1:41 PM, Danyang Su <danyang...@gmail.com> wrote:

On 18-02-19 03:30 PM, Matthew Knepley wrote:

On Mon, Feb 19, 2018 at 3:11 PM, Danyang Su <danyang...@gmail.com> wrote:
Hi Matt,

Would you please let me know how to check if a cell is local owned? When overlap 
is 0 in DMPlexDistribute, all the cells are local owned. How about overlap > 0? 
It sounds like impossible to check by node because a cell can be local owned even 
if none of the nodes in this cell is local owned.

If a cell is in the PetscSF, then it is not locally owned. The local nodes in 
the SF are sorted, so I use
PetscFindInt 
(http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html).

Hi Matt,

Would you please give me a little more about how to mark the ghost cells when 
overlap > 0? What do you mean a cell is in the PetscSF? I use PetscSFView to 
export the graph (original mesh file pile.vtk) and it exports all the cells, 
including the ghost cells (PETScSFView.txt).

Yes, I will send you some sample code when I get time. The first problem is 
that you are looking at a different PetscSF. This looks like the
one returned by DMPlexDistribute(). This is mapping the serial mesh to the 
parallel mesh. You want

   
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetPointSF.html

Then you can look at

   
https://bitbucket.org/petsc/petsc/src/1788fc36644e622df8cb1a0de85676ccc5af0239/src/dm/impls/plex/plexsubmesh.c?at=master=file-view-default#plexsubmesh.c-683

I get the pointSF, get out the list of leaves, and find points in it using 
PetscFindInt()

Hi Matt,
By using the local dm, I can get the PetscSF I want, as shown below. Now I need to get the number 
of ghost cells or local cells (here 4944) or number of leaves (here 825) for each processor. I try 
to use PetscSFGetGraph to get number of leaves in Fortran. After including "petscsf.h", I 
got compilation error saying "You need a ISO C conforming compiler to use the glibc

Re: [petsc-users] [petsc-maint] how to check if cell is local owned in DMPlex

2018-03-05 Thread Danyang Su


Hi Barry and Matt,

The compiling problem should be caused by the PETSc version installed on 
my computer. When updated to PETSc-Dev version, the ex1f example works 
fine. However, I cannot compile this example under PETSc-3.8.3 version.


After updating to PETSc-Dev verison, I encounter another compiling 
problem in my code.


    MatSolverPackage :: solver_pkg_flow
    1
Error: Unclassifiable statement at (1)

Including petscmat.h or petscpc.h does not help to solve this problem. I 
can rewrite this part to get rid of this. But I would rather to keep 
this if there is alternative way to go. What is the head file should I 
include in order to use MatSolverPackage?


Thanks,

Danyang


On 18-03-04 11:15 AM, Smith, Barry F. wrote:

   See src/vec/is/sf/examples/tutorials/ex1f.F90 in the master branch of the 
PETSc git repository

   BTW:

git grep -i petscsfgetgraph

will show every use of the function in the source code. Very useful tool

Barry



On Mar 4, 2018, at 1:05 PM, Danyang Su <danyang...@gmail.com> wrote:



On 18-03-04 08:08 AM, Matthew Knepley wrote:

On Fri, Mar 2, 2018 at 3:22 PM, Danyang Su <danyang...@gmail.com> wrote:
Hi Matt,
I use the latest Fortran style in PETSc 3.8. Enclosed are the PETSc 
configuration, code compiling log and the function that causes compiling error. 
The compiling error happens after I include petscsf.h in the following section. 
I didn't find petscsf.h in petsc/finclude/ folder so I use the head file in the 
'include' folder and this seems not allowed.

I apologize for taking so long. The PetscSF definitions are in

#include 

Hi Matt,

After including
#include 
   use petscis

I still get error saying undefined reference to `petscsfgetgraph_'

Did I miss any other head file?

Thanks,

Danyang

You are correct that they should be moved out.

   Thanks,

  Matt

#ifdef PETSC_V3_8_X

#include 
#include 
#include 
#include 
   use petscsys
   use petscdmplex
   use petscsf

#endif

Thanks,

Danyang


On 18-03-02 12:08 PM, Matthew Knepley wrote:

On Fri, Mar 2, 2018 at 3:00 PM, Danyang Su <danyang...@gmail.com> wrote:
On 18-03-02 10:58 AM, Matthew Knepley wrote:

On Fri, Mar 2, 2018 at 1:41 PM, Danyang Su <danyang...@gmail.com> wrote:

On 18-02-19 03:30 PM, Matthew Knepley wrote:

On Mon, Feb 19, 2018 at 3:11 PM, Danyang Su <danyang...@gmail.com> wrote:
Hi Matt,

Would you please let me know how to check if a cell is local owned? When overlap 
is 0 in DMPlexDistribute, all the cells are local owned. How about overlap > 0? 
It sounds like impossible to check by node because a cell can be local owned even 
if none of the nodes in this cell is local owned.

If a cell is in the PetscSF, then it is not locally owned. The local nodes in 
the SF are sorted, so I use
PetscFindInt 
(http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html).

Hi Matt,

Would you please give me a little more about how to mark the ghost cells when 
overlap > 0? What do you mean a cell is in the PetscSF? I use PetscSFView to 
export the graph (original mesh file pile.vtk) and it exports all the cells, 
including the ghost cells (PETScSFView.txt).

Yes, I will send you some sample code when I get time. The first problem is 
that you are looking at a different PetscSF. This looks like the
one returned by DMPlexDistribute(). This is mapping the serial mesh to the 
parallel mesh. You want

   
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetPointSF.html

Then you can look at

   
https://bitbucket.org/petsc/petsc/src/1788fc36644e622df8cb1a0de85676ccc5af0239/src/dm/impls/plex/plexsubmesh.c?at=master=file-view-default#plexsubmesh.c-683

I get the pointSF, get out the list of leaves, and find points in it using 
PetscFindInt()

Hi Matt,
By using the local dm, I can get the PetscSF I want, as shown below. Now I need to get the number 
of ghost cells or local cells (here 4944) or number of leaves (here 825) for each processor. I try 
to use PetscSFGetGraph to get number of leaves in Fortran. After including "petscsf.h", I 
got compilation error saying "You need a ISO C conforming compiler to use the glibc 
headers". Is there any alternative way to do this? I do not need the ghost-neighbor mapping, 
but just the number of local owned cells.

Also, make sure you are using the latest Fortran style for PETSc:

   
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/UsingFortran.html
  
   [0] Number of roots=11449, leaves=825, remote ranks=1

   [0] 4944 <- (1,0)
   [0] 4945 <- (1,28)
   [0] 4946 <- (1,56)
...
   [1] Number of roots=11695, leaves=538, remote ranks=1
   [1] 5056 <- (0,21)
   [1] 5057 <- (0,43)
   [1] 5058 <- (0,65)
   [1] 5059 <- (0,87)

In file included from /usr/include/features.h:375:0,
  from /usr/include/stdio.h:28,
  from /home/dsu/Soft/PETSc/petsc-3.8.3/include/petscsys.h:175

Re: [petsc-users] how to check if cell is local owned in DMPlex

2018-03-04 Thread Danyang Su




On 18-03-04 08:08 AM, Matthew Knepley wrote:
On Fri, Mar 2, 2018 at 3:22 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi Matt,

I use the latest Fortran style in PETSc 3.8. Enclosed are the
PETSc configuration, code compiling log and the function that
causes compiling error. The compiling error happens after I
include petscsf.h in the following section. I didn't find
petscsf.h in petsc/finclude/ folder so I use the head file in the
'include' folder and this seems not allowed.


I apologize for taking so long. The PetscSF definitions are in

#include 

Hi Matt,

After including
#include 
  use petscis

I still get error saying undefined reference to `petscsfgetgraph_'

Did I miss any other head file?

Thanks,

Danyang


You are correct that they should be moved out.

  Thanks,

 Matt

#ifdef PETSC_V3_8_X

#include 
#include 
#include 
#include 
  use petscsys
  use petscdmplex
  use petscsf

#endif

Thanks,

Danyang


On 18-03-02 12:08 PM, Matthew Knepley wrote:

On Fri, Mar 2, 2018 at 3:00 PM, Danyang Su <danyang...@gmail.com
<mailto:danyang...@gmail.com>> wrote:

On 18-03-02 10:58 AM, Matthew Knepley wrote:


On Fri, Mar 2, 2018 at 1:41 PM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:


On 18-02-19 03:30 PM, Matthew Knepley wrote:

    On Mon, Feb 19, 2018 at 3:11 PM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:

Hi Matt,

Would you please let me know how to check if a cell
is local owned? When overlap is 0 in
DMPlexDistribute, all the cells are local owned.
How about overlap > 0? It sounds like impossible to
check by node because a cell can be local owned
even if none of the nodes in this cell is local owned.


If a cell is in the PetscSF, then it is not locally
owned. The local nodes in the SF are sorted, so I use
PetscFindInt

(http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html

<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html>).

Hi Matt,

Would you please give me a little more about how to mark
the ghost cells when overlap > 0? What do you mean a
cell is in the PetscSF? I use PetscSFView to export the
graph (original mesh file pile.vtk) and it exports all
the cells, including the ghost cells (PETScSFView.txt).


Yes, I will send you some sample code when I get time. The
first problem is that you are looking at a different
PetscSF. This looks like the
one returned by DMPlexDistribute(). This is mapping the
serial mesh to the parallel mesh. You want


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetPointSF.html

<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetPointSF.html>

Then you can look at


https://bitbucket.org/petsc/petsc/src/1788fc36644e622df8cb1a0de85676ccc5af0239/src/dm/impls/plex/plexsubmesh.c?at=master=file-view-default#plexsubmesh.c-683

<https://bitbucket.org/petsc/petsc/src/1788fc36644e622df8cb1a0de85676ccc5af0239/src/dm/impls/plex/plexsubmesh.c?at=master=file-view-default#plexsubmesh.c-683>

I get the pointSF, get out the list of leaves, and find
points in it using PetscFindInt()

Hi Matt,
By using the local dm, I can get the PetscSF I want, as shown
below. Now I need to get the number of ghost cells or local
cells (here 4944) or number of leaves (here 825) for each
processor. I try to use PetscSFGetGraph to get number of
leaves in Fortran. After including "petscsf.h", I got
compilation error saying "You need a ISO C conforming
compiler to use the glibc headers". Is there any alternative
way to do this? I do not need the ghost-neighbor mapping, but
just the number of local owned cells.


Also, make sure you are using the latest Fortran style for PETSc:


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/UsingFortran.html

<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/UsingFortran.html>

  [0] Number of roots=11449, leaves=825, remote ranks=1
  [0] 4944 <- (1,0)
  [0] 4945 <- (1,28)
  [0] 4946 <- (1,56)
...
  [1] Number of roots=11695, leaves=538, remote ranks=1
  [1] 5056 <- (0,21)
  [1] 5057 <- (0,43)
  [1] 5058 <- (0,65)
  [1] 5059 <- (0,87)

In file included from /usr/inclu

Re: [petsc-users] how to check if cell is local owned in DMPlex

2018-03-02 Thread Danyang Su


Hi Matt,

I use the latest Fortran style in PETSc 3.8. Enclosed are the PETSc 
configuration, code compiling log and the function that causes compiling 
error. The compiling error happens after I include petscsf.h in the 
following section. I didn't find petscsf.h in petsc/finclude/ folder so 
I use the head file in the 'include' folder and this seems not allowed.


#ifdef PETSC_V3_8_X

#include 
#include 
#include 
#include 
  use petscsys
  use petscdmplex
  use petscsf

#endif

Thanks,

Danyang


On 18-03-02 12:08 PM, Matthew Knepley wrote:
On Fri, Mar 2, 2018 at 3:00 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


On 18-03-02 10:58 AM, Matthew Knepley wrote:


On Fri, Mar 2, 2018 at 1:41 PM, Danyang Su <danyang...@gmail.com
<mailto:danyang...@gmail.com>> wrote:


On 18-02-19 03:30 PM, Matthew Knepley wrote:

On Mon, Feb 19, 2018 at 3:11 PM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:

Hi Matt,

Would you please let me know how to check if a cell is
local owned? When overlap is 0 in DMPlexDistribute, all
the cells are local owned. How about overlap > 0? It
sounds like impossible to check by node because a cell
can be local owned even if none of the nodes in this
cell is local owned.


If a cell is in the PetscSF, then it is not locally owned.
The local nodes in the SF are sorted, so I use
PetscFindInt

(http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html

<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html>).

Hi Matt,

Would you please give me a little more about how to mark the
ghost cells when overlap > 0? What do you mean a cell is in
the PetscSF? I use PetscSFView to export the graph (original
mesh file pile.vtk) and it exports all the cells, including
the ghost cells (PETScSFView.txt).


Yes, I will send you some sample code when I get time. The first
problem is that you are looking at a different PetscSF. This
looks like the
one returned by DMPlexDistribute(). This is mapping the serial
mesh to the parallel mesh. You want


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetPointSF.html

<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetPointSF.html>

Then you can look at


https://bitbucket.org/petsc/petsc/src/1788fc36644e622df8cb1a0de85676ccc5af0239/src/dm/impls/plex/plexsubmesh.c?at=master=file-view-default#plexsubmesh.c-683

<https://bitbucket.org/petsc/petsc/src/1788fc36644e622df8cb1a0de85676ccc5af0239/src/dm/impls/plex/plexsubmesh.c?at=master=file-view-default#plexsubmesh.c-683>

I get the pointSF, get out the list of leaves, and find points in
it using PetscFindInt()

Hi Matt,
By using the local dm, I can get the PetscSF I want, as shown
below. Now I need to get the number of ghost cells or local cells
(here 4944) or number of leaves (here 825) for each processor. I
try to use PetscSFGetGraph to get number of leaves in Fortran.
After including "petscsf.h", I got compilation error saying "You
need a ISO C conforming compiler to use the glibc headers". Is
there any alternative way to do this? I do not need the
ghost-neighbor mapping, but just the number of local owned cells.


Also, make sure you are using the latest Fortran style for PETSc:

http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/UsingFortran.html

  [0] Number of roots=11449, leaves=825, remote ranks=1
  [0] 4944 <- (1,0)
  [0] 4945 <- (1,28)
  [0] 4946 <- (1,56)
...
  [1] Number of roots=11695, leaves=538, remote ranks=1
  [1] 5056 <- (0,21)
  [1] 5057 <- (0,43)
  [1] 5058 <- (0,65)
  [1] 5059 <- (0,87)

In file included from /usr/include/features.h:375:0,
 from /usr/include/stdio.h:28,
 from
/home/dsu/Soft/PETSc/petsc-3.8.3/include/petscsys.h:175,
 from
/home/dsu/Soft/PETSc/petsc-3.8.3/include/petscsf.h:7,
 from ../../solver/solver_ddmethod.F90:4837:
/usr/include/x86_64-linux-gnu/sys/cdefs.h:30:3: error: #error "You
need a ISO C conforming compile\
r to use the glibc headers"
 # error "You need a ISO C conforming compiler to use the glibc
headers"


Can you send this to petsc-ma...@mcs.anl.gov 
<mailto:petsc-ma...@mcs.anl.gov>? It looks like a build problem that 
can be fixed.


  Thanks,

    Matt

Thanks,

Danyang


  Thanks,

    Matt

Thanks,

Danyang


  Thanks,

    Matt

Thanks,

Danyang

--

Re: [petsc-users] how to check if cell is local owned in DMPlex

2018-03-02 Thread Danyang Su




On 18-03-02 10:58 AM, Matthew Knepley wrote:
On Fri, Mar 2, 2018 at 1:41 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:



On 18-02-19 03:30 PM, Matthew Knepley wrote:

On Mon, Feb 19, 2018 at 3:11 PM, Danyang Su <danyang...@gmail.com
<mailto:danyang...@gmail.com>> wrote:

Hi Matt,

Would you please let me know how to check if a cell is local
owned? When overlap is 0 in DMPlexDistribute, all the cells
are local owned. How about overlap > 0? It sounds like
impossible to check by node because a cell can be local owned
even if none of the nodes in this cell is local owned.


If a cell is in the PetscSF, then it is not locally owned. The
local nodes in the SF are sorted, so I use
PetscFindInt

(http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html

<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html>).

Hi Matt,

Would you please give me a little more about how to mark the ghost
cells when overlap > 0? What do you mean a cell is in the PetscSF?
I use PetscSFView to export the graph (original mesh file
pile.vtk) and it exports all the cells, including the ghost cells
(PETScSFView.txt).


Yes, I will send you some sample code when I get time. The first 
problem is that you are looking at a different PetscSF. This looks 
like the
one returned by DMPlexDistribute(). This is mapping the serial mesh to 
the parallel mesh. You want


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetPointSF.html

Then you can look at

https://bitbucket.org/petsc/petsc/src/1788fc36644e622df8cb1a0de85676ccc5af0239/src/dm/impls/plex/plexsubmesh.c?at=master=file-view-default#plexsubmesh.c-683

I get the pointSF, get out the list of leaves, and find points in it 
using PetscFindInt()

Hi Matt,
By using the local dm, I can get the PetscSF I want, as shown below. Now 
I need to get the number of ghost cells or local cells (here 4944) or 
number of leaves (here 825) for each processor. I try to use 
PetscSFGetGraph to get number of leaves in Fortran. After including 
"petscsf.h", I got compilation error saying "You need a ISO C conforming 
compiler to use the glibc headers". Is there any alternative way to do 
this? I do not need the ghost-neighbor mapping, but just the number of 
local owned cells.


  [0] Number of roots=11449, leaves=825, remote ranks=1
  [0] 4944 <- (1,0)
  [0] 4945 <- (1,28)
  [0] 4946 <- (1,56)
...
  [1] Number of roots=11695, leaves=538, remote ranks=1
  [1] 5056 <- (0,21)
  [1] 5057 <- (0,43)
  [1] 5058 <- (0,65)
  [1] 5059 <- (0,87)

In file included from /usr/include/features.h:375:0,
 from /usr/include/stdio.h:28,
 from 
/home/dsu/Soft/PETSc/petsc-3.8.3/include/petscsys.h:175,

 from /home/dsu/Soft/PETSc/petsc-3.8.3/include/petscsf.h:7,
 from ../../solver/solver_ddmethod.F90:4837:
/usr/include/x86_64-linux-gnu/sys/cdefs.h:30:3: error: #error "You need 
a ISO C conforming compile\

r to use the glibc headers"
 # error "You need a ISO C conforming compiler to use the glibc headers"

Thanks,

Danyang


  Thanks,

    Matt

Thanks,

Danyang


  Thanks,

    Matt

Thanks,

Danyang

-- 
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.caam.rice.edu/%7Emk51/>





--
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.caam.rice.edu/%7Emk51/>

Re: [petsc-users] how to check if cell is local owned in DMPlex

2018-03-02 Thread Danyang Su

On 18-03-02 10:58 AM, Matthew Knepley wrote:
On Fri, Mar 2, 2018 at 1:41 PM, Danyang Su <danyang...@gmail.com
<mailto:danyang...@gmail.com>> wrote:

On 18-02-19 03:30 PM, Matthew Knepley wrote:

On Mon, Feb 19, 2018 at 3:11 PM, Danyang Su <danyang...@gmail.com
<mailto:danyang...@gmail.com>> wrote:

Hi Matt,

Would you please let me know how to check if a cell is local
owned? When overlap is 0 in DMPlexDistribute, all the cells
are local owned. How about overlap > 0? It sounds like
impossible to check by node because a cell can be local owned
even if none of the nodes in this cell is local owned.

If a cell is in the PetscSF, then it is not locally owned. The
local nodes in the SF are sorted, so I use
PetscFindInt

(http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html

<http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscFindInt.html>).

Hi Matt,

Would you please give me a little more about how to mark the ghost
cells when overlap > 0? What do you mean a cell is in the PetscSF?
I use PetscSFView to export the graph (original mesh file
pile.vtk) and it exports all the cells, including the ghost cells
(PETScSFView.txt).

Yes, I will send you some sample code when I get time. The first
problem is that you are looking at a different PetscSF. This looks
like the
one returned by DMPlexDistribute(). This is mapping the serial mesh to
the parallel mesh. You want

http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetPointSF.html

Then you can look at

https://bitbucket.org/petsc/petsc/src/1788fc36644e622df8cb1a0de85676ccc5af0239/src/dm/impls/plex/plexsubmesh.c?at=master=file-view-default#plexsubmesh.c-683

I get the pointSF, get out the list of leaves, and find points in it
using PetscFindInt()
Thanks Matt. I will try to figure it out based on your provided link and
will let you know if I get it work.

Danyang

Thanks,

Matt

Thanks,

Danyang

Thanks,

Matt

Thanks,

Danyang

--
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.caam.rice.edu/%7Emk51/>

--
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.caam.rice.edu/%7Emk51/>

Re: [petsc-users] object name overwritten in VecView

2018-02-28 Thread Danyang Su


Hi Barry and Matt,

Thanks for your quick response. Considering the output performance, as 
well as the long-term plan of PETSc development, which format would you 
suggest? I personally prefer the data format that can be post-processed 
by Paraview as our sequential code (written without PETSc) is also using 
Paraview compatible data format. XDMF sounds promising as suggested by Matt.


Thanks,

Danyang

On 18-02-28 08:17 AM, Smith, Barry F. wrote:


  It turns out the fix is really easy. Here is a patch.

  Apply it with

    patch -p1 < barry-vtk.patch

  then do

    make gnumake

  all in $PETSC_DIR



> On Feb 28, 2018, at 9:07 AM, Danyang Su <danyang...@gmail.com> wrote:
>
> Hi Matt，
>
> Thanks for your suggestion and I will use xmf instead.
>
> Regards,
>
> Danyang
>
> On February 28, 2018 3:58:08 AM PST, Matthew Knepley 
<knep...@gmail.com> wrote:
> On Wed, Feb 28, 2018 at 12:39 AM, Smith, Barry F. 
<bsm...@mcs.anl.gov> wrote:

>
>   Matt,
>
>   I have confirmed this is reproducible and a bug. The problem 
arises because

>
> frame #0: 0x00010140625a 
libpetsc.3.8.dylib`PetscViewerVTKAddField_VTK(viewer=0x7fe66760c750, 
dm=0x7fe668810820, 
PetscViewerVTKWriteFunction=(libpetsc.3.8.dylib`DMPlexVTKWriteAll at 
plexvtk.c:633), fieldtype=PETSC_VTK_POINT_FIELD, 
vec=0x7fe66880ee20) at vtkv.c:140
> frame #1: 0x000101404e6e 
libpetsc.3.8.dylib`PetscViewerVTKAddField(viewer=0x7fe66760c750, 
dm=0x7fe668810820, 
PetscViewerVTKWriteFunction=(libpetsc.3.8.dylib`DMPlexVTKWriteAll at 
plexvtk.c:633), fieldtype=PETSC_VTK_POINT_FIELD, 
vec=0x7fe66880ee20) at vtkv.c:46
> frame #2: 0x000101e0b7c3 
libpetsc.3.8.dylib`VecView_Plex_Local(v=0x7fe66880ee20, 
viewer=0x7fe66760c750) at plex.c:301
> frame #3: 0x000101e0ead7 
libpetsc.3.8.dylib`VecView_Plex(v=0x7fe66880e820, 
viewer=0x7fe66760c750) at plex.c:348

>
> keeps a linked list of vectors that are to be viewed and the vectors 
are the same Vec because they are obtained with DMGetLocalVector().

>
> The safest fix is to have PetscViewerVTKAddField_VTK() do a 
VecDuplicate() on the vector passed in and store that in the linked 
list instead of just storing a pointer to the passed in vector (which 
might and can be overwritten before all the linked vectors are 
actually stored).

>
> Danyang,
>
> Barry is right, and the bug can be fixed the way he says. However, 
this points out why VTK is bad format. I think a better choice is

> to use HDF5 and XDMF. For example, in my code now I always use
>
>   DMVIewFromOptions(dm, NULL, "-dm_view");
>
> and then later (perhaps several times)
>
>   VecViewFromOptions(u, NULL, "-u_vec_view")
>   VecViewFromOptions(v, NULL, "-v_vec_view")
>
> and then on the command line
>
>   -dm_view hdf5:test.h5 -u_vec_view hdf5:test.h5::append -v_vec_view 
hdf5:test.h5::append

>
> which produces a file
>
>   test.h5
>
> Then I run
>
>   $PETSC_DIR/bin/petsc_gen_xdmf.py test.h5
>
> which produces another file
>
>   test.xmf
>
> This can be loaded by Paraview for visualization.
>
>   Thanks,
>
>  Matt
>
>
>
>   Barry
>
>
>
> > On Feb 27, 2018, at 10:44 PM, Danyang Su <danyang...@gmail.com> wrote:
> >
> > Hi All,
> >
> > How to set different object names when using multiple VecView? I 
try to use PetscObjectSetName with multiple output, but the object 
name is overwritten by the last one.

> >
> > As shown below, as well as the enclosed files as example, the 
vector name in sol.vtk is vec_v for both vector u and v.

> >
> >  call PetscViewerCreate(PETSC_COMM_WORLD, viewer, 
ierr);CHKERRA(ierr)
> >  call PetscViewerSetType(viewer, PETSCVIEWERVTK, 
ierr);CHKERRA(ierr)
> >  call PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK, 
ierr);CHKERRA(ierr)
> >  call PetscViewerFileSetName(viewer, 'sol.vtk', 
ierr);CHKERRA(ierr)

> >
> >  call PetscObjectSetName(u, 'vec_u', ierr);CHKERRA(ierr)
> >  call VecView(u, viewer, ierr);CHKERRA(ierr)
> >
> >  call PetscObjectSetName(v, 'vec_v', ierr);CHKERRA(ierr)
> >  call VecView(v, viewer, ierr);CHKERRA(ierr)
> >
> >  call PetscViewerDestroy(viewer, ierr);CHKERRA(ierr)
> >
> >  call DMRestoreGlobalVector(dm, u, ierr);CHKERRA(ierr)
> >  call DMRestoreGlobalVector(dm, v, ierr);CHKERRA(ierr)
> >
> > Thanks,
> >
> > Danyang
> >
> > 
>
>
>
>
> --
> 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/ 
<https://www.cse.buffalo.edu/%7Eknepley/>

>
> --
> Sent from my Android device with K-9 Mail. Please excuse my brevity.

Re: [petsc-users] object name overwritten in VecView

2018-02-28 Thread Danyang Su

Hi Matt，

Thanks for your suggestion and I will use xmf instead. 

Regards,

Danyang

On February 28, 2018 3:58:08 AM PST, Matthew Knepley <knep...@gmail.com> wrote:
>On Wed, Feb 28, 2018 at 12:39 AM, Smith, Barry F. <bsm...@mcs.anl.gov>
>wrote:
>
>>
>>   Matt,
>>
>>   I have confirmed this is reproducible and a bug. The problem arises
>> because
>>
>> frame #0: 0x00010140625a libpetsc.3.8.dylib`
>> PetscViewerVTKAddField_VTK(viewer=0x7fe66760c750,
>> dm=0x7fe668810820,
>PetscViewerVTKWriteFunction=(libpetsc.3.8.dylib`DMPlexVTKWriteAll
>> at plexvtk.c:633), fieldtype=PETSC_VTK_POINT_FIELD,
>> vec=0x7fe66880ee20) at vtkv.c:140
>> frame #1: 0x000101404e6e libpetsc.3.8.dylib`
>> PetscViewerVTKAddField(viewer=0x7fe66760c750,
>dm=0x7fe668810820,
>> PetscViewerVTKWriteFunction=(libpetsc.3.8.dylib`DMPlexVTKWriteAll at
>> plexvtk.c:633), fieldtype=PETSC_VTK_POINT_FIELD,
>vec=0x7fe66880ee20)
>> at vtkv.c:46
>> frame #2: 0x000101e0b7c3
>libpetsc.3.8.dylib`VecView_Plex_Local(v=0x7fe66880ee20,
>> viewer=0x7fe66760c750) at plex.c:301
>> frame #3: 0x000101e0ead7
>libpetsc.3.8.dylib`VecView_Plex(v=0x7fe66880e820,
>> viewer=0x7fe66760c750) at plex.c:348
>>
>> keeps a linked list of vectors that are to be viewed and the vectors
>are
>> the same Vec because they are obtained with DMGetLocalVector().
>>
>> The safest fix is to have PetscViewerVTKAddField_VTK() do a
>VecDuplicate()
>> on the vector passed in and store that in the linked list instead of
>just
>> storing a pointer to the passed in vector (which might and can be
>> overwritten before all the linked vectors are actually stored).
>>
>
>Danyang,
>
>Barry is right, and the bug can be fixed the way he says. However, this
>points out why VTK is bad format. I think a better choice is
>to use HDF5 and XDMF. For example, in my code now I always use
>
>  DMVIewFromOptions(dm, NULL, "-dm_view");
>
>and then later (perhaps several times)
>
>  VecViewFromOptions(u, NULL, "-u_vec_view")
>  VecViewFromOptions(v, NULL, "-v_vec_view")
>
>and then on the command line
>
>  -dm_view hdf5:test.h5 -u_vec_view hdf5:test.h5::append -v_vec_view
>hdf5:test.h5::append
>
>which produces a file
>
>  test.h5
>
>Then I run
>
>  $PETSC_DIR/bin/petsc_gen_xdmf.py test.h5
>
>which produces another file
>
>  test.xmf
>
>This can be loaded by Paraview for visualization.
>
>  Thanks,
>
> Matt
>
>
>
>>
>>   Barry
>>
>>
>>
>> > On Feb 27, 2018, at 10:44 PM, Danyang Su <danyang...@gmail.com>
>wrote:
>> >
>> > Hi All,
>> >
>> > How to set different object names when using multiple VecView? I
>try to
>> use PetscObjectSetName with multiple output, but the object name is
>> overwritten by the last one.
>> >
>> > As shown below, as well as the enclosed files as example, the
>vector
>> name in sol.vtk is vec_v for both vector u and v.
>> >
>> >  call PetscViewerCreate(PETSC_COMM_WORLD, viewer,
>> ierr);CHKERRA(ierr)
>> >  call PetscViewerSetType(viewer, PETSCVIEWERVTK,
>ierr);CHKERRA(ierr)
>> >  call PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK,
>> ierr);CHKERRA(ierr)
>> >  call PetscViewerFileSetName(viewer, 'sol.vtk',
>ierr);CHKERRA(ierr)
>> >
>> >  call PetscObjectSetName(u, 'vec_u', ierr);CHKERRA(ierr)
>> >  call VecView(u, viewer, ierr);CHKERRA(ierr)
>> >
>> >  call PetscObjectSetName(v, 'vec_v', ierr);CHKERRA(ierr)
>> >  call VecView(v, viewer, ierr);CHKERRA(ierr)
>> >
>> >  call PetscViewerDestroy(viewer, ierr);CHKERRA(ierr)
>> >
>> >  call DMRestoreGlobalVector(dm, u, ierr);CHKERRA(ierr)
>> >  call DMRestoreGlobalVector(dm, v, ierr);CHKERRA(ierr)
>> >
>> > Thanks,
>> >
>> > Danyang
>> >
>> > 
>>
>>
>
>
>-- 
>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.caam.rice.edu/~mk51/>

-- 
Sent from my Android device with K-9 Mail. Please excuse my brevity.

[petsc-users] object name overwritten in VecView

2018-02-27 Thread Danyang Su


Hi All,

How to set different object names when using multiple VecView? I try to 
use PetscObjectSetName with multiple output, but the object name is 
overwritten by the last one.


As shown below, as well as the enclosed files as example, the vector 
name in sol.vtk is vec_v for both vector u and v.


  call PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr);CHKERRA(ierr)
  call PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr);CHKERRA(ierr)
  call PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK, 
ierr);CHKERRA(ierr)

  call PetscViewerFileSetName(viewer, 'sol.vtk', ierr);CHKERRA(ierr)

  call PetscObjectSetName(u, 'vec_u', ierr);CHKERRA(ierr)
  call VecView(u, viewer, ierr);CHKERRA(ierr)

  call PetscObjectSetName(v, 'vec_v', ierr);CHKERRA(ierr)
  call VecView(v, viewer, ierr);CHKERRA(ierr)

  call PetscViewerDestroy(viewer, ierr);CHKERRA(ierr)

  call DMRestoreGlobalVector(dm, u, ierr);CHKERRA(ierr)
  call DMRestoreGlobalVector(dm, v, ierr);CHKERRA(ierr)

Thanks,

Danyang

  program DMPlexTestField
#include "petsc/finclude/petscdmplex.h"
#include "petsc/finclude/petscdmlabel.h"
  use petscdmplex
  implicit none

  DM :: dm
  DMLabel :: label
  Vec :: u, v
  PetscViewer :: viewer
  PetscSection :: section
  PetscInt :: dim,numCells,numFields,numBC
  PetscInt :: i,val
  PetscInt, target, dimension(3) ::  numComp
  PetscInt, pointer :: pNumComp(:)
  PetscInt, target, dimension(12) ::  numDof
  PetscInt, pointer :: pNumDof(:)
  PetscInt, target, dimension(1) ::  bcField
  PetscInt, pointer :: pBcField(:)
  PetscInt :: zero,eight
  IS, target, dimension(1) ::   bcCompIS
  IS, target, dimension(1) ::   bcPointIS
  IS, pointer :: pBcCompIS(:)
  IS, pointer :: pBcPointIS(:)
  PetscBool :: interpolate
  PetscErrorCode :: ierr

  call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
  if (ierr .ne. 0) then
print*,'Unable to initialize PETSc'
stop
  endif
  dim = 2
  call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-dim', dim,PETSC_NULL_BOOL, ierr);CHKERRA(ierr)
  interpolate = PETSC_TRUE
! Create a mesh
  if (dim .eq. 2) then
 numCells = 2
  else
 numCells = 1
  endif
  call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, numCells,interpolate, dm, ierr);CHKERRA(ierr)
! Create a scalar field u, a vector field v, and a surface vector field w
  numFields  = 3
  numComp(1) = 1
  numComp(2) = 1
  numComp(3) = 1
  pNumComp => numComp
  do i = 1, numFields*(dim+1)
 numDof(i) = 0
  end do
! Let u be defined on vertices
  numDof(0*(dim+1)+1) = 1
! Let v be defined on cells
  numDof(1*(dim+1)+1) = 1
! Let v be defined on faces
  numDof(2*(dim+1)+1)   = 1
  pNumDof => numDof
! Setup boundary conditions
  numBC = 1
! Test label retrieval
  call DMGetLabel(dm, 'marker', label, ierr);CHKERRA(ierr)
  zero = 0
  call DMLabelGetValue(label, zero, val, ierr);CHKERRA(ierr)
  if (val .ne. -1) then
CHKERRA(1)
  endif
  eight = 8
  call DMLabelGetValue(label, eight, val, ierr);CHKERRA(ierr)
  if (val .ne. 1) then
CHKERRA(1)
  endif
! Prescribe a Dirichlet condition on u on the boundary
!   Label "marker" is made by the mesh creation routine
  bcField(1) = 0
  pBcField => bcField
  call ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, bcCompIS(1), ierr);CHKERRA(ierr)
  pBcCompIS => bcCompIS
  call DMGetStratumIS(dm, 'marker', 1, bcPointIS(1),ierr);CHKERRA(ierr)
  pBcPointIS => bcPointIS
! Create a PetscSection with this data layout
  call DMPlexCreateSection(dm,dim,numFields,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr)
  CHKERRA(ierr)
  call ISDestroy(bcCompIS(1), ierr);CHKERRA(ierr)
  call ISDestroy(bcPointIS(1), ierr);CHKERRA(ierr)
! Name the Field variables
  call PetscSectionSetFieldName(section, 0, 'u', ierr);CHKERRA(ierr)
  call PetscSectionSetFieldName(section, 1, 'v', ierr);CHKERRA(ierr)
  call PetscSectionSetFieldName(section, 2, 'w', ierr);CHKERRA(ierr)
  call PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr);CHKERRA(ierr)
! Tell the DM to use this data layout
  call DMSetDefaultSection(dm, section, ierr);CHKERRA(ierr)
! Create a Vec with this layout and view it
  call DMGetGlobalVector(dm, u, ierr);CHKERRA(ierr)
  call VecDuplicate(u,v,ierr);CHKERRA(ierr)

  call PetscViewerCreate(PETSC_COMM_WORLD, viewer, ierr);CHKERRA(ierr)
  call PetscViewerSetType(viewer, PETSCVIEWERVTK, ierr);CHKERRA(ierr)
  call PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK, ierr);CHKERRA(ierr)
  call PetscViewerFileSetName(viewer, 'sol.vtk', ierr);CHKERRA(ierr)

  call PetscObjectSetName(u, 'vec_u', ierr);CHKERRA(ierr)
  call VecView(u,

Re: [petsc-users] Cell type for DMPlexCreateFromCellList

2018-02-23 Thread Danyang Su



On 18-02-23 03:04 AM, Matthew Knepley wrote:
On Fri, Feb 23, 2018 at 1:33 AM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi All,

What cell types does DMPlexCreateFromCellList support? I test this
with triangle, tetrahedron and prism. Both triangle and
tetrahedron work but prism mesh throws error saying "Cone size 6
not supported for dimension 3".

Could anyone tell me all the supported cell types?


The limitation occurs in two places:

  1) Calculating edges and faces: I only know how to do this for tri, 
tet, quad, and hex. Give PETSC_FALSE to CreateFromCellList() and this 
error will go away.
Passing PETSC_FALSE to CreateFromCellList() works. There is no problem 
in the distributed nodes and cells. Is there any plan to add 
prism(wedge) in the near future?


      You could also provide the information for interpolating prisms, 
which would depend on how they are ordered when you read in cells.


  2) Even if you read them in, I have no geometric routines for prisms.
That's fine for me. I just need to create dmplex based on the given 
cell/vertex list and distribute over all the processors. Then follow the 
general routines that have been used in the structured grid.


Thanks,

Danyang


  Thanks,

     Matt

Thanks,

Danyang




--
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.caam.rice.edu/%7Emk51/>

[petsc-users] Cell type for DMPlexCreateFromCellList

2018-02-22 Thread Danyang Su


Hi All,

What cell types does DMPlexCreateFromCellList support? I test this with 
triangle, tetrahedron and prism. Both triangle and tetrahedron work but 
prism mesh throws error saying "Cone size 6 not supported for dimension 3".


Could anyone tell me all the supported cell types?

Thanks,

Danyang

Re: [petsc-users] Question on DMPlexCreateSection for Fortran

2018-02-22 Thread Danyang Su


Hi Matt,

Just to let you know that after updating to PETSc 3.8.3 version, the 
DMPlexCreateSection in my code now works.


One more question, what is the PETSC_NULL_XXX for IS pointer, as shown 
below, in C code, it just pass NULL, but in fortran, what is the name of 
null object for pBcCompIS and pBcPointIS.


    call DMPlexCreateSection(dmda_flow%da,dmda_flow%dim,   &
numFields,pNumComp,pNumDof,   &
numBC,pBcField,   &
pBcCompIS,pBcPointIS, &
 PETSC_NULL_IS,section,ierr)
    CHKERRQ(ierr)

Thanks,

Danyang

On 18-02-21 09:22 AM, Danyang Su wrote:


Hi Matt,

To test the Segmentation Violation problem in my code, I modified the 
example ex1f90.F to reproduce the problem I have in my own code.


If use DMPlexCreateBoxMesh to generate the mesh, the code works fine. 
However, if I use DMPlexCreateGmshFromFile, using the same mesh 
exported from "DMPlexCreateBoxMesh", it gives Segmentation Violation 
error.


Did I miss something in the input mesh file? My first guess is the 
label "marker" used in the code, but I couldn't find any place to set 
this label.


Would you please let me know how to solve this problem. My code is 
done in a similar way as ex1f90, it reads mesh from external file or 
creates from cell list, distributes the mesh (these already work), and 
then creates sections and sets ndof to the nodes.


Thanks,

Danyang


On 18-02-20 10:07 AM, Danyang Su wrote:

On 18-02-20 09:52 AM, Matthew Knepley wrote:
On Tue, Feb 20, 2018 at 12:30 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi All,

I tried to compile the DMPlexCreateSection code but got error
information as shown below.

Error: Symbol 'petsc_null_is' at (1) has no IMPLICIT type

I tried to use PETSC_NULL_OBJECT instead of PETSC_NULL_IS, then
the code can be compiled but run into Segmentation Violation
error in DMPlexCreateSection.

From the webpage

http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexCreateSection.html 



The F90 version is DMPlexCreateSectionF90. Doing this with F77 
arrays would have been too painful.

Hi Matt,

Sorry, I still cannot compile the code if use DMPlexCreateSectionF90 
instead of DMPlexCreateSection. Would you please tell me in more 
details?


undefined reference to `dmplexcreatesectionf90_'

then I #include , but this throws 
more error during compilation.



    Included at 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/petscdmplex.h90:6:

    Included at ../../solver/solver_ddmethod.F90:62:

  PETSCSECTION_HIDE section
  1
Error: Unclassifiable statement at (1)
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/ftn-custom/petscdmplex.h90:167.10:
    Included at 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/petscdmplex.h90:6:

    Included at ../../solver/solver_ddmethod.F90:62:

  PETSCSECTION_HIDE section
  1
Error: Unclassifiable statement at (1)
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/ftn-custom/petscdmplex.h90:179.10:
    Included at 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/petscdmplex.h90:6:

    Included at ../../solver/solver_ddmethod.F90:62:



  Thanks,

     Matt

dmda_flow%da is distributed dm object that works fine.

The fortran example I follow is

http://www.mcs.anl.gov/petsc/petsc-dev/src/dm/impls/plex/examples/tutorials/ex1f90.F90

<http://www.mcs.anl.gov/petsc/petsc-dev/src/dm/impls/plex/examples/tutorials/ex1f90.F90>.


What parameters should I use if passing null to bcField,
bcComps, bcPoints and perm.

PetscErrorCode

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscErrorCode.html#PetscErrorCode>
  DMPlexCreateSection

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMPLEX/DMPlexCreateSection.html#DMPlexCreateSection>(DM
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DM/DM.html#DM> 
 dm,PetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  dim,PetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numFields,constPetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numComp[],constPetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numDof[],PetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numBC,constPetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  bcField[],
constIS
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/IS/IS.html#IS> 
 bcComps[], constIS
<https://www.mcs.anl.gov/pe

Re: [petsc-users] Question on DMPlexCreateSection for Fortran

2018-02-21 Thread Danyang Su


Hi Matt,

To test the Segmentation Violation problem in my code, I modified the 
example ex1f90.F to reproduce the problem I have in my own code.


If use DMPlexCreateBoxMesh to generate the mesh, the code works fine. 
However, if I use DMPlexCreateGmshFromFile, using the same mesh exported 
from "DMPlexCreateBoxMesh", it gives Segmentation Violation error.


Did I miss something in the input mesh file? My first guess is the label 
"marker" used in the code, but I couldn't find any place to set this label.


Would you please let me know how to solve this problem. My code is done 
in a similar way as ex1f90, it reads mesh from external file or creates 
from cell list, distributes the mesh (these already work), and then 
creates sections and sets ndof to the nodes.


Thanks,

Danyang


On 18-02-20 10:07 AM, Danyang Su wrote:

On 18-02-20 09:52 AM, Matthew Knepley wrote:
On Tue, Feb 20, 2018 at 12:30 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi All,

I tried to compile the DMPlexCreateSection code but got error
information as shown below.

Error: Symbol 'petsc_null_is' at (1) has no IMPLICIT type

I tried to use PETSC_NULL_OBJECT instead of PETSC_NULL_IS, then
the code can be compiled but run into Segmentation Violation
error in DMPlexCreateSection.

From the webpage

http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexCreateSection.html 



The F90 version is DMPlexCreateSectionF90. Doing this with F77 arrays 
would have been too painful.

Hi Matt,

Sorry, I still cannot compile the code if use DMPlexCreateSectionF90 
instead of DMPlexCreateSection. Would you please tell me in more details?


undefined reference to `dmplexcreatesectionf90_'

then I #include , but this throws more 
error during compilation.



    Included at 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/petscdmplex.h90:6:

    Included at ../../solver/solver_ddmethod.F90:62:

  PETSCSECTION_HIDE section
  1
Error: Unclassifiable statement at (1)
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/ftn-custom/petscdmplex.h90:167.10:
    Included at 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/petscdmplex.h90:6:

    Included at ../../solver/solver_ddmethod.F90:62:

  PETSCSECTION_HIDE section
  1
Error: Unclassifiable statement at (1)
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/ftn-custom/petscdmplex.h90:179.10:
    Included at 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/petscdmplex.h90:6:

    Included at ../../solver/solver_ddmethod.F90:62:



  Thanks,

     Matt

dmda_flow%da is distributed dm object that works fine.

The fortran example I follow is

http://www.mcs.anl.gov/petsc/petsc-dev/src/dm/impls/plex/examples/tutorials/ex1f90.F90

<http://www.mcs.anl.gov/petsc/petsc-dev/src/dm/impls/plex/examples/tutorials/ex1f90.F90>.


What parameters should I use if passing null to bcField, bcComps,
bcPoints and perm.

PetscErrorCode

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscErrorCode.html#PetscErrorCode>
  DMPlexCreateSection

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMPLEX/DMPlexCreateSection.html#DMPlexCreateSection>(DM
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DM/DM.html#DM> 
 dm,PetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  dim,PetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numFields,constPetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numComp[],constPetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numDof[],PetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numBC,constPetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  bcField[],
constIS
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/IS/IS.html#IS> 
 bcComps[], constIS
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/IS/IS.html#IS> 
 bcPoints[],IS
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/IS/IS.html#IS> 
 perm,PetscSection

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/IS/PetscSection.html#PetscSection>
  *section)

#include 
#include 
#include 

...

#ifdef USG
    numFields = 1
    numComp(1) = 1
    pNumComp => numComp

    do i = 1, numFields*(dmda_flow%dim+1)
  numDof(i) = 0
    end do
    numDof(0*(dmda_flow%dim+1)+1) = dmda_flow%dof
    pNumDof => numDof

    numBC = 0

Re: [petsc-users] Question on DMPlexCreateSection for Fortran

2018-02-20 Thread Danyang Su


On 18-02-20 09:52 AM, Matthew Knepley wrote:
On Tue, Feb 20, 2018 at 12:30 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi All,

I tried to compile the DMPlexCreateSection code but got error
information as shown below.

Error: Symbol 'petsc_null_is' at (1) has no IMPLICIT type

I tried to use PETSC_NULL_OBJECT instead of PETSC_NULL_IS, then
the code can be compiled but run into Segmentation Violation error
in DMPlexCreateSection.

From the webpage

http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexCreateSection.html 



The F90 version is DMPlexCreateSectionF90. Doing this with F77 arrays 
would have been too painful.

Hi Matt,

Sorry, I still cannot compile the code if use DMPlexCreateSectionF90 
instead of DMPlexCreateSection. Would you please tell me in more details?


undefined reference to `dmplexcreatesectionf90_'

then I #include , but this throws more 
error during compilation.



    Included at 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/petscdmplex.h90:6:

    Included at ../../solver/solver_ddmethod.F90:62:

  PETSCSECTION_HIDE section
  1
Error: Unclassifiable statement at (1)
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/ftn-custom/petscdmplex.h90:167.10:
    Included at 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/petscdmplex.h90:6:

    Included at ../../solver/solver_ddmethod.F90:62:

  PETSCSECTION_HIDE section
  1
Error: Unclassifiable statement at (1)
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/ftn-custom/petscdmplex.h90:179.10:
    Included at 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/finclude/petscdmplex.h90:6:

    Included at ../../solver/solver_ddmethod.F90:62:



  Thanks,

     Matt

dmda_flow%da is distributed dm object that works fine.

The fortran example I follow is

http://www.mcs.anl.gov/petsc/petsc-dev/src/dm/impls/plex/examples/tutorials/ex1f90.F90

<http://www.mcs.anl.gov/petsc/petsc-dev/src/dm/impls/plex/examples/tutorials/ex1f90.F90>.


What parameters should I use if passing null to bcField, bcComps,
bcPoints and perm.

PetscErrorCode

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscErrorCode.html#PetscErrorCode>
  DMPlexCreateSection

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMPLEX/DMPlexCreateSection.html#DMPlexCreateSection>(DM
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DM/DM.html#DM> 
 dm,PetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  dim,PetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numFields,constPetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numComp[],constPetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numDof[],PetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  numBC,constPetscInt

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/PetscInt.html#PetscInt>
  bcField[],
constIS
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/IS/IS.html#IS> 
 bcComps[], constIS
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/IS/IS.html#IS> 
 bcPoints[],IS
<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/IS/IS.html#IS> 
 perm,PetscSection

<https://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/IS/PetscSection.html#PetscSection>
  *section)

#include 
#include 
#include 

...

#ifdef USG
    numFields = 1
    numComp(1) = 1
    pNumComp => numComp

    do i = 1, numFields*(dmda_flow%dim+1)
  numDof(i) = 0
    end do
    numDof(0*(dmda_flow%dim+1)+1) = dmda_flow%dof
    pNumDof => numDof

    numBC = 0

    call DMPlexCreateSection(dmda_flow%da,dmda_flow%dim, &
numFields,pNumComp,pNumDof, &
numBC,PETSC_NULL_INTEGER, &
PETSC_NULL_IS,PETSC_NULL_IS, & !Error here
PETSC_NULL_IS,section,ierr)
    CHKERRQ(ierr)

    call PetscSectionSetFieldName(section,0,'flow',ierr)
    CHKERRQ(ierr)

    call DMSetDefaultSection(dmda_flow%da,section,ierr)
    CHKERRQ(ierr)

    call PetscSectionDestroy(section,ierr)
    CHKERRQ(ierr)
#endif

Thanks,

Danyang




--
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.caam.rice.edu/%7Emk51/>

[petsc-users] Question on DMPlexCreateSection for Fortran

2018-02-20 Thread Danyang Su


Hi All,

I tried to compile the DMPlexCreateSection code but got error 
information as shown below.


Error: Symbol 'petsc_null_is' at (1) has no IMPLICIT type

I tried to use PETSC_NULL_OBJECT instead of PETSC_NULL_IS, then the code 
can be compiled but run into Segmentation Violation error in 
DMPlexCreateSection.


dmda_flow%da is distributed dm object that works fine.

The fortran example I follow is 
http://www.mcs.anl.gov/petsc/petsc-dev/src/dm/impls/plex/examples/tutorials/ex1f90.F90. 



What parameters should I use if passing null to bcField, bcComps, 
bcPoints and perm.


PetscErrorCode 
  DMPlexCreateSection 
(DM 
  dm,PetscInt 
  dim,PetscInt 
  numFields,constPetscInt 
  numComp[],constPetscInt 
  numDof[],PetscInt 
  numBC,constPetscInt 
  bcField[],
constIS 
  bcComps[], constIS 
  bcPoints[],IS 
  perm,PetscSection 
  *section)


#include 
#include 
#include 

...

#ifdef USG
    numFields = 1
    numComp(1) = 1
    pNumComp => numComp

    do i = 1, numFields*(dmda_flow%dim+1)
  numDof(i) = 0
    end do
    numDof(0*(dmda_flow%dim+1)+1) = dmda_flow%dof
    pNumDof => numDof

    numBC = 0

    call DMPlexCreateSection(dmda_flow%da,dmda_flow%dim,  &
numFields,pNumComp,pNumDof,   &
numBC,PETSC_NULL_INTEGER,   &
PETSC_NULL_IS,PETSC_NULL_IS, & !Error here
 PETSC_NULL_IS,section,ierr)
    CHKERRQ(ierr)

    call PetscSectionSetFieldName(section,0,'flow',ierr)
    CHKERRQ(ierr)

    call DMSetDefaultSection(dmda_flow%da,section,ierr)
    CHKERRQ(ierr)

    call PetscSectionDestroy(section,ierr)
    CHKERRQ(ierr)
#endif

Thanks,

Danyang

[petsc-users] how to check if cell is local owned in DMPlex

2018-02-19 Thread Danyang Su


Hi Matt,

Would you please let me know how to check if a cell is local owned? When 
overlap is 0 in DMPlexDistribute, all the cells are local owned. How 
about overlap > 0? It sounds like impossible to check by node because a 
cell can be local owned even if none of the nodes in this cell is local 
owned.


Thanks,

Danyang

Re: [petsc-users] Error when use DMPlexGetVertexNumbering


On 18-02-16 10:50 AM, Matthew Knepley wrote:
On Fri, Feb 16, 2018 at 1:45 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi Matt,

I try to get the global vertex index and cell index from local
mesh and run into problem. What I need is local to global index
(the original index used in DMPlexCreateFromCellList is best, as
user know exactly where the node/cell is) for vertices and cells,
which will be used to assign material properties and some
parameters to the specified cell/vertex.


I would recommend doing this before you distribute the mesh. Just set 
these properties using a DMLabel and it will be automatically distributed.

I will try this. Thanks.


I can use coordinates to select vertex/cell which has already
included, but still want to keep this feature. This is pretty
straightforward when using structured grid. For the unstructured
grid, I just got compiling error saying "You need a ISO C
conforming compiler to use the glibc headers"


For any compilation problem, you have to send the configure.log and 
make.log. However, it appears that you are not using the same compiler 
that you configured with.
Sorry for confusion. There are linux-gnu-dbg (debug version) and 
linux-gnu-opt (optimized version) configuration, the one I used to 
compile is debug version and the attached one is optimized version. I 
will try you recommendation first.


Thanks,

Danyang


   Matt

Would you please let me know if I need to change the configuration
of PETSc or is there any alternative ways to avoid using
DMPlexGetVertexNumbering and DMPlexGetCellNumbering but get local
to global index?

The error information during compilation is shown below, followed
by PETSc configuration.

 -o ../../solver/solver_ddmethod.o ../../solver/solver_ddmethod.F90

In file included from /usr/include/features.h:375:0,
 from /usr/include/stdio.h:28,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscsys.h:161,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscis.h:8,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscvec.h:10,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscmat.h:7,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/private/dmpleximpl.h:5,
 from ../../solver/solver_ddmethod.F90:4122:
/usr/include/x86_64-linux-gnu/sys/cdefs.h:30:3: error: #error "You
need a ISO C conforming compiler to us\
e the glibc headers"
 # error "You need a ISO C conforming compiler to use the glibc
headers"
   ^
In file included from /usr/include/features.h:399:0,
 from /usr/include/stdio.h:28,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscsys.h:161,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscis.h:8,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscvec.h:10,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscmat.h:7,
 from
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/private/dmpleximpl.h:5,
 from ../../solver/solver_ddmethod.F90:4122:
/usr/include/x86_64-linux-gnu/gnu/stubs.h:7:0: fatal error:
gnu/stubs-32.h: No such file or directory
 # include 
 ^
compilation terminated.
make: [../../solver/solver_ddmethod.o] Error 1 (ignored)


PETSc configuration
--with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-mumps
--download-scalapack --download-parmetis --download-metis
--download-ptscotch --download-fblaslapack --download-mpich
--download-hypre --download-superlu_dist --download-hdf5=yes
--with-debugging=0 COPTFLAGS="-O3 -march=native -mtune=native"
CXXOPTFLAGS="-O3 -march=native -mtune=native" FOPTFLAGS="-O3
-march=native -mtune=native"

Thanks and regards,

Danyang




--
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.caam.rice.edu/%7Emk51/>

[petsc-users] Error when use DMPlexGetVertexNumbering


Hi Matt,

I try to get the global vertex index and cell index from local mesh and 
run into problem. What I need is local to global index (the original 
index used in DMPlexCreateFromCellList is best, as user know exactly 
where the node/cell is) for vertices and cells, which will be used to 
assign material properties and some parameters to the specified 
cell/vertex. I can use coordinates to select vertex/cell which has 
already included, but still want to keep this feature. This is pretty 
straightforward when using structured grid. For the unstructured grid, I 
just got compiling error saying "You need a ISO C conforming compiler to 
use the glibc headers"


Would you please let me know if I need to change the configuration of 
PETSc or is there any alternative ways to avoid using 
DMPlexGetVertexNumbering and DMPlexGetCellNumbering but get local to 
global index?


The error information during compilation is shown below, followed by 
PETSc configuration.


 -o ../../solver/solver_ddmethod.o ../../solver/solver_ddmethod.F90

In file included from /usr/include/features.h:375:0,
 from /usr/include/stdio.h:28,
 from 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscsys.h:161,

 from /home/dsu/Soft/PETSc/petsc-3.7.5/include/petscis.h:8,
 from 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscvec.h:10,
 from 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscmat.h:7,
 from 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/private/dmpleximpl.h:5,

 from ../../solver/solver_ddmethod.F90:4122:
/usr/include/x86_64-linux-gnu/sys/cdefs.h:30:3: error: #error "You need 
a ISO C conforming compiler to us\

e the glibc headers"
 # error "You need a ISO C conforming compiler to use the glibc headers"
   ^
In file included from /usr/include/features.h:399:0,
 from /usr/include/stdio.h:28,
 from 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscsys.h:161,

 from /home/dsu/Soft/PETSc/petsc-3.7.5/include/petscis.h:8,
 from 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscvec.h:10,
 from 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petscmat.h:7,
 from 
/home/dsu/Soft/PETSc/petsc-3.7.5/include/petsc/private/dmpleximpl.h:5,

 from ../../solver/solver_ddmethod.F90:4122:
/usr/include/x86_64-linux-gnu/gnu/stubs.h:7:0: fatal error: 
gnu/stubs-32.h: No such file or directory

 # include 
 ^
compilation terminated.
make: [../../solver/solver_ddmethod.o] Error 1 (ignored)


PETSc configuration
--with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-mumps 
--download-scalapack --download-parmetis --download-metis 
--download-ptscotch --download-fblaslapack --download-mpich 
--download-hypre --download-superlu_dist --download-hdf5=yes 
--with-debugging=0 COPTFLAGS="-O3 -march=native -mtune=native" 
CXXOPTFLAGS="-O3 -march=native -mtune=native" FOPTFLAGS="-O3 
-march=native -mtune=native"


Thanks and regards,

Danyang

Re: [petsc-users] Question on DMPlexCreateFromCellList and DMPlexCreateFromFile


On 18-02-16 10:13 AM, Matthew Knepley wrote:
On Fri, Feb 16, 2018 at 11:36 AM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


On 18-02-15 05:57 PM, Matthew Knepley wrote:


On Thu, Feb 15, 2018 at 7:40 PM, Danyang Su <danyang...@gmail.com
<mailto:danyang...@gmail.com>> wrote:

Hi Matt,

I have a question on DMPlexCreateFromCellList and
DMPlexCreateFromFile. When use DMPlexCreateFromFile with Gmsh
file input, it works fine and each processor gets its own
part. However, when use DMPlexCreateFromCellList, all the
processors have the same global mesh. To my understand, I
should put the global mesh as input, right?


No. Each process should get part of the mesh in
CreateFromCellList(), but the most common thing to do is to
feed the whole mesh in on proc 0, and nothing in on the other procs.

Thanks for the explanation. It works now.


Great. Also feel free to suggest improvements, examples, or better 
documentation.
Thanks, I will bother you a lot recently while porting the code from 
structured grid version to unstructured grid version. Thanks in advance.

Danyang


  Thanks,

     Matt

Danyang


  Thanks,

    Matt

Otherwise, I should use DMPlexCreateFromCellListParallel
instead if the input is local mesh.

Below is the test code I use, results from method 1 is wrong
and that from method 2 is correct. Would you please help to
check if I did anything wrong with DMPlexCreateFromCellList
input?

!test with 4 processor, global num_cells = 8268, global
num_nodes = 4250

!correct results

 check rank    2  istart 2034 iend 3116
 check rank    3  istart 2148 iend 3293
 check rank    1  istart 2044 iend 3133
 check rank    0  istart 2042 iend 3131

!wrong results

  check rank    0  istart 8268  iend    12518
  check rank    1  istart 8268  iend    12518
  check rank    2  istart 8268  iend    12518
  check rank    3  istart 8268  iend    12518


  !c *    test part *
  !c method 1: create DMPlex from cell list, same
duplicated global meshes over all processors
  !c the input parameters num_cells, num_nodes,
dmplex_cells, dmplex_verts are all global parameters (global
mesh data)
  call
DMPlexCreateFromCellList(Petsc_Comm_World,ndim,num_cells, &
num_nodes,num_nodes_per_cell,  &
Petsc_True,dmplex_cells,ndim,  &
dmplex_verts,dmda_flow%da,ierr)
  CHKERRQ(ierr)


  !c method 2: create DMPlex from Gmsh file, for test
purpose, this works fine, each processor gets its own part
  call DMPlexCreateFromFile(Petsc_Comm_World, &
prefix(:l_prfx)//'.msh',0, &
dmda_flow%da,ierr)
  CHKERRQ(ierr)

  !c *end of test part*


  distributedMesh = PETSC_NULL_OBJECT

  !c distribute mesh over processes
  call DMPlexDistribute(dmda_flow%da,0,PETSC_NULL_OBJECT, &
distributedMesh,ierr)
  CHKERRQ(ierr)

  !c destroy original global mesh after distribution
  if (distributedMesh /= PETSC_NULL_OBJECT) then
    call DMDestroy(dmda_flow%da,ierr)
    CHKERRQ(ierr)
    !c set the global mesh as distributed mesh
    dmda_flow%da = distributedMesh
  end if

  !c get coordinates
  call DMGetCoordinatesLocal(dmda_flow%da,gc,ierr)
  CHKERRQ(ierr)

  call DMGetCoordinateDM(dmda_flow%da,cda,ierr)
  CHKERRQ(ierr)

  call DMGetDefaultSection(cda,cs,ierr)
  CHKERRQ(ierr)

  call PetscSectionGetChart(cs,istart,iend,ierr)
  CHKERRQ(ierr)

#ifdef DEBUG
    if(info_debug > 0) then
  write(*,*) "check rank ",rank," istart ",istart,"
iend ",iend
    end if
#endif


Thanks and regards,

Danyang




-- 
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.caam.rice.edu/%7Emk51/>





--
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.caam.rice.edu/%7Emk51/>

Re: [petsc-users] Question on DMPlexCreateFromCellList and DMPlexCreateFromFile




On 18-02-15 05:57 PM, Matthew Knepley wrote:
On Thu, Feb 15, 2018 at 7:40 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi Matt,

I have a question on DMPlexCreateFromCellList and
DMPlexCreateFromFile. When use DMPlexCreateFromFile with Gmsh file
input, it works fine and each processor gets its own part.
However, when use DMPlexCreateFromCellList, all the processors
have the same global mesh. To my understand, I should put the
global mesh as input, right?


No. Each process should get part of the mesh in CreateFromCellList(), 
but the most common thing to do is to

feed the whole mesh in on proc 0, and nothing in on the other procs.

Thanks for the explanation. It works now.

Danyang


  Thanks,

    Matt

Otherwise, I should use DMPlexCreateFromCellListParallel instead
if the input is local mesh.

Below is the test code I use, results from method 1 is wrong and
that from method 2 is correct. Would you please help to check if I
did anything wrong with DMPlexCreateFromCellList input?

!test with 4 processor, global num_cells = 8268, global num_nodes
= 4250

!correct results

 check rank    2  istart 2034 iend 3116
 check rank    3  istart 2148 iend 3293
 check rank    1  istart 2044 iend 3133
 check rank    0  istart 2042 iend 3131

!wrong results

  check rank    0  istart 8268 iend    12518
  check rank    1  istart 8268 iend    12518
  check rank    2  istart 8268 iend    12518
  check rank    3  istart 8268 iend    12518


  !c *    test part *
  !c method 1: create DMPlex from cell list, same duplicated
global meshes over all processors
  !c the input parameters num_cells, num_nodes, dmplex_cells,
dmplex_verts are all global parameters (global mesh data)
  call DMPlexCreateFromCellList(Petsc_Comm_World,ndim,num_cells, &
num_nodes,num_nodes_per_cell,  &
Petsc_True,dmplex_cells,ndim,  &
dmplex_verts,dmda_flow%da,ierr)
  CHKERRQ(ierr)


  !c method 2: create DMPlex from Gmsh file, for test purpose,
this works fine, each processor gets its own part
  call DMPlexCreateFromFile(Petsc_Comm_World, &
prefix(:l_prfx)//'.msh',0,  &
  dmda_flow%da,ierr)
  CHKERRQ(ierr)

  !c *end of test part*


  distributedMesh = PETSC_NULL_OBJECT

  !c distribute mesh over processes
  call DMPlexDistribute(dmda_flow%da,0,PETSC_NULL_OBJECT, &
    distributedMesh,ierr)
  CHKERRQ(ierr)

  !c destroy original global mesh after distribution
  if (distributedMesh /= PETSC_NULL_OBJECT) then
    call DMDestroy(dmda_flow%da,ierr)
    CHKERRQ(ierr)
    !c set the global mesh as distributed mesh
    dmda_flow%da = distributedMesh
  end if

  !c get coordinates
  call DMGetCoordinatesLocal(dmda_flow%da,gc,ierr)
  CHKERRQ(ierr)

  call DMGetCoordinateDM(dmda_flow%da,cda,ierr)
  CHKERRQ(ierr)

  call DMGetDefaultSection(cda,cs,ierr)
  CHKERRQ(ierr)

  call PetscSectionGetChart(cs,istart,iend,ierr)
  CHKERRQ(ierr)

#ifdef DEBUG
    if(info_debug > 0) then
  write(*,*) "check rank ",rank," istart ",istart," iend
",iend
    end if
#endif


Thanks and regards,

Danyang




--
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.caam.rice.edu/%7Emk51/>

Re: [petsc-users] Is OpenMP still available for PETSc?

2017-06-30 Thread Danyang Su


Hi Barry,

Thanks for the quick response. What I want to test is to check if OpenMP 
has any benefit when total degrees of freedoms per processor drops below 
5k. When using pure MPI my code shows good speedup if total degrees of 
freedoms per processor is above 10k. But below this value, the parallel 
efficiency decreases.


The petsc 3.6 change log indicates

 * Removed all threadcomm support including --with-pthreadclasses and
   --with-openmpclasses configure arguments

I guess petsc 3.5 version is the last version I can test, right?

Thanks,

Danyang


On 17-06-30 03:49 PM, Barry Smith wrote:

   The current version of PETSc does not use OpenMP, you are free to use OpenMP 
in your portions of the code of course. If you want PETSc using OpenMP you have 
to use the old, unsupported version of PETSc. We never found any benefit to 
using OpenMP.

Barry


On Jun 30, 2017, at 5:40 PM, Danyang Su <danyang...@gmail.com> wrote:

Dear All,

I recalled there was OpenMP available for PETSc for the old development version. When google 
"petsc hybrid mpi openmp", there returned some papers about this feature. My code was 
first parallelized using OpenMP and then redeveloped using PETSc, with OpenMP kept but not used 
together with MPI. Before retesting the code using hybrid mpi-openmp, I picked one PETSc example 
ex10 by adding "omp_set_num_threads(max_threads);" under PetscInitialize.

The PETSc is the current development version configured as follows

--with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-debugging=0 --CFLAGS=-fopenmp --CXXFLAGS=-fopenmp 
--FFLAGS=-fopenmp COPTFLAGS="-O3 -march=native -mtune=native" CXXOPTFLAGS="-O3 -march=native 
-mtune=native" FOPTFLAGS="-O3 -march=native -mtune=native" --with-large-file-io=1 
--download-cmake=yes --download-mumps --download-scalapack --download-parmetis --download-metis 
--download-ptscotch --download-fblaslapack --download-mpich --download-hypre --download-superlu_dist 
--download-hdf5=yes --with-openmp --with-threadcomm --with-pthreadclasses --with-openmpclasses

The code can be successfully compiled. However, when I run the code with 
OpenMP, it does not work, the time shows no change in performance if 1 or 2 
threads per processor is used. Also, the CPU/Threads usage indicates that no 
thread is used.

I just wonder if OpenMP is still available in the latest version, though it is 
not recommended to use.

mpiexec -n 2 ./ex10 -f0 mat_rhs_pc_nonzero/a_react_in_2.bin -rhs 
mat_rhs_pc_nonzero/b_react_in_2.bin -ksp_rtol 1.0e-20 -ksp_monitor 
-ksp_error_if_not_converged -sub_pc_factor_shift_type nonzero -mat_view 
ascii::ascii_info -log_view -max_threads 1 -threadcomm_type openmp 
-threadcomm_nthreads 1

KSPSolve   1 1.0 8.9934e-01 1.0 1.03e+09 1.0 7.8e+01 3.6e+04 
7.8e+01 69 97 89  6 76  89 97 98 98 96  2290
PCSetUp2 1.0 8.9590e-02 1.0 2.91e+07 1.0 0.0e+00 0.0e+00 
0.0e+00  7  3  0  0  0   9  3  0  0  0   648
PCSetUpOnBlocks2 1.0 8.9465e-02 1.0 2.91e+07 1.0 0.0e+00 0.0e+00 
0.0e+00  7  3  0  0  0   9  3  0  0  0   649
PCApply   40 1.0 3.1993e-01 1.0 2.70e+08 1.0 0.0e+00 0.0e+00 
0.0e+00 24 25  0  0  0  32 25  0  0  0  1686

mpiexec -n 2 ./ex10 -f0 mat_rhs_pc_nonzero/a_react_in_2.bin -rhs 
mat_rhs_pc_nonzero/b_react_in_2.bin -ksp_rtol 1.0e-20 -ksp_monitor 
-ksp_error_if_not_converged -sub_pc_factor_shift_type nonzero -mat_view 
ascii::ascii_info -log_view -max_threads 2 -threadcomm_type openmp 
-threadcomm_nthreads 2

KSPSolve   1 1.0 8.9701e-01 1.0 1.03e+09 1.0 7.8e+01 3.6e+04 
7.8e+01 69 97 89  6 76  89 97 98 98 96  2296
PCSetUp2 1.0 8.7635e-02 1.0 2.91e+07 1.0 0.0e+00 0.0e+00 
0.0e+00  7  3  0  0  0   9  3  0  0  0   663
PCSetUpOnBlocks2 1.0 8.7511e-02 1.0 2.91e+07 1.0 0.0e+00 0.0e+00 
0.0e+00  7  3  0  0  0   9  3  0  0  0   664
PCApply   40 1.0 3.1878e-01 1.0 2.70e+08 1.0 0.0e+00 0.0e+00 
0.0e+00 24 25  0  0  0  32 25  0  0  0  1692

Thanks and regards,

Danyang

Re: [petsc-users] PCFactorSetShiftType does not work in code but -pc_factor_set_shift_type works

2017-05-25 Thread Danyang Su


Hi Hong,

It works like a charm. I really appreciate your help.

Regards,

Danyang


On 17-05-25 07:49 AM, Hong wrote:

Danyang:
You must access inner pc, then set shift. See
petsc/src/ksp/ksp/examples/tutorials/ex7.c

For example, I add following to 
petsc/src/ksp/ksp/examples/tutorials/ex2.c, line 191:

  PetscBool isbjacobi;
  PCpc;
  ierr = KSPGetPC(ksp,);CHKERRQ(ierr);
  ierr = 
PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,);CHKERRQ(ierr);

  if (isbjacobi) {
PetscInt nlocal;
KSP  *subksp;
PC   subpc;

ierr = KSPSetUp(ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,);CHKERRQ(ierr);

/* Extract the array of KSP contexts for the local blocks */
ierr = PCBJacobiGetSubKSP(pc,,NULL,);CHKERRQ(ierr);
printf("isbjacobi, nlocal %D, set option to subpc...\n",nlocal);
for (i=0; i<nlocal; i++) {
  ierr = KSPGetPC(subksp[i],);CHKERRQ(ierr);
  ierr = PCFactorSetShiftType(subpc,MAT_SHIFT_NONZERO);CHKERRQ(ierr);
}
  }


Dear Hong and Barry,

I have implemented this option in the code, as we also need to use
configuration from file for convenience. When I run the code using
options, it works fine, however, when I run the code using
configuration file, it does not work. The code has two set of
equations, flow and reactive, with prefix been set to "flow_" and
"react_". When I run the code using

mpiexec -n 4 ../executable -flow_sub_pc_factor_shift_type nonzero
-react_sub_pc_factor_shift_type nonzero

it works. However, if I run using

mpiexec -n 4 ../executable

and let the executable file read the options from file, it just
does not work at "call
PCFactorSetShiftType(pc_flow,MAT_SHIFT_NONZERO, ierr)  or none,
positive_definite ...". Do I miss something here?

Below is the pseudo code I have used for flow equations, similar
for reactive equations.

  call MatCreateAIJ(Petsc_Comm_World,nndof,nndof,nngbldof, &
nngbldof,d_nz,PETSC_NULL_INTEGER,o_nz, &
PETSC_NULL_INTEGER,a_flow,ierr)
  CHKERRQ(ierr)

call MatSetFromOptions(a_flow,ierr)
CHKERRQ(ierr)

call KSPCreate(Petsc_Comm_World, ksp_flow, ierr)
CHKERRQ(ierr)

call KSPAppendOptionsPrefix(ksp_flow,"flow_",ierr)
CHKERRQ(ierr)

call KSPSetInitialGuessNonzero(ksp_flow, &
b_initial_guess_nonzero_flow, ierr)
CHKERRQ(ierr)

call KSPSetInitialGuessNonzero(ksp_flow, &
b_initial_guess_nonzero_flow, ierr)
CHKERRQ(ierr)

call KSPSetDM(ksp_flow,dmda_flow%da,ierr)
CHKERRQ(ierr)
call KSPSetDMActive(ksp_flow,PETSC_FALSE,ierr)
CHKERRQ(ierr)

*CHECK IF READ OPTION FROM FILE*
if (read_option_from_file) then

  call KSPSetType(ksp_flow, KSPGMRES, ierr) !or KSPBCGS or
others...
  CHKERRQ(ierr)

  call KSPGetPC(ksp_flow, pc_flow, ierr)
  CHKERRQ(ierr)

  call PCSetType(pc_flow,PCBJACOBI, ierr)   !or PCILU
or PCJACOBI or PCHYPRE ...
  CHKERRQ(ierr)

  call PCFactorSetShiftType(pc_flow,MAT_SHIFT_NONZERO,
ierr)  or none, positive_definite ...
  CHKERRQ(ierr)

end if

call
PCFactorGetMatSolverPackage(pc_flow,solver_pkg_flow,ierr)
CHKERRQ(ierr)

call compute_jacobian(rank,dmda_flow%da, &
a_flow,a_in,ia_in,ja_in,nngl_in, &
row_idx_l2pg,col_idx_l2pg,&
  b_non_interlaced)
call KSPSetFromOptions(ksp_flow,ierr)
CHKERRQ(ierr)

call KSPSetUp(ksp_flow,ierr)
CHKERRQ(ierr)

call KSPSetUpOnBlocks(ksp_flow,ierr)
CHKERRQ(ierr)

call KSPSolve(ksp_flow,b_flow,x_flow,ierr)
CHKERRQ(ierr)


Thanks and Regards,

Danyang

On 17-05-24 06:32 PM, Hong wrote:

Remove your option '-vecload_block_size 10'.
Hong

On Wed, May 24, 2017 at 3:06 PM, Danyang Su <danyang...@gmail.com
<mailto:danyang...@gmail.com>> wrote:

Dear Hong,

I just tested with different number of processors for the
same matrix. It sometimes got "ERROR: Arguments are
incompatible" for different number of processors. It works
fine using 4, 8, or 24 processors, but failed with "ERROR:
Arguments are incompatible" using 16 or 48 processors. The
error information is attached. I tested this on my local
computer with 6 cores 12 threads. Any suggestion on this?

Thanks,

Danyang


On 17-05-24 12:28 PM, Danyang Su wrote:


Hi Hong,

[petsc-users] PCFactorSetShiftType does not work in code but -pc_factor_set_shift_type works

2017-05-25 Thread Danyang Su


Dear Hong and Barry,

I have implemented this option in the code, as we also need to use 
configuration from file for convenience. When I run the code using 
options, it works fine, however, when I run the code using configuration 
file, it does not work. The code has two set of equations, flow and 
reactive, with prefix been set to "flow_" and "react_". When I run the 
code using


mpiexec -n 4 ../executable -flow_sub_pc_factor_shift_type nonzero 
-react_sub_pc_factor_shift_type nonzero


it works. However, if I run using

mpiexec -n 4 ../executable

and let the executable file read the options from file, it just does not 
work at "call PCFactorSetShiftType(pc_flow,MAT_SHIFT_NONZERO, ierr)  or 
none, positive_definite ...". Do I miss something here?


Below is the pseudo code I have used for flow equations, similar for 
reactive equations.


  call MatCreateAIJ(Petsc_Comm_World,nndof,nndof,nngbldof, &
nngbldof,d_nz,PETSC_NULL_INTEGER,o_nz, &
PETSC_NULL_INTEGER,a_flow,ierr)
  CHKERRQ(ierr)

call MatSetFromOptions(a_flow,ierr)
CHKERRQ(ierr)

call KSPCreate(Petsc_Comm_World, ksp_flow, ierr)
CHKERRQ(ierr)

call KSPAppendOptionsPrefix(ksp_flow,"flow_",ierr)
CHKERRQ(ierr)

call KSPSetInitialGuessNonzero(ksp_flow,   &
b_initial_guess_nonzero_flow, ierr)
CHKERRQ(ierr)

call KSPSetInitialGuessNonzero(ksp_flow,   &
b_initial_guess_nonzero_flow, ierr)
CHKERRQ(ierr)

call KSPSetDM(ksp_flow,dmda_flow%da,ierr)
CHKERRQ(ierr)
call KSPSetDMActive(ksp_flow,PETSC_FALSE,ierr)
CHKERRQ(ierr)

*CHECK IF READ OPTION FROM FILE*
if (read_option_from_file) then

  call KSPSetType(ksp_flow, KSPGMRES, ierr)!or KSPBCGS or 
others...

  CHKERRQ(ierr)

  call KSPGetPC(ksp_flow, pc_flow, ierr)
  CHKERRQ(ierr)

  call PCSetType(pc_flow,PCBJACOBI, ierr)   !or PCILU or 
PCJACOBI or PCHYPRE ...

  CHKERRQ(ierr)

  call PCFactorSetShiftType(pc_flow,MAT_SHIFT_NONZERO, ierr)  
or none, positive_definite ...

  CHKERRQ(ierr)

end if

call PCFactorGetMatSolverPackage(pc_flow,solver_pkg_flow,ierr)
CHKERRQ(ierr)

call compute_jacobian(rank,dmda_flow%da,   &
a_flow,a_in,ia_in,ja_in,nngl_in, &
row_idx_l2pg,col_idx_l2pg,   &
  b_non_interlaced)
call KSPSetFromOptions(ksp_flow,ierr)
CHKERRQ(ierr)

call KSPSetUp(ksp_flow,ierr)
CHKERRQ(ierr)

call KSPSetUpOnBlocks(ksp_flow,ierr)
CHKERRQ(ierr)

call KSPSolve(ksp_flow,b_flow,x_flow,ierr)
CHKERRQ(ierr)


Thanks and Regards,

Danyang

On 17-05-24 06:32 PM, Hong wrote:

Remove your option '-vecload_block_size 10'.
Hong

On Wed, May 24, 2017 at 3:06 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Dear Hong,

I just tested with different number of processors for the same
matrix. It sometimes got "ERROR: Arguments are incompatible" for
different number of processors. It works fine using 4, 8, or 24
processors, but failed with "ERROR: Arguments are incompatible"
using 16 or 48 processors. The error information is attached. I
tested this on my local computer with 6 cores 12 threads. Any
suggestion on this?

Thanks,

Danyang


On 17-05-24 12:28 PM, Danyang Su wrote:


Hi Hong,

Awesome. Thanks for testing the case. I will try your options for
the code and get back to you later.

Regards,

Danyang


On 17-05-24 12:21 PM, Hong wrote:

Danyang :
I tested your data.
Your matrices encountered zero pivots, e.g.
petsc/src/ksp/ksp/examples/tutorials (master)
$ mpiexec -n 24 ./ex10 -f0 a_react_in_2.bin -rhs
b_react_in_2.bin -ksp_monitor -ksp_error_if_not_converged

[15]PETSC ERROR: Zero pivot in LU factorization:
http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
<http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot>
[15]PETSC ERROR: Zero pivot row 1249 value 2.05808e-14 tolerance
2.22045e-14
...

Adding option '-sub_pc_factor_shift_type nonzero', I got
mpiexec -n 24 ./ex10 -f0 a_react_in_2.bin -rhs b_react_in_2.bin
-ksp_monitor -ksp_error_if_not_converged
-sub_pc_factor_shift_type nonzero -mat_view ascii::ascii_info

Mat Object: 24 MPI processes
  type: mpiaij
  rows=45, cols=45
  total: nonzeros=6991400, allocated nonzeros=6991400
  total number of mallocs used during MatSetValues calls =0
not using I-node (on process 0) routines
  0 KSP Residual norm 5.84911755e+01
  1 KSP Residual norm

Re: [petsc-users] Question on incomplete factorization level and fill


Hi All,

I just delete the .info file and it works without problem now.

Thanks,

Danyang


On 17-05-24 06:32 PM, Hong wrote:

Remove your option '-vecload_block_size 10'.
Hong

On Wed, May 24, 2017 at 3:06 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Dear Hong,

I just tested with different number of processors for the same
matrix. It sometimes got "ERROR: Arguments are incompatible" for
different number of processors. It works fine using 4, 8, or 24
processors, but failed with "ERROR: Arguments are incompatible"
using 16 or 48 processors. The error information is attached. I
tested this on my local computer with 6 cores 12 threads. Any
suggestion on this?

Thanks,

Danyang


    On 17-05-24 12:28 PM, Danyang Su wrote:


Hi Hong,

Awesome. Thanks for testing the case. I will try your options for
the code and get back to you later.

Regards,

Danyang


On 17-05-24 12:21 PM, Hong wrote:

Danyang :
I tested your data.
Your matrices encountered zero pivots, e.g.
petsc/src/ksp/ksp/examples/tutorials (master)
$ mpiexec -n 24 ./ex10 -f0 a_react_in_2.bin -rhs
b_react_in_2.bin -ksp_monitor -ksp_error_if_not_converged

[15]PETSC ERROR: Zero pivot in LU factorization:
http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
<http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot>
[15]PETSC ERROR: Zero pivot row 1249 value 2.05808e-14 tolerance
2.22045e-14
...

Adding option '-sub_pc_factor_shift_type nonzero', I got
mpiexec -n 24 ./ex10 -f0 a_react_in_2.bin -rhs b_react_in_2.bin
-ksp_monitor -ksp_error_if_not_converged
-sub_pc_factor_shift_type nonzero -mat_view ascii::ascii_info

Mat Object: 24 MPI processes
  type: mpiaij
  rows=45, cols=45
  total: nonzeros=6991400, allocated nonzeros=6991400
  total number of mallocs used during MatSetValues calls =0
not using I-node (on process 0) routines
  0 KSP Residual norm 5.84911755e+01
  1 KSP Residual norm 6.824179430230e-01
  2 KSP Residual norm 3.994483555787e-02
  3 KSP Residual norm 6.085841461433e-03
  4 KSP Residual norm 8.876162583511e-04
  5 KSP Residual norm 9.407780665278e-05
Number of iterations =   5
Residual norm 0.00542891

Hong

Hi Matt,

Yes. The matrix is 45x45 sparse. The hypre takes
hundreds of iterates, not for all but in most of the
timesteps. The matrix is not well conditioned, with nonzero
entries range from 1.0e-29 to 1.0e2. I also made double
check if there is anything wrong in the parallel version,
however, the matrix is the same with sequential version
except some round error which is relatively very small.
Usually for those not well conditioned matrix, direct solver
should be faster than iterative solver, right? But when I
use the sequential iterative solver with ILU prec developed
almost 20 years go by others, the solver converge fast with
appropriate factorization level. In other words, when I use
24 processor using hypre, the speed is almost the same as as
the old sequential iterative solver using 1 processor.

I use most of the default configuration for the general case
with pretty good speedup. And I am not sure if I miss
something for this problem.

Thanks,

Danyang


On 17-05-24 11:12 AM, Matthew Knepley wrote:

On Wed, May 24, 2017 at 12:50 PM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:

Hi Matthew and Barry,

Thanks for the quick response.

I also tried superlu and mumps, both work but it is
about four times slower than ILU(dt) prec through
hypre, with 24 processors I have tested.

You mean the total time is 4x? And you are taking hundreds
of iterates? That seems hard to believe, unless you are
dropping
a huge number of elements.

When I look into the convergence information, the
method using ILU(dt) still takes 200 to 3000 linear
iterations for each newton iteration. One reason is
this equation is hard to solve. As for the general
cases, the same method works awesome and get very good
speedup.

I do not understand what you mean here.

I also doubt if I use hypre correctly for this case. Is
there anyway to check this problem, or is it possible
to increase the factorization level through hypre?

I don't know.

  Matt

Thanks,

Danyang


On 17-05-24 04:59 AM, Matthew Knepley wrote:

        On Wed, May 24, 2017 at 2:21 AM, Danyang Su
<danyang...@gmail

Re: [petsc-users] Question on incomplete factorization level and fill


Dear Hong,

I just tested with different number of processors for the same matrix. 
It sometimes got "ERROR: Arguments are incompatible" for different 
number of processors. It works fine using 4, 8, or 24 processors, but 
failed with "ERROR: Arguments are incompatible" using 16 or 48 
processors. The error information is attached. I tested this on my local 
computer with 6 cores 12 threads. Any suggestion on this?


Thanks,

Danyang


On 17-05-24 12:28 PM, Danyang Su wrote:


Hi Hong,

Awesome. Thanks for testing the case. I will try your options for the 
code and get back to you later.


Regards,

Danyang


On 17-05-24 12:21 PM, Hong wrote:

Danyang :
I tested your data.
Your matrices encountered zero pivots, e.g.
petsc/src/ksp/ksp/examples/tutorials (master)
$ mpiexec -n 24 ./ex10 -f0 a_react_in_2.bin -rhs b_react_in_2.bin 
-ksp_monitor -ksp_error_if_not_converged


[15]PETSC ERROR: Zero pivot in LU factorization: 
http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
[15]PETSC ERROR: Zero pivot row 1249 value 2.05808e-14 tolerance 
2.22045e-14

...

Adding option '-sub_pc_factor_shift_type nonzero', I got
mpiexec -n 24 ./ex10 -f0 a_react_in_2.bin -rhs b_react_in_2.bin 
-ksp_monitor -ksp_error_if_not_converged -sub_pc_factor_shift_type 
nonzero -mat_view ascii::ascii_info


Mat Object: 24 MPI processes
  type: mpiaij
  rows=45, cols=45
  total: nonzeros=6991400, allocated nonzeros=6991400
  total number of mallocs used during MatSetValues calls =0
not using I-node (on process 0) routines
  0 KSP Residual norm 5.84911755e+01
  1 KSP Residual norm 6.824179430230e-01
  2 KSP Residual norm 3.994483555787e-02
  3 KSP Residual norm 6.085841461433e-03
  4 KSP Residual norm 8.876162583511e-04
  5 KSP Residual norm 9.407780665278e-05
Number of iterations =   5
Residual norm 0.00542891

Hong

Hi Matt,

Yes. The matrix is 45x45 sparse. The hypre takes hundreds
of iterates, not for all but in most of the timesteps. The matrix
is not well conditioned, with nonzero entries range from 1.0e-29
to 1.0e2. I also made double check if there is anything wrong in
the parallel version, however, the matrix is the same with
sequential version except some round error which is relatively
very small. Usually for those not well conditioned matrix, direct
solver should be faster than iterative solver, right? But when I
use the sequential iterative solver with ILU prec developed
almost 20 years go by others, the solver converge fast with
appropriate factorization level. In other words, when I use 24
processor using hypre, the speed is almost the same as as the old
sequential iterative solver using 1 processor.

I use most of the default configuration for the general case with
pretty good speedup. And I am not sure if I miss something for
this problem.

Thanks,

Danyang


On 17-05-24 11:12 AM, Matthew Knepley wrote:

On Wed, May 24, 2017 at 12:50 PM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:

Hi Matthew and Barry,

Thanks for the quick response.

I also tried superlu and mumps, both work but it is about
four times slower than ILU(dt) prec through hypre, with 24
processors I have tested.

You mean the total time is 4x? And you are taking hundreds of
iterates? That seems hard to believe, unless you are dropping
a huge number of elements.

When I look into the convergence information, the method
using ILU(dt) still takes 200 to 3000 linear iterations for
each newton iteration. One reason is this equation is hard
to solve. As for the general cases, the same method works
awesome and get very good speedup.

I do not understand what you mean here.

I also doubt if I use hypre correctly for this case. Is
there anyway to check this problem, or is it possible to
increase the factorization level through hypre?

I don't know.

  Matt

Thanks,

Danyang


On 17-05-24 04:59 AM, Matthew Knepley wrote:

On Wed, May 24, 2017 at 2:21 AM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:

Dear All,

I use PCFactorSetLevels for ILU and PCFactorSetFill for
other preconditioning in my code to help solve the
problems that the default option is hard to solve.
However, I found the latter one, PCFactorSetFill does
not take effect for my problem. The matrices and rhs as
well as the solutions are attached from the link below.
I obtain the solution using hypre preconditioner and it
takes 7 and 38 iterations for matrix 1 and matrix 2.
However, if I use other preconditioner, the solver just
failed at the first matrix. I have tested this matrix
u

Re: [petsc-users] Question on incomplete factorization level and fill


Hi Hong,

Awesome. Thanks for testing the case. I will try your options for the 
code and get back to you later.


Regards,

Danyang


On 17-05-24 12:21 PM, Hong wrote:

Danyang :
I tested your data.
Your matrices encountered zero pivots, e.g.
petsc/src/ksp/ksp/examples/tutorials (master)
$ mpiexec -n 24 ./ex10 -f0 a_react_in_2.bin -rhs b_react_in_2.bin 
-ksp_monitor -ksp_error_if_not_converged


[15]PETSC ERROR: Zero pivot in LU factorization: 
http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
[15]PETSC ERROR: Zero pivot row 1249 value 2.05808e-14 tolerance 
2.22045e-14

...

Adding option '-sub_pc_factor_shift_type nonzero', I got
mpiexec -n 24 ./ex10 -f0 a_react_in_2.bin -rhs b_react_in_2.bin 
-ksp_monitor -ksp_error_if_not_converged -sub_pc_factor_shift_type 
nonzero -mat_view ascii::ascii_info


Mat Object: 24 MPI processes
  type: mpiaij
  rows=45, cols=45
  total: nonzeros=6991400, allocated nonzeros=6991400
  total number of mallocs used during MatSetValues calls =0
not using I-node (on process 0) routines
  0 KSP Residual norm 5.84911755e+01
  1 KSP Residual norm 6.824179430230e-01
  2 KSP Residual norm 3.994483555787e-02
  3 KSP Residual norm 6.085841461433e-03
  4 KSP Residual norm 8.876162583511e-04
  5 KSP Residual norm 9.407780665278e-05
Number of iterations =   5
Residual norm 0.00542891

Hong

Hi Matt,

Yes. The matrix is 45x45 sparse. The hypre takes hundreds
of iterates, not for all but in most of the timesteps. The matrix
is not well conditioned, with nonzero entries range from 1.0e-29
to 1.0e2. I also made double check if there is anything wrong in
the parallel version, however, the matrix is the same with
sequential version except some round error which is relatively
very small. Usually for those not well conditioned matrix, direct
solver should be faster than iterative solver, right? But when I
use the sequential iterative solver with ILU prec developed almost
20 years go by others, the solver converge fast with appropriate
factorization level. In other words, when I use 24 processor using
hypre, the speed is almost the same as as the old sequential
iterative solver using 1 processor.

I use most of the default configuration for the general case with
pretty good speedup. And I am not sure if I miss something for
this problem.

Thanks,

Danyang


On 17-05-24 11:12 AM, Matthew Knepley wrote:

On Wed, May 24, 2017 at 12:50 PM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:

Hi Matthew and Barry,

Thanks for the quick response.

I also tried superlu and mumps, both work but it is about
four times slower than ILU(dt) prec through hypre, with 24
processors I have tested.

You mean the total time is 4x? And you are taking hundreds of
iterates? That seems hard to believe, unless you are dropping
a huge number of elements.

When I look into the convergence information, the method
using ILU(dt) still takes 200 to 3000 linear iterations for
each newton iteration. One reason is this equation is hard to
solve. As for the general cases, the same method works
awesome and get very good speedup.

I do not understand what you mean here.

I also doubt if I use hypre correctly for this case. Is there
anyway to check this problem, or is it possible to increase
the factorization level through hypre?

I don't know.

  Matt

Thanks,

Danyang


On 17-05-24 04:59 AM, Matthew Knepley wrote:

On Wed, May 24, 2017 at 2:21 AM, Danyang Su
<danyang...@gmail.com <mailto:danyang...@gmail.com>> wrote:

Dear All,

I use PCFactorSetLevels for ILU and PCFactorSetFill for
other preconditioning in my code to help solve the
problems that the default option is hard to solve.
However, I found the latter one, PCFactorSetFill does
not take effect for my problem. The matrices and rhs as
well as the solutions are attached from the link below.
I obtain the solution using hypre preconditioner and it
takes 7 and 38 iterations for matrix 1 and matrix 2.
However, if I use other preconditioner, the solver just
failed at the first matrix. I have tested this matrix
using the native sequential solver (not PETSc) with ILU
preconditioning. If I set the incomplete factorization
level to 0, this sequential solver will take more than
100 iterations. If I increase the factorization level to
1 or more, it just takes several iterations. This remind
me that the PC factor for this matrices should be
increased. However, when I tried it in PETSc, it just
does not work.

Matrix

Re: [petsc-users] Question on incomplete factorization level and fill


Hi Matt,

Yes. The matrix is 45x45 sparse. The hypre takes hundreds of 
iterates, not for all but in most of the timesteps. The matrix is not 
well conditioned, with nonzero entries range from 1.0e-29 to 1.0e2. I 
also made double check if there is anything wrong in the parallel 
version, however, the matrix is the same with sequential version except 
some round error which is relatively very small. Usually for those not 
well conditioned matrix, direct solver should be faster than iterative 
solver, right? But when I use the sequential iterative solver with ILU 
prec developed almost 20 years go by others, the solver converge fast 
with appropriate factorization level. In other words, when I use 24 
processor using hypre, the speed is almost the same as as the old 
sequential iterative solver using 1 processor.


I use most of the default configuration for the general case with pretty 
good speedup. And I am not sure if I miss something for this problem.


Thanks,

Danyang


On 17-05-24 11:12 AM, Matthew Knepley wrote:
On Wed, May 24, 2017 at 12:50 PM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Hi Matthew and Barry,

Thanks for the quick response.

I also tried superlu and mumps, both work but it is about four
times slower than ILU(dt) prec through hypre, with 24 processors I
have tested.

You mean the total time is 4x? And you are taking hundreds of 
iterates? That seems hard to believe, unless you are dropping

a huge number of elements.

When I look into the convergence information, the method using
ILU(dt) still takes 200 to 3000 linear iterations for each newton
iteration. One reason is this equation is hard to solve. As for
the general cases, the same method works awesome and get very good
speedup.

I do not understand what you mean here.

I also doubt if I use hypre correctly for this case. Is there
anyway to check this problem, or is it possible to increase the
factorization level through hypre?

I don't know.

  Matt

Thanks,

Danyang


On 17-05-24 04:59 AM, Matthew Knepley wrote:

On Wed, May 24, 2017 at 2:21 AM, Danyang Su <danyang...@gmail.com
<mailto:danyang...@gmail.com>> wrote:

Dear All,

I use PCFactorSetLevels for ILU and PCFactorSetFill for other
preconditioning in my code to help solve the problems that
the default option is hard to solve. However, I found the
latter one, PCFactorSetFill does not take effect for my
problem. The matrices and rhs as well as the solutions are
attached from the link below. I obtain the solution using
hypre preconditioner and it takes 7 and 38 iterations for
matrix 1 and matrix 2. However, if I use other
preconditioner, the solver just failed at the first matrix. I
have tested this matrix using the native sequential solver
(not PETSc) with ILU preconditioning. If I set the incomplete
factorization level to 0, this sequential solver will take
more than 100 iterations. If I increase the factorization
level to 1 or more, it just takes several iterations. This
remind me that the PC factor for this matrices should be
increased. However, when I tried it in PETSc, it just does
not work.

Matrix and rhs can be obtained from the link below.

https://eilinator.eos.ubc.ca:8443/index.php/s/CalUcq9CMeblk4R
<https://eilinator.eos.ubc.ca:8443/index.php/s/CalUcq9CMeblk4R>

Would anyone help to check if you can make this work by
increasing the PC factor level or fill?


We have ILU(k) supported in serial. However ILU(dt) which takes a
tolerance only works through Hypre

http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html
<http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html>

I recommend you try SuperLU or MUMPS, which can both be
downloaded automatically by configure, and
do a full sparse LU.

  Thanks,

Matt

Thanks and regards,

Danyang





-- 
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

http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/%7Emk51/>





--
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

http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/%7Emk51/>

Re: [petsc-users] Question on incomplete factorization level and fill


Hi Matthew and Barry,

Thanks for the quick response.

I also tried superlu and mumps, both work but it is about four times 
slower than ILU(dt) prec through hypre, with 24 processors I have 
tested. When I look into the convergence information, the method using 
ILU(dt) still takes 200 to 3000 linear iterations for each newton 
iteration. One reason is this equation is hard to solve. As for the 
general cases, the same method works awesome and get very good speedup. 
I also doubt if I use hypre correctly for this case. Is there anyway to 
check this problem, or is it possible to increase the factorization 
level through hypre?


Thanks,

Danyang


On 17-05-24 04:59 AM, Matthew Knepley wrote:
On Wed, May 24, 2017 at 2:21 AM, Danyang Su <danyang...@gmail.com 
<mailto:danyang...@gmail.com>> wrote:


Dear All,

I use PCFactorSetLevels for ILU and PCFactorSetFill for other
preconditioning in my code to help solve the problems that the
default option is hard to solve. However, I found the latter one,
PCFactorSetFill does not take effect for my problem. The matrices
and rhs as well as the solutions are attached from the link below.
I obtain the solution using hypre preconditioner and it takes 7
and 38 iterations for matrix 1 and matrix 2. However, if I use
other preconditioner, the solver just failed at the first matrix.
I have tested this matrix using the native sequential solver (not
PETSc) with ILU preconditioning. If I set the incomplete
factorization level to 0, this sequential solver will take more
than 100 iterations. If I increase the factorization level to 1 or
more, it just takes several iterations. This remind me that the PC
factor for this matrices should be increased. However, when I
tried it in PETSc, it just does not work.

Matrix and rhs can be obtained from the link below.

https://eilinator.eos.ubc.ca:8443/index.php/s/CalUcq9CMeblk4R
<https://eilinator.eos.ubc.ca:8443/index.php/s/CalUcq9CMeblk4R>

Would anyone help to check if you can make this work by increasing
the PC factor level or fill?


We have ILU(k) supported in serial. However ILU(dt) which takes a 
tolerance only works through Hypre


http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html

I recommend you try SuperLU or MUMPS, which can both be downloaded 
automatically by configure, and

do a full sparse LU.

  Thanks,

Matt

Thanks and regards,

Danyang





--
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

http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/%7Emk51/>

[petsc-users] Question on incomplete factorization level and fill