Re: [petsc-users] Petsc Section in DMPlex

[Nicholas Arnold-Medabalimi]

Hi

Thank you for your help with this barrage of questions. The AddField, while
nice for visualization, isn't the critical path for my development.

I've gotten the Filtered DM with its corresponding state vector sorted out.
Now I'm moving on to the mesh distribution. From your 2nd email in this
thread I've added
call DMPlexDistribute(dmplex_filtered, 1, distrib_sf, dmplex_distrib,ierr)
! adds section to dmplex_filtered and allocates vec_distrib using
DMCreateGlobalVector
call addSectionToDMPlex(dmplex_distrib,vec_distrib)
call DMGetGlobalSection(dmplex_distrib,distributedSection,ierr)
call DMPlexDistributeField(dmplex_disb, distrib_sf, filteredfieldSection,
vec_filtered, distributedSection, vec_distrib,ierr)

I'm not entirely clear if the DM I should feed into Distribute field should
be the starting or end DM but based on doc I think it should be the
destination. But I'm getting this error on run.
[1]PETSC ERROR: - Error Message
--
[1]PETSC ERROR: No support for this operation for this object type
[1]PETSC ERROR: PetscSectionSetUp is currently unsupported for
includesConstraints = PETSC_TRUE
[1]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.

I suspect this is occurring because I am incorrectly handling the 3
sets(full, filtered, distributed) of DMPLEX, Section, and Vec's.

Thanks
Nicholas

On Thu, Dec 8, 2022 at 12:50 PM Nicholas Arnold-Medabalimi <
narno...@umich.edu> wrote:

> I think I've figured out the issue. In previous efforts, I used
> DMAddField, which I think was key for the output to work properly.
>
>
> On Thu, Dec 8, 2022 at 10:48 AM Nicholas Arnold-Medabalimi <
> narno...@umich.edu> wrote:
>
>> Hi Matt
>>
>> Thanks. Found the issue just messed up the Fortran to C indexing.
>>
>> Another question. I have been using the Petsc VTK output to view things.
>> In some previous efforts, I used the PetscFVM object to set up my section
>> data. When I output vectors using that method in ParaView, I could view the
>> rank information and the state vector to visualize. As far as I can tell
>> when I do the same with Vectors that were created with my manually created
>> section the VecView using PETSCVIEWERVTK only mesh is output with the Rank
>> distribution is viewable (basically as if I had done a DM output instead of
>> a Vec output). My guess is this is because I'm not setting something in my
>> section Field setup properly for the VTK viewer to output it?
>>
>> For my section set up, I am calling
>> PetscSectionCreate
>> PetscSectionSetChart
>> PetscSectionSetFieldName
>> and then the appropriate PetscSectionSetDof
>> PetscSectionSetup
>> DMSetLocalSection
>>
>> Thanks again for all the help.
>>
>> Sincerely
>> Nicholas
>>
>>
>> On Thu, Dec 8, 2022 at 7:27 AM Matthew Knepley  wrote:
>>
>>> On Thu, Dec 8, 2022 at 3:04 AM Nicholas Arnold-Medabalimi <
>>> narno...@umich.edu> wrote:
>>>
 Hi Matt

 I think I've gotten it just about there. I'm just having an issue with
 the VecISCopy. I have an IS built that matches size correctly to map from
 the full state to the filtered state. The core issue I think, is should the
 expanded IS the ownership range of the vector subtracted out. Looking at
 the implementation, it looks like VecISCopy takes care of that for me.
 (Line 573 in src/vec/vec/utils/projection.c) But I could be mistaken.

>>>
>>> It is a good question. We have tried to give guidance on the manpage:
>>>
>>>   The index set identifies entries in the global vector. Negative
>>> indices are skipped; indices outside the ownership range of vfull will
>>> raise an error.
>>>
>>> which means that it expects _global_ indices, and you have retrieved the
>>> global section, so that matches.
>>> The calculation of the index size looks right to me, and so does the
>>> index calculation.
>>>
>>> I would put a check in the loop, making sure that the calculated indices
>>> lie within [oStart, oEnd). The global
>>> section is designed to ensure that. It is not clear why one would lie
>>> outside.
>>>
>>> When I am debugging, I run a very small problem, and print out all the
>>> sections.
>>>
>>>   Thanks,
>>>
>>>  Matt
>>>
>>>
 The error I am getting is:
 [0]PETSC ERROR: - Error Message
 --
 [0]PETSC ERROR: No support for this operation for this object type
 [0]PETSC ERROR: Only owned values supported


 Here is what I am currently doing.

 call DMPlexFilter(dmplex_full, iBlankLabel, 1, dmplex_filtered,ierr)
 call DMPlexGetSubpointIS(dmplex_filtered, subpointsIS,ierr)

 ! adds section to dmplex_filtered and allocates vec_filtered using
 DMCreateGlobalVector
 call addSectionToDMPlex(dmplex_filtered,vec_filtered)

 ! Get Sections for dmplex_filtered and dmplex_full
 call 

[petsc-users] Modifying Entries Tied to Specific Points in Finite Element Stiffness Matrix using DMPlex

[Hongrui Yu]

Hello! I'm trying to adapt a serial Finite Element code using PETSc. In this
code it reads in special stiffness terms between the boundary DoFs from an
input file, and add them to corresponding locations in the global Jacobian
matrix. I currently use a DM Plex object to store the mesh information. My
understanding is that once the DM is distributed its points are renumbered
across different ranks. I wonder if there is a good way to find the
corresponding entries that needs to be modified in the global Jacobian
matrix? 

 

For Vectors I'm currently creating a Natural Vector and simply do
DMPlexNaturalToGlobal. Is there a way to create a "Natural Mat" just like
"Natural Vector" and then do some sort of NaturalToGlobal for this Mat? 

 

Any help would be highly appreciated!

 

Kevin

Re: [petsc-users] Fortran Interface NULL object / Casting

[Barry Smith]

   You would use PETSC_NULL_DMLABEL but Matt needs to customize the PETSc 
Fortran stub for DMAddField() for you to handle accepting the NULL from PETSc.

   Barry





> On Dec 8, 2022, at 1:05 PM, Nicholas Arnold-Medabalimi  
> wrote:
> 
> Hi Petsc Users
> 
> I am trying to use DMAddField in a Fortran code. I had some questions on 
> casting/passing NULL. I follow how to pass NULL for standard types (INT, 
> CHAR, etc). 
> 
> Is there a method/best practice for passing NULL for Petsc type arguments? 
> (In this case DMAddLabel I'd want to pass NULL to a DMLabel Parameter not an 
> intrinsic type)
> 
> For the 2nd question, In some cases in C we would cast a more specific object 
> to Petsc Object 
> DMAddField(dm, NULL, (PetscObject)fvm); (where fvm is a PetscFVM type)
> Is there a method or practice to do the same in Fortran.
> 
> Thanks
> Nicholas
> 
> -- 
> Nicholas Arnold-Medabalimi
> 
> Ph.D. Candidate
> Computational Aeroscience Lab
> University of Michigan

Re: [petsc-users] prevent linking to multithreaded BLAS?

[Jed Brown]

Barry Smith  writes:

>> We could test at runtime whether child threads exist/are created when 
>> calling BLAS and deliver a warning. 
>
>   How does one test for this? Some standard Unix API for checking this?

I'm not sure, the ids of child threads are in /proc/$pid/task/ and (when opened 
by a process) /proc/self/task/. See man procfs for details.

[petsc-users] Fortran Interface NULL object / Casting

[Nicholas Arnold-Medabalimi]

Hi Petsc Users

I am trying to use DMAddField in a Fortran code. I had some questions on
casting/passing NULL. I follow how to pass NULL for standard types (INT,
CHAR, etc).

Is there a method/best practice for passing NULL for Petsc type arguments?
(In this case DMAddLabel I'd want to pass NULL to a DMLabel Parameter not
an intrinsic type)

For the 2nd question, In some cases in C we would cast a more specific
object to Petsc Object
DMAddField(dm, NULL, (PetscObject)fvm); (where fvm is a PetscFVM type)
Is there a method or practice to do the same in Fortran.

Thanks
Nicholas

-- 
Nicholas Arnold-Medabalimi

Ph.D. Candidate
Computational Aeroscience Lab
University of Michigan

Re: [petsc-users] Petsc Section in DMPlex

[Nicholas Arnold-Medabalimi]

I think I've figured out the issue. In previous efforts, I used DMAddField,
which I think was key for the output to work properly.


On Thu, Dec 8, 2022 at 10:48 AM Nicholas Arnold-Medabalimi <
narno...@umich.edu> wrote:

> Hi Matt
>
> Thanks. Found the issue just messed up the Fortran to C indexing.
>
> Another question. I have been using the Petsc VTK output to view things.
> In some previous efforts, I used the PetscFVM object to set up my section
> data. When I output vectors using that method in ParaView, I could view the
> rank information and the state vector to visualize. As far as I can tell
> when I do the same with Vectors that were created with my manually created
> section the VecView using PETSCVIEWERVTK only mesh is output with the Rank
> distribution is viewable (basically as if I had done a DM output instead of
> a Vec output). My guess is this is because I'm not setting something in my
> section Field setup properly for the VTK viewer to output it?
>
> For my section set up, I am calling
> PetscSectionCreate
> PetscSectionSetChart
> PetscSectionSetFieldName
> and then the appropriate PetscSectionSetDof
> PetscSectionSetup
> DMSetLocalSection
>
> Thanks again for all the help.
>
> Sincerely
> Nicholas
>
>
> On Thu, Dec 8, 2022 at 7:27 AM Matthew Knepley  wrote:
>
>> On Thu, Dec 8, 2022 at 3:04 AM Nicholas Arnold-Medabalimi <
>> narno...@umich.edu> wrote:
>>
>>> Hi Matt
>>>
>>> I think I've gotten it just about there. I'm just having an issue with
>>> the VecISCopy. I have an IS built that matches size correctly to map from
>>> the full state to the filtered state. The core issue I think, is should the
>>> expanded IS the ownership range of the vector subtracted out. Looking at
>>> the implementation, it looks like VecISCopy takes care of that for me.
>>> (Line 573 in src/vec/vec/utils/projection.c) But I could be mistaken.
>>>
>>
>> It is a good question. We have tried to give guidance on the manpage:
>>
>>   The index set identifies entries in the global vector. Negative indices
>> are skipped; indices outside the ownership range of vfull will raise an
>> error.
>>
>> which means that it expects _global_ indices, and you have retrieved the
>> global section, so that matches.
>> The calculation of the index size looks right to me, and so does the
>> index calculation.
>>
>> I would put a check in the loop, making sure that the calculated indices
>> lie within [oStart, oEnd). The global
>> section is designed to ensure that. It is not clear why one would lie
>> outside.
>>
>> When I am debugging, I run a very small problem, and print out all the
>> sections.
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> The error I am getting is:
>>> [0]PETSC ERROR: - Error Message
>>> --
>>> [0]PETSC ERROR: No support for this operation for this object type
>>> [0]PETSC ERROR: Only owned values supported
>>>
>>>
>>> Here is what I am currently doing.
>>>
>>> call DMPlexFilter(dmplex_full, iBlankLabel, 1, dmplex_filtered,ierr)
>>> call DMPlexGetSubpointIS(dmplex_filtered, subpointsIS,ierr)
>>>
>>> ! adds section to dmplex_filtered and allocates vec_filtered using
>>> DMCreateGlobalVector
>>> call addSectionToDMPlex(dmplex_filtered,vec_filtered)
>>>
>>> ! Get Sections for dmplex_filtered and dmplex_full
>>> call DMGetGlobalSection(dmplex_filtered,filteredfieldSection,ierr)
>>> call DMGetGlobalSection(dmplex_full,fullfieldSection,ierr)
>>>
>>>
>>>
>>> call ISGetIndicesF90(subpointsIS, subPointKey,ierr)
>>> ExpandedIndexSize = 0
>>> do i = 1, size(subPointKey)
>>> call PetscSectionGetDof(fullfieldSection, subPointKey(i), dof,ierr)
>>> ExpandedIndexSize = ExpandedIndexSize + dof
>>> enddo
>>>
>>>
>>> !Create expandedIS from offset sections of full and filtered sections
>>> allocate(ExpandedIndex(ExpandedIndexSize))
>>> call VecGetOwnershipRange(vec_full,oStart,oEnd,ierr)
>>> do i = 1, size(subPointKey)
>>> call PetscSectionGetOffset(fullfieldSection, subPointKey(i),
>>> offset,ierr)
>>> call PetscSectionGetDof(fullfieldSection, subPointKey(i), dof,ierr)
>>> !offset=offset-oStart   !looking at VecIScopy it takes care of this
>>> subtraction (not sure)
>>> do j = 1, (dof)
>>> ExpandedIndex((i-1)*dof+j) = offset+j
>>> end do
>>> enddo
>>>
>>> call ISCreateGeneral(PETSC_COMM_WORLD, ExpandedIndexSize,
>>> ExpandedIndex, PETSC_COPY_VALUES, expandedIS,ierr)
>>> call ISRestoreIndicesF90(subpointsIS, subPointKey,ierr)
>>> deallocate(ExpandedIndex)
>>>
>>>
>>> call VecGetLocalSize(vec_full,sizeVec,ierr)
>>> write(*,*) sizeVec
>>> call VecGetLocalSize(vec_filtered,sizeVec,ierr)
>>> write(*,*) sizeVec
>>> call ISGetLocalSize(expandedIS,sizeVec,ierr)
>>> write(*,*) sizeVec
>>> call PetscSynchronizedFlush(PETSC_COMM_WORLD,ierr)
>>>
>>>
>>> call VecISCopy(vec_full,expandedIS,SCATTER_REVERSE,vec_filtered,ierr)
>>>
>>>
>>> Thanks again for the great help.
>>>
>>> Sincerely
>>> Nicholas
>>>
>>>
>>> On Wed, Dec 7, 2022 

Re: [petsc-users] prevent linking to multithreaded BLAS?

[Barry Smith]

> On Dec 7, 2022, at 11:56 PM, Jed Brown  wrote:
> 
> It isn't always wrong to link threaded BLAS. For example, a user might need 
> to call threaded BLAS on the side (but the application can only link one) or 
> a sparse direct solver might want threading for the supernode.

  Indeed, the user asked specifically for their work flow if configure could, 
based on additional configure argument, ensure that they did not get a threaded 
BLAS; they did not ask that configure never give a threaded BLAS or even that 
it give a non-threaded BLAS by default.


> We could test at runtime whether child threads exist/are created when calling 
> BLAS and deliver a warning. 

  How does one test for this? Some standard Unix API for checking this?

> 
> Barry Smith  writes:
> 
>>   There would need to be, for example, some symbol in all the threaded BLAS 
>> libraries that is not in the unthreaded libraries. Of at least in some of 
>> the threaded libraries but never in the unthreaded. 
>> 
>>   BlasLapack.py could check for the special symbol(s) to determine.
>> 
>>  Barry
>> 
>> 
>>> On Dec 7, 2022, at 4:47 PM, Mark Lohry  wrote:
>>> 
>>> Thanks, yes, I figured out the OMP_NUM_THREADS=1 way while triaging it, and 
>>> the --download-fblaslapack way occurred to me. 
>>> 
>>> I was hoping for something that "just worked" (refuse to build in this 
>>> case) but I don't know if it's programmatically possible for petsc to tell 
>>> whether or not it's linking to a threaded BLAS?
>>> 
>>> On Wed, Dec 7, 2022 at 4:35 PM Satish Balay >> > wrote:
 If you don't specify a blas to use - petsc configure will guess and use 
 what it can find.
 
 So only way to force it use a particular blas is to specify one [one way 
 is --download-fblaslapack]
 
 Wrt multi-thread openblas -  you can force it run single threaded [by one 
 of these 2 env variables]
 
# Use single thread openblas
export OPENBLAS_NUM_THREADS=1
export OMP_NUM_THREADS=1
 
 Satish
 
 
 On Wed, 7 Dec 2022, Mark Lohry wrote:
 
> I ran into an unexpected issue -- on an NP-core machine, each MPI rank of
> my application was launching NP threads, such that when running with
> multiple ranks the machine was quickly oversubscribed and performance
> tanked.
> 
> The root cause of this was petsc linking against the system-provided
> library (libopenblas0-pthread in this case) set by the update-alternatives
> in ubuntu. At some point this machine got updated to using the threaded
> blas implementation instead of serial; not sure how, and I wouldn't have
> noticed if I weren't running interactively.
> 
> Is there any mechanism in petsc or its build system to prevent linking
> against an inappropriate BLAS, or do I need to be diligent about manually
> setting the BLAS library in the configuration stage?
> 
> Thanks,
> Mark
> 
 

Re: [petsc-users] Petsc Section in DMPlex

[Nicholas Arnold-Medabalimi]

Hi Matt

Thanks. Found the issue just messed up the Fortran to C indexing.

Another question. I have been using the Petsc VTK output to view things. In
some previous efforts, I used the PetscFVM object to set up my section
data. When I output vectors using that method in ParaView, I could view the
rank information and the state vector to visualize. As far as I can tell
when I do the same with Vectors that were created with my manually created
section the VecView using PETSCVIEWERVTK only mesh is output with the Rank
distribution is viewable (basically as if I had done a DM output instead of
a Vec output). My guess is this is because I'm not setting something in my
section Field setup properly for the VTK viewer to output it?

For my section set up, I am calling
PetscSectionCreate
PetscSectionSetChart
PetscSectionSetFieldName
and then the appropriate PetscSectionSetDof
PetscSectionSetup
DMSetLocalSection

Thanks again for all the help.

Sincerely
Nicholas


On Thu, Dec 8, 2022 at 7:27 AM Matthew Knepley  wrote:

> On Thu, Dec 8, 2022 at 3:04 AM Nicholas Arnold-Medabalimi <
> narno...@umich.edu> wrote:
>
>> Hi Matt
>>
>> I think I've gotten it just about there. I'm just having an issue with
>> the VecISCopy. I have an IS built that matches size correctly to map from
>> the full state to the filtered state. The core issue I think, is should the
>> expanded IS the ownership range of the vector subtracted out. Looking at
>> the implementation, it looks like VecISCopy takes care of that for me.
>> (Line 573 in src/vec/vec/utils/projection.c) But I could be mistaken.
>>
>
> It is a good question. We have tried to give guidance on the manpage:
>
>   The index set identifies entries in the global vector. Negative indices
> are skipped; indices outside the ownership range of vfull will raise an
> error.
>
> which means that it expects _global_ indices, and you have retrieved the
> global section, so that matches.
> The calculation of the index size looks right to me, and so does the index
> calculation.
>
> I would put a check in the loop, making sure that the calculated indices
> lie within [oStart, oEnd). The global
> section is designed to ensure that. It is not clear why one would lie
> outside.
>
> When I am debugging, I run a very small problem, and print out all the
> sections.
>
>   Thanks,
>
>  Matt
>
>
>> The error I am getting is:
>> [0]PETSC ERROR: - Error Message
>> --
>> [0]PETSC ERROR: No support for this operation for this object type
>> [0]PETSC ERROR: Only owned values supported
>>
>>
>> Here is what I am currently doing.
>>
>> call DMPlexFilter(dmplex_full, iBlankLabel, 1, dmplex_filtered,ierr)
>> call DMPlexGetSubpointIS(dmplex_filtered, subpointsIS,ierr)
>>
>> ! adds section to dmplex_filtered and allocates vec_filtered using
>> DMCreateGlobalVector
>> call addSectionToDMPlex(dmplex_filtered,vec_filtered)
>>
>> ! Get Sections for dmplex_filtered and dmplex_full
>> call DMGetGlobalSection(dmplex_filtered,filteredfieldSection,ierr)
>> call DMGetGlobalSection(dmplex_full,fullfieldSection,ierr)
>>
>>
>>
>> call ISGetIndicesF90(subpointsIS, subPointKey,ierr)
>> ExpandedIndexSize = 0
>> do i = 1, size(subPointKey)
>> call PetscSectionGetDof(fullfieldSection, subPointKey(i), dof,ierr)
>> ExpandedIndexSize = ExpandedIndexSize + dof
>> enddo
>>
>>
>> !Create expandedIS from offset sections of full and filtered sections
>> allocate(ExpandedIndex(ExpandedIndexSize))
>> call VecGetOwnershipRange(vec_full,oStart,oEnd,ierr)
>> do i = 1, size(subPointKey)
>> call PetscSectionGetOffset(fullfieldSection, subPointKey(i),
>> offset,ierr)
>> call PetscSectionGetDof(fullfieldSection, subPointKey(i), dof,ierr)
>> !offset=offset-oStart   !looking at VecIScopy it takes care of this
>> subtraction (not sure)
>> do j = 1, (dof)
>> ExpandedIndex((i-1)*dof+j) = offset+j
>> end do
>> enddo
>>
>> call ISCreateGeneral(PETSC_COMM_WORLD, ExpandedIndexSize, ExpandedIndex,
>> PETSC_COPY_VALUES, expandedIS,ierr)
>> call ISRestoreIndicesF90(subpointsIS, subPointKey,ierr)
>> deallocate(ExpandedIndex)
>>
>>
>> call VecGetLocalSize(vec_full,sizeVec,ierr)
>> write(*,*) sizeVec
>> call VecGetLocalSize(vec_filtered,sizeVec,ierr)
>> write(*,*) sizeVec
>> call ISGetLocalSize(expandedIS,sizeVec,ierr)
>> write(*,*) sizeVec
>> call PetscSynchronizedFlush(PETSC_COMM_WORLD,ierr)
>>
>>
>> call VecISCopy(vec_full,expandedIS,SCATTER_REVERSE,vec_filtered,ierr)
>>
>>
>> Thanks again for the great help.
>>
>> Sincerely
>> Nicholas
>>
>>
>> On Wed, Dec 7, 2022 at 9:29 PM Matthew Knepley  wrote:
>>
>>> On Wed, Dec 7, 2022 at 9:21 PM Nicholas Arnold-Medabalimi <
>>> narno...@umich.edu> wrote:
>>>
 Thank you for the help.

 I think the last piece of the puzzle is how do I create the "expanded
 IS" from the subpoint IS using the section?

>>>
>>> Loop over the points in the IS. For each point, get the dof and offset

Re: [petsc-users] [EXTERNAL] Re: Kokkos backend for Mat and Vec diverging when running on CUDA device.

[Fackler, Philip via petsc-users]

Great! Thank you!

Philip Fackler
Research Software Engineer, Application Engineering Group
Advanced Computing Systems Research Section
Computer Science and Mathematics Division
Oak Ridge National Laboratory

From: Junchao Zhang 
Sent: Wednesday, December 7, 2022 18:47
To: Fackler, Philip 
Cc: xolotl-psi-developm...@lists.sourceforge.net 
; petsc-users@mcs.anl.gov 
; Blondel, Sophie ; Roth, Philip 

Subject: Re: [EXTERNAL] Re: [petsc-users] Kokkos backend for Mat and Vec 
diverging when running on CUDA device.

Hi, Philip,
 I could reproduce the error. I need to find a  way to debug it.  Thanks.

/home/jczhang/xolotl/test/system/SystemTestCase.cpp(317): fatal error: in 
"System/PSI_1": absolute value of diffNorm{0.19704848134353209} exceeds 1e-10
*** 1 failure is detected in the test module "Regression"

--Junchao Zhang


On Tue, Dec 6, 2022 at 10:10 AM Fackler, Philip 
mailto:fackle...@ornl.gov>> wrote:
I think it would be simpler to use the develop branch for this issue. But you 
can still just build the SystemTester. Then (if you changed the PSI_1 case) run:

 ./test/system/SystemTester -t System/PSI_1 -- -v​

(No need for multiple MPI ranks)

Thanks,

Philip Fackler
Research Software Engineer, Application Engineering Group
Advanced Computing Systems Research Section
Computer Science and Mathematics Division
Oak Ridge National Laboratory

From: Junchao Zhang mailto:junchao.zh...@gmail.com>>
Sent: Monday, December 5, 2022 15:40
To: Fackler, Philip mailto:fackle...@ornl.gov>>
Cc: 
xolotl-psi-developm...@lists.sourceforge.net
 
mailto:xolotl-psi-developm...@lists.sourceforge.net>>;
 petsc-users@mcs.anl.gov 
mailto:petsc-users@mcs.anl.gov>>; Blondel, Sophie 
mailto:sblon...@utk.edu>>; Roth, Philip 
mailto:rot...@ornl.gov>>
Subject: Re: [EXTERNAL] Re: [petsc-users] Kokkos backend for Mat and Vec 
diverging when running on CUDA device.

I configured with xolotl branch feature-petsc-kokkos, and typed `make` under 
~/xolotl-build/.  Though there were errors,  a lot of *Tester were built.
[ 62%] Built target xolotlViz
[ 63%] Linking CXX executable TemperatureProfileHandlerTester
[ 64%] Linking CXX executable TemperatureGradientHandlerTester
[ 64%] Built target TemperatureProfileHandlerTester
[ 64%] Built target TemperatureConstantHandlerTester
[ 64%] Built target TemperatureGradientHandlerTester
[ 65%] Linking CXX executable HeatEquationHandlerTester
[ 65%] Built target HeatEquationHandlerTester
[ 66%] Linking CXX executable FeFitFluxHandlerTester
[ 66%] Linking CXX executable W111FitFluxHandlerTester
[ 67%] Linking CXX executable FuelFitFluxHandlerTester
[ 67%] Linking CXX executable W211FitFluxHandlerTester
Which Tester should I use to run with the parameter file 
benchmarks/params_system_PSI_2.txt? And how many ranks should I use?  Could you 
give an example command line?
Thanks.

--Junchao Zhang


On Mon, Dec 5, 2022 at 2:22 PM Junchao Zhang 
mailto:junchao.zh...@gmail.com>> wrote:
Hello, Philip,
   Do I still need to use the feature-petsc-kokkos branch?
--Junchao Zhang


On Mon, Dec 5, 2022 at 11:08 AM Fackler, Philip 
mailto:fackle...@ornl.gov>> wrote:
Junchao,

Thank you for working on this. If you open the parameter file for, say, the 
PSI_2 system test case (benchmarks/params_system_PSI_2.txt), simply add 
-dm_mat_type aijkokkos -dm_vec_type kokkos​` to the "petscArgs=" field (or the 
corresponding cusparse/cuda option).

Thanks,

Philip Fackler
Research Software Engineer, Application Engineering Group
Advanced Computing Systems Research Section
Computer Science and Mathematics Division
Oak Ridge National Laboratory

From: Junchao Zhang mailto:junchao.zh...@gmail.com>>
Sent: Thursday, December 1, 2022 17:05
To: Fackler, Philip mailto:fackle...@ornl.gov>>
Cc: 
xolotl-psi-developm...@lists.sourceforge.net
 
mailto:xolotl-psi-developm...@lists.sourceforge.net>>;
 petsc-users@mcs.anl.gov 
mailto:petsc-users@mcs.anl.gov>>; Blondel, Sophie 
mailto:sblon...@utk.edu>>; Roth, Philip 
mailto:rot...@ornl.gov>>
Subject: Re: [EXTERNAL] Re: [petsc-users] Kokkos backend for Mat and Vec 
diverging when running on CUDA device.

Hi, Philip,
  Sorry for the long delay.  I could not get something useful from the 
-log_view output.  Since I have already built xolotl, could you give me 
instructions on how to do a xolotl test to reproduce the divergence with petsc 
GPU backends (but fine on CPU)?
  Thank you.
--Junchao Zhang


On Wed, Nov 16, 2022 at 1:38 PM Fackler, Philip 
mailto:fackle...@ornl.gov>> wrote:
-- PETSc 
Performance Summary: 
--

Unknown Name on a  named PC0115427 with 1 processor, by 4pf Wed Nov 16 14:36:46 
2022
Using 

Re: [petsc-users] Petsc Section in DMPlex

[Matthew Knepley]

On Thu, Dec 8, 2022 at 3:04 AM Nicholas Arnold-Medabalimi <
narno...@umich.edu> wrote:

> Hi Matt
>
> I think I've gotten it just about there. I'm just having an issue with the
> VecISCopy. I have an IS built that matches size correctly to map from the
> full state to the filtered state. The core issue I think, is should the
> expanded IS the ownership range of the vector subtracted out. Looking at
> the implementation, it looks like VecISCopy takes care of that for me.
> (Line 573 in src/vec/vec/utils/projection.c) But I could be mistaken.
>

It is a good question. We have tried to give guidance on the manpage:

  The index set identifies entries in the global vector. Negative indices
are skipped; indices outside the ownership range of vfull will raise an
error.

which means that it expects _global_ indices, and you have retrieved the
global section, so that matches.
The calculation of the index size looks right to me, and so does the index
calculation.

I would put a check in the loop, making sure that the calculated indices
lie within [oStart, oEnd). The global
section is designed to ensure that. It is not clear why one would lie
outside.

When I am debugging, I run a very small problem, and print out all the
sections.

  Thanks,

 Matt


> The error I am getting is:
> [0]PETSC ERROR: - Error Message
> --
> [0]PETSC ERROR: No support for this operation for this object type
> [0]PETSC ERROR: Only owned values supported
>
>
> Here is what I am currently doing.
>
> call DMPlexFilter(dmplex_full, iBlankLabel, 1, dmplex_filtered,ierr)
> call DMPlexGetSubpointIS(dmplex_filtered, subpointsIS,ierr)
>
> ! adds section to dmplex_filtered and allocates vec_filtered using
> DMCreateGlobalVector
> call addSectionToDMPlex(dmplex_filtered,vec_filtered)
>
> ! Get Sections for dmplex_filtered and dmplex_full
> call DMGetGlobalSection(dmplex_filtered,filteredfieldSection,ierr)
> call DMGetGlobalSection(dmplex_full,fullfieldSection,ierr)
>
>
>
> call ISGetIndicesF90(subpointsIS, subPointKey,ierr)
> ExpandedIndexSize = 0
> do i = 1, size(subPointKey)
> call PetscSectionGetDof(fullfieldSection, subPointKey(i), dof,ierr)
> ExpandedIndexSize = ExpandedIndexSize + dof
> enddo
>
>
> !Create expandedIS from offset sections of full and filtered sections
> allocate(ExpandedIndex(ExpandedIndexSize))
> call VecGetOwnershipRange(vec_full,oStart,oEnd,ierr)
> do i = 1, size(subPointKey)
> call PetscSectionGetOffset(fullfieldSection, subPointKey(i),
> offset,ierr)
> call PetscSectionGetDof(fullfieldSection, subPointKey(i), dof,ierr)
> !offset=offset-oStart   !looking at VecIScopy it takes care of this
> subtraction (not sure)
> do j = 1, (dof)
> ExpandedIndex((i-1)*dof+j) = offset+j
> end do
> enddo
>
> call ISCreateGeneral(PETSC_COMM_WORLD, ExpandedIndexSize, ExpandedIndex,
> PETSC_COPY_VALUES, expandedIS,ierr)
> call ISRestoreIndicesF90(subpointsIS, subPointKey,ierr)
> deallocate(ExpandedIndex)
>
>
> call VecGetLocalSize(vec_full,sizeVec,ierr)
> write(*,*) sizeVec
> call VecGetLocalSize(vec_filtered,sizeVec,ierr)
> write(*,*) sizeVec
> call ISGetLocalSize(expandedIS,sizeVec,ierr)
> write(*,*) sizeVec
> call PetscSynchronizedFlush(PETSC_COMM_WORLD,ierr)
>
>
> call VecISCopy(vec_full,expandedIS,SCATTER_REVERSE,vec_filtered,ierr)
>
>
> Thanks again for the great help.
>
> Sincerely
> Nicholas
>
>
> On Wed, Dec 7, 2022 at 9:29 PM Matthew Knepley  wrote:
>
>> On Wed, Dec 7, 2022 at 9:21 PM Nicholas Arnold-Medabalimi <
>> narno...@umich.edu> wrote:
>>
>>> Thank you for the help.
>>>
>>> I think the last piece of the puzzle is how do I create the "expanded
>>> IS" from the subpoint IS using the section?
>>>
>>
>> Loop over the points in the IS. For each point, get the dof and offset
>> from the Section. Make a new IS that has all the
>> dogs, namely each run [offset, offset+dof).
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> Sincerely
>>> Nicholas
>>>
>>>
>>> On Wed, Dec 7, 2022 at 7:06 AM Matthew Knepley 
>>> wrote:
>>>
 On Wed, Dec 7, 2022 at 6:51 AM Nicholas Arnold-Medabalimi <
 narno...@umich.edu> wrote:

> Hi
>
> Thank you so much for your patience. One thing to note: I don't have
> any need to go back from the filtered distributed mapping back to the full
> but it is good to know.
>
> One aside question.
> 1) Is natural and global ordering the same in this context?
>

 No.


> As far as implementing what you have described.
>
> When I call ISView on the generated SubpointIS, I get an unusual error
> which I'm not sure how to interpret. (this case is running on 2 ranks and
> the filter label has points located on both ranks of the original DM.
> However, if I manually get the indices (the commented lines), it seems to
> not have any issues.
> call DMPlexFilter(dmplex_full, iBlankLabel, 1, dmplex_filtered,ierr)
> call 

Re: [petsc-users] Petsc Section in DMPlex

[Nicholas Arnold-Medabalimi]

Hi Matt

I think I've gotten it just about there. I'm just having an issue with the
VecISCopy. I have an IS built that matches size correctly to map from the
full state to the filtered state. The core issue I think, is should the
expanded IS the ownership range of the vector subtracted out. Looking at
the implementation, it looks like VecISCopy takes care of that for me.
(Line 573 in src/vec/vec/utils/projection.c) But I could be mistaken.

The error I am getting is:
[0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: Only owned values supported


Here is what I am currently doing.

call DMPlexFilter(dmplex_full, iBlankLabel, 1, dmplex_filtered,ierr)
call DMPlexGetSubpointIS(dmplex_filtered, subpointsIS,ierr)

! adds section to dmplex_filtered and allocates vec_filtered using
DMCreateGlobalVector
call addSectionToDMPlex(dmplex_filtered,vec_filtered)

! Get Sections for dmplex_filtered and dmplex_full
call DMGetGlobalSection(dmplex_filtered,filteredfieldSection,ierr)
call DMGetGlobalSection(dmplex_full,fullfieldSection,ierr)



call ISGetIndicesF90(subpointsIS, subPointKey,ierr)
ExpandedIndexSize = 0
do i = 1, size(subPointKey)
call PetscSectionGetDof(fullfieldSection, subPointKey(i), dof,ierr)
ExpandedIndexSize = ExpandedIndexSize + dof
enddo


!Create expandedIS from offset sections of full and filtered sections
allocate(ExpandedIndex(ExpandedIndexSize))
call VecGetOwnershipRange(vec_full,oStart,oEnd,ierr)
do i = 1, size(subPointKey)
call PetscSectionGetOffset(fullfieldSection, subPointKey(i),
offset,ierr)
call PetscSectionGetDof(fullfieldSection, subPointKey(i), dof,ierr)
!offset=offset-oStart   !looking at VecIScopy it takes care of this
subtraction (not sure)
do j = 1, (dof)
ExpandedIndex((i-1)*dof+j) = offset+j
end do
enddo

call ISCreateGeneral(PETSC_COMM_WORLD, ExpandedIndexSize, ExpandedIndex,
PETSC_COPY_VALUES, expandedIS,ierr)
call ISRestoreIndicesF90(subpointsIS, subPointKey,ierr)
deallocate(ExpandedIndex)


call VecGetLocalSize(vec_full,sizeVec,ierr)
write(*,*) sizeVec
call VecGetLocalSize(vec_filtered,sizeVec,ierr)
write(*,*) sizeVec
call ISGetLocalSize(expandedIS,sizeVec,ierr)
write(*,*) sizeVec
call PetscSynchronizedFlush(PETSC_COMM_WORLD,ierr)


call VecISCopy(vec_full,expandedIS,SCATTER_REVERSE,vec_filtered,ierr)


Thanks again for the great help.

Sincerely
Nicholas


On Wed, Dec 7, 2022 at 9:29 PM Matthew Knepley  wrote:

> On Wed, Dec 7, 2022 at 9:21 PM Nicholas Arnold-Medabalimi <
> narno...@umich.edu> wrote:
>
>> Thank you for the help.
>>
>> I think the last piece of the puzzle is how do I create the "expanded IS"
>> from the subpoint IS using the section?
>>
>
> Loop over the points in the IS. For each point, get the dof and offset
> from the Section. Make a new IS that has all the
> dogs, namely each run [offset, offset+dof).
>
>   Thanks,
>
>  Matt
>
>
>> Sincerely
>> Nicholas
>>
>>
>> On Wed, Dec 7, 2022 at 7:06 AM Matthew Knepley  wrote:
>>
>>> On Wed, Dec 7, 2022 at 6:51 AM Nicholas Arnold-Medabalimi <
>>> narno...@umich.edu> wrote:
>>>
 Hi

 Thank you so much for your patience. One thing to note: I don't have
 any need to go back from the filtered distributed mapping back to the full
 but it is good to know.

 One aside question.
 1) Is natural and global ordering the same in this context?

>>>
>>> No.
>>>
>>>
 As far as implementing what you have described.

 When I call ISView on the generated SubpointIS, I get an unusual error
 which I'm not sure how to interpret. (this case is running on 2 ranks and
 the filter label has points located on both ranks of the original DM.
 However, if I manually get the indices (the commented lines), it seems to
 not have any issues.
 call DMPlexFilter(dmplex_full, iBlankLabel, 1, dmplex_filtered,ierr)
 call DMPlexGetSubpointIS(dmplex_filtered, subpointsIS,ierr)
 !call ISGetIndicesF90(subpointsIS, subPointKey,ierr)
 !write(*,*) subPointKey
 !call ISRestoreIndicesF90(subpointsIS, subPointKey,ierr)
 call ISView(subpointsIS,PETSC_VIEWER_STDOUT_WORLD,ierr)

 [1]PETSC ERROR: - Error Message
 --
 [1]PETSC ERROR: Arguments must have same communicators
 [1]PETSC ERROR: Different communicators in the two objects: Argument #
 1 and 2 flag 3
 [1]PETSC ERROR: See https://petsc.org/release/faq/ for trouble
 shooting.
 [1]PETSC ERROR: Petsc Development GIT revision:
 v3.18.1-320-g7810d690132  GIT Date: 2022-11-20 20:25:41 -0600
 [1]PETSC ERROR: Configure options with-fc=mpiifort
 with-mpi-f90=mpiifort --download-triangle --download-parmetis
 --download-metis --with-debugging=1 --download-hdf5
 --prefix=/home/narnoldm/packages/petsc_install
 [1]PETSC