from:"Blaise Bourdin"

Re: [petsc-users] DMPlex and Gmsh

2023-11-08 Thread Blaise Bourdin




Hi,


I think that you need to use the magical keyword “-dm_plex_gmsh_mark_vertices” for that


Blaise



On Nov 8, 2023, at 1:13 PM, Sharan Roongta  wrote:







Caution: External email. 








Dear Petsc team,

 

I want to load a .msh file generated using Gmsh software into the DMPlex object. There are several things I would want to clarify, but I would like to start with “Physical tags”.

 

If I have defined “Physical Points”, “Physical Surface”, and “Physical Volume” in my .geo file, I get the physical tags in the “.msh” file.
When I load this mesh in DMPlex, and view the DM:

 

call DMView(globalMesh, PETSC_VIEWER_STDOUT_WORLD,err_PETSc)

  CHKERRQ(err_PETSc)

 

This is the output I get:

 

DM Object: n/a 1 MPI process

  type: plex

n/a in 3 dimensions:

  Number of 0-cells per rank: 14

  Number of 1-cells per rank: 49

  Number of 2-cells per rank: 60

  Number of 3-cells per rank: 24

Labels:

  celltype: 4 strata with value/size (0 (14), 6 (24), 3 (60), 1 (49))

  depth: 4 strata with value/size (0 (14), 1 (49), 2 (60), 3 (24))

  Cell Sets: 1 strata with value/size (8 (24))

  Face Sets: 6 strata with value/size (2 (4), 3 (4), 4 (4), 5 (4), 6 (4), 7 (4))

I was expecting to get the “Node Sets” or “Vertex Sets” also. Is my assumption wrong?

If yes, then how can one figure out the boundary nodes and their tags where I want to apply certain boundary conditions?
Currently we apply boundary conditions on faces, therefore “Face Sets” was enough. But now we want to apply displacements on certain boundary nodes.

 

I have also attached the .geo and .msh file (hope you can open it)

The Petsc version I am using is 3.18.6.


Thanks and Regards,

Sharan Roongta

 

---
Max-Planck-Institut für Eisenforschung GmbH 
Max-Planck-Straße 1 
D-40237 Düsseldorf 

Handelsregister B 2533 
Amtsgericht Düsseldorf 
   
Geschäftsführung 
Prof. Dr. Gerhard Dehm 
Prof. Dr. Jörg Neugebauer 
Prof. Dr. Dierk Raabe 
Dr. Kai de Weldige 

Ust.-Id.-Nr.: DE 11 93 58 514 
Steuernummer: 105 5891 1000 
- 
Please consider that invitations and e-mails of our institute are 
only valid if they end with …@mpie.de. 
If you are not sure of the validity please contact r...@mpie.de

Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails 
aus unserem Haus nur mit der Endung …@mpie.de gültig sind. 
In Zweifelsfällen wenden Sie sich bitte an r...@mpie.de

 




-
Stay up to date and follow us on LinkedIn, Twitter and YouTube.

Max-Planck-Institut für Eisenforschung GmbH
Max-Planck-Straße 1
D-40237 Düsseldorf
 
Handelsregister B 2533 
Amtsgericht Düsseldorf
 
Geschäftsführung
Prof. Dr. Gerhard Dehm
Prof. Dr. Jörg Neugebauer
Prof. Dr. Dierk Raabe
Dr. Kai de Weldige
 
Ust.-Id.-Nr.: DE 11 93 58 514 
Steuernummer: 105 5891 1000


Please consider that invitations and e-mails of our institute are 
only valid if they end with …@mpie.de. 
If you are not sure of the validity please contact r...@mpie.de

Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails
aus unserem Haus nur mit der Endung …@mpie.de gültig sind. 
In Zweifelsfällen wenden Sie sich bitte an r...@mpie.de
-





















— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] Face Set not appearing when loading GMSH mesh to DMPlex

2023-09-21 Thread Blaise Bourdin




Can you send me the geo or msh file, or instruction on how to build your example? I am not familiar with the gmsh API


Blaise



On Sep 20, 2023, at 8:16 AM, Anna Dalklint via petsc-users  wrote:







Caution: External email. 









Hello,


I am trying to load a mesh, created via the gmsh api, into a DMPlex. The code is seen below:


  6     gmsh::initialize();

  7     elsize0 = 1.0;
  8 double a = 20.0;   
  9 double b = 20.0;   
 10 double c = 10.0;   
 11 double d = 30.0;   
 12 double e = 70.0;   
 13     
 14 gmsh::model::geo::addPoint(0.0, 0.0, 0.0, elsize0, 1);
 15 gmsh::model::geo::addPoint(a, 0.0, 0.0, elsize0, 2);
 16 gmsh::model::geo::addPoint(0.0, d, 0.0, elsize0, 3);
 17 gmsh::model::geo::addPoint(a, c, 0.0, elsize0, 4);
 18 gmsh::model::geo::addPoint(e, d, 0.0, elsize0, 5);
 19 gmsh::model::geo::addPoint(e, c, 0.0, elsize0, 6);
 20 gmsh::model::geo::addPoint(a, 0.0, b, elsize0, 7);
 21 gmsh::model::geo::addPoint(0.0, 0.0, b, elsize0, 8);
 22 gmsh::model::geo::addPoint(0.0, d, b, elsize0, 9);
 23 gmsh::model::geo::addPoint(e, d, b, elsize0, 10);
 24 gmsh::model::geo::addPoint(e, c, b, elsize0, 11);
 25 gmsh::model::geo::addPoint(a, c, b, elsize0, 12);
 26 
 27 // Define lines connecting points
 28 gmsh::model::geo::addLine(2, 1, 1);
 29 gmsh::model::geo::addLine(1, 3, 2);
 30 gmsh::model::geo::addLine(3, 5, 3);
 31 gmsh::model::geo::addLine(5, 6, 4);
 32 gmsh::model::geo::addLine(6, 4, 5);
 33 gmsh::model::geo::addLine(4, 2, 6);
 34 gmsh::model::geo::addLine(2, 7, 7);
 35 gmsh::model::geo::addLine(7, 8, 8);
 36 gmsh::model::geo::addLine(8, 1, 9);
 37 gmsh::model::geo::addLine(8, 9, 10);
 38 gmsh::model::geo::addLine(9, 3, 11);
 39 gmsh::model::geo::addLine(9, 10, 12);
 40 gmsh::model::geo::addLine(10, 11, 13);
 41 gmsh::model::geo::addLine(11, 12, 14);
 42 gmsh::model::geo::addLine(12, 7, 15);
 43 gmsh::model::geo::addLine(12, 4, 16);
 44 gmsh::model::geo::addLine(11, 6, 17);
 45 gmsh::model::geo::addLine(10, 5, 18);
 46 
 47 // Define surfaces through their boundaries
 48 gmsh::model::geo::addCurveLoop({3, 4, 5, 6, 1, 2}, 1);
 49 gmsh::model::geo::addPlaneSurface({1}, 1);
 50 gmsh::model::geo::addCurveLoop({18, 4, -17, -13}, 2);
 51 gmsh::model::geo::addPlaneSurface({2}, 2);
 52 gmsh::model::geo::addCurveLoop({5, -16, -14, 17}, 3);
 53 gmsh::model::geo::addPlaneSurface({3}, 3);
 54 gmsh::model::geo::addCurveLoop({16, 6, 7, -15}, 4);
 55 gmsh::model::geo::addPlaneSurface({4}, 4);
 56 gmsh::model::geo::addCurveLoop({12, 13, 14, 15, 8, 10}, 5);
 57 gmsh::model::geo::addPlaneSurface({5}, 5);
 58 gmsh::model::geo::addCurveLoop({11, -2, -9, 10}, 6);
 59 gmsh::model::geo::addPlaneSurface({6}, 6);
 60 gmsh::model::geo::addCurveLoop({1, -9, -8, -7}, 7);
 61 gmsh::model::geo::addPlaneSurface({7}, 7);
 62 gmsh::model::geo::addCurveLoop({11, 3, -18, -12}, 8); 
 63 gmsh::model::geo::addPlaneSurface({8}, 8); 
 64  
 65 gmsh::model::geo::addSurfaceLoop({8, 1, 2, 3, 4, 7, 6, 5}, 1);
 66 gmsh::model::geo::addVolume({1}, 1);
 67 
 68 gmsh::model::geo::synchronize();
 70 
 72 gmsh::model::addPhysicalGroup(3, {1}, 1);
 74 gmsh::model::addPhysicalGroup(2, {7}, 101, "dir_y");


  98 gmsh::option::setNumber("Mesh.Algorithm", 6);
  99     gmsh::option::setNumber("Mesh.ElementOrder", 1);
  100   gmsh::model::mesh::generate(3);

  101   gmsh::write("filename.msh");




The mesh is generated, but the Physical Group is not correctly loaded into the DMPlex. If I run a DMView on my dm, I get:



171 DM Object: Parallel Mesh 4 MPI processes
172   type: plex
173 Parallel Mesh in 3 dimensions:
174   Number of 0-cells per rank: 7138 7192 7328 7171
175   Number of 1-cells per rank: 45825 45992 46817 45824
176   Number of 2-cells per rank: 74791 74901 76203 74568
177   Number of 3-cells per rank: 36103 36100 36713 35914
178 Labels:
179   depth: 4 strata with value/size (0 (7138), 1 (45825), 2 (74791), 3 (36103))
180   celltype: 4 strata with value/size (0 (7138), 1 (45825), 3 (74791), 6 (36103))
181   Cell Sets: 1 strata with value/size (1 (36103))
182   Face Sets: 0 strata with value/size ()
183   Vertex Sets: 1 strata with value/size (1 (5200))



Why are the elements of the physicalGroup "dir_y" (Plane Surface 7) not loaded as a face set? If I change the Plane Surface to another surface, e.g. 8, it marks the elements as the correct face set. Another
 weird thing is that the number of vertices in the above Vertex set is not correct. 



I have tried changing the numbering of Plane Surface 7 (to e.g. 70) but the issue remains.



I run the code with the -dm_plex_gmsh_mark_vertices and -dm_plex_gmsh_multiple_tags flags.

Re: [petsc-users] PETSc with Xcode 15

2023-09-21 Thread Blaise Bourdin




FWIW, CLT 15.0 also seems to include changes to the linker, with incompatible options etc… I was able to rebuild mpich and petsc but I get many linker warnings and have not fully tested my build


Before CLT 15.0 update


SiMini:mef90-dmplex (dmplex)$ ld -v

@(#)PROGRAM:ld  PROJECT:ld64-857.1

BUILD 23:13:29 May  7 2023

configured to support archs: armv6 armv7 armv7s arm64 arm64e arm64_32 i386 x86_64 x86_64h armv6m armv7k armv7m armv7em

LTO support using: LLVM version 14.0.3, (clang-1403.0.22.14.1) (static support for 29, runtime is 29)



TAPI support using: Apple TAPI version 14.0.3 (tapi-1403.0.5.1)


vs after:

bblaptop:mef90-dmplex (dmplex)$ ld -v

@(#)PROGRAM:ld  PROJECT:dyld-1015.7

BUILD 18:48:48 Aug 22 2023

configured to support archs: armv6 armv7 armv7s arm64 arm64e arm64_32 i386 x86_64 x86_64h armv6m armv7k armv7m armv7em

will use ld-classic for: armv6 armv7 armv7s arm64_32 i386 armv6m armv7k armv7m armv7em

LTO support using: LLVM version 15.0.0 (static support for 29, runtime is 29)

TAPI support using: Apple TAPI version 15.0.0 (tapi-1500.0.12.3)

Library search paths:

Framework search paths:






bblaptop:tutorials (main)$ pwd

/opt/HPC/petsc-main/src/vec/vec/tutorials

bblaptop:tutorials (main)$ make ex1

mpicc -Wl,-bind_at_load -Wl,-multiply_defined,suppress -Wl,-multiply_defined -Wl,suppress -Wl,-search_paths_first -Wl,-no_compact_unwind -Wimplicit-function-declaration -Wunused -Wuninitialized -fPIC
 -g3 -O0  -I/opt/HPC/petsc-main/include -I/opt/HPC/petsc-main/ventura-gcc13.2-arm64-g/include -I/opt/X11/include      ex1.c  -Wl,-rpath,/opt/HPC/petsc-main/ventura-gcc13.2-arm64-g/lib -L/opt/HPC/petsc-main/ventura-gcc13.2-arm64-g/lib -Wl,-rpath,/opt/HPC/petsc-main/ventura-gcc13.2-arm64-g/lib
 -L/opt/HPC/petsc-main/ventura-gcc13.2-arm64-g/lib -Wl,-rpath,/opt/X11/lib -L/opt/X11/lib -Wl,-rpath,/opt/homebrew/Cellar/mpich/4.1.2/lib -L/opt/homebrew/Cellar/mpich/4.1.2/lib -Wl,-rpath,/opt/homebrew/Cellar/gcc/13.2.0/lib/gcc/current/gcc/aarch64-apple-darwin22/13
 -L/opt/homebrew/Cellar/gcc/13.2.0/lib/gcc/current/gcc/aarch64-apple-darwin22/13 -Wl,-rpath,/opt/homebrew/Cellar/gcc/13.2.0/lib/gcc/current/gcc -L/opt/homebrew/Cellar/gcc/13.2.0/lib/gcc/current/gcc -Wl,-rpath,/opt/homebrew/Cellar/gcc/13.2.0/lib/gcc/current
 -L/opt/homebrew/Cellar/gcc/13.2.0/lib/gcc/current -lpetsc -lHYPRE -ldmumps -lmumps_common -lpord -lpthread -lscalapack -lsuperlu -lml -llapack -lblas -lparmetis -lmetis -lexoIIv2for32 -lexodus -lnetcdf -lpnetcdf -lhdf5_hl -lhdf5 -lchaco -ltriangle -lz -lctetgen
 -lX11 -lmpifort -lmpi -lpmpi -lgfortran -lemutls_w -lquadmath -lc++ -o ex1

ld: warning: -multiply_defined is obsolete

ld: warning: -multiply_defined is obsolete

ld: warning: duplicate -rpath '/opt/HPC/petsc-main/ventura-gcc13.2-arm64-g/lib' ignored

ld: warning: -bind_at_load is deprecated on macOS

ld: warning: ignoring duplicate libraries: '-lmpi', '-lpmpi'









That’s quite a curveball a week ahead of a major software update.


Blaise



On Sep 20, 2023, at 10:43 PM, Ju Liu  wrote:








Caution: External email. 






Hi PETSc team:


I recently got my Xcode command line tools upgraded to version 15. When installing PETSc, the configure command returns an error. My
 configure command is 




./configure --with-cc=gcc --with-fc=0 --with-cxx=0 --download-f2cblaslapack --download-mpich




and the message is "Cannot locate
 all the standard C header files needed by PETSc".




The configure.log file is attached.




How shall I fix this? 




Thanks,




Ju
























— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

[petsc-users] MatSetSizes: C vs python

2023-06-06 Thread Blaise Bourdin





Hi,

Does anybody understand why MatSetSizes seem to behave differently in C and python?

I would expect the attached examples to be strictly equivalent but the python version fails in parallel. It may be that the python interface is different, but I don’t see any mention of this in the python docs.

Regards,
Blaise


SiBookPro:test (master)$ mpirun -np 2 python3 testL2G2.py  nrowsLoc: 10 ncolsLoc: 20
Traceback (most recent call last):
  File "/Users/blaise/Development/ccG_CR/test/testL2G2.py", line 20, in 
nrowsLoc: 11 ncolsLoc: 21
Traceback (most recent call last):
  File "/Users/blaise/Development/ccG_CR/test/testL2G2.py", line 20, in 
    sys.exit(main())
    sys.exit(main())
 ^^
  File "/Users/blaise/Development/ccG_CR/test/testL2G2.py", line 12, in main
 ^^
  File "/Users/blaise/Development/ccG_CR/test/testL2G2.py", line 12, in main
    P.setSizes([nrowsLoc,ncolsLoc],1)
    P.setSizes([nrowsLoc,ncolsLoc],1)
  File "petsc4py/PETSc/Mat.pyx", line 323, in petsc4py.PETSc.Mat.setSizes
petsc4py.PETSc.Error: error code 62
[1] MatSetSizes() at /opt/HPC/petsc-release/src/mat/utils/gcreate.c:161
[1] Invalid argument
[1] Int value must be same on all processes, argument # 4
  File "petsc4py/PETSc/Mat.pyx", line 323, in petsc4py.PETSc.Mat.setSizes
petsc4py.PETSc.Error: error code 62
[0] MatSetSizes() at /opt/HPC/petsc-release/src/mat/utils/gcreate.c:161
[0] Invalid argument
[0] Int value must be same on all processes, argument # 4







— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243







testL2G2.c
Description: testL2G2.c
#!/usr/bin/env python3
import sys
from petsc4py import PETSc

def main():
comm = PETSc.COMM_WORLD
nrowsLoc = 10 + comm.rank
ncolsLoc = 20 + comm.rank

P = PETSc.Mat().create(comm=comm)
print("nrowsLoc: {} ncolsLoc: {}".format(nrowsLoc,ncolsLoc))
P.setSizes([nrowsLoc,ncolsLoc],1)
P.setType('mpiaij')
P.setUp()
P.assemblyBegin()
P.assemblyEnd()
P.view()

if __name__ == "__main__":
sys.exit(main())

Re: [petsc-users] Composing different Field Components into single Vector (or data array)

2023-04-17 Thread Blaise Bourdin

Hi,

You can simply is VecISCopy.
Have a look at the construction of isU and friends at ex26.c:291, and calls to VecISCopy around line 357.

Blaise

On Apr 17, 2023, at 6:54 PM, James Wright wrote:

Thanks for the reply! The `Vec`s are created through DMs, not sections, but DMPlex uses sections, so I think this still applies. Super/Sub `DM`s sounds like a good solution. The only thing that I'm
not seeing in the examples is being able to take `Vec` data from one a sub `DM`'s `Vec` to a super `DM`'s `Vec`.

To be a bit more explicit, let's say I have 3 fields, A, B, C. I'm imagining creating a "main" `DM` with those fields and then creating sub `DM`s, dmA, dmB, dmC. I need to have, in a single `Vec` fields B and C, A and C together. So I'd create two super DMs
dm_AC and dm_BC.

If I have a global `Vec` on dm_BC, is it possible that I can get the C data onto a dm_AC `Vec`? I think I could using SubVectors and VecISCopy, but am not entirely sure. There doesn't appear to be an example in the tutorials anywhere that puts these pieces
together in this way.

Thanks,

James Wright
Graduate Research Assistant, PhD
University of Colorado Boulder

Cell: (864) 498 8869
Email:
ja...@jameswright.xyz
Website:
jameswright.xyz

On Mon, Apr 17, 2023 at 12:42 PM Blaise Bourdin <bour...@mcmaster.ca> wrote:

Hi,

If you have created your vectors using sections, all you need is to call DMCreateSuperDM. See src/dm/impls/plex/tests/ex26.c:300 for an example.

Blaise

On Apr 17, 2023, at 1:30 PM, James Wright <jrwrigh@gmail.com> wrote:

Hello,

I currently have two DMPlex objects, both `DMClone`ed from the same "original" DMPlex object. They have different fields with different numbers of components on each of them. I would like to compose certain components from each into a single contiguous array.
I'd also like to do this from a DMGlobal vector. What is the best/most robust way to do that?

Obviously I can just `VecGetArrayRead` each of the vectors and loop through them manually, but that relies on knowledge of the array data layout. There's `DMPlexGetLocalOffsets`, but since I want to use the global Vec, there isn't a corresponding `DMPlexGetGlobalOffsets`
(that I can see anyways). This manually looping would also require that the underlying arrangement of the DOFs is the same between the two DMPlex objects, but I assume this is true from the `DMClone` operation.

Thanks,

James Wright
Graduate Research Assistant, PhD
University of Colorado Boulder

Cell: (864) 498 8869
Email:
ja...@jameswright.xyz
Website:
jameswright.xyz

—
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] Composing different Field Components into single Vector (or data array)

2023-04-17 Thread Blaise Bourdin

Hi,

If you have created your vectors using sections, all you need is to call DMCreateSuperDM. See src/dm/impls/plex/tests/ex26.c:300 for an example.

Blaise

On Apr 17, 2023, at 1:30 PM, James Wright wrote:

Hello,

Thanks,

James Wright
Graduate Research Assistant, PhD
University of Colorado Boulder

Cell: (864) 498 8869
Email:
ja...@jameswright.xyz
Website:
jameswright.xyz

Re: [petsc-users] GAMG failure

2023-03-28 Thread Blaise Bourdin

On Mar 27, 2023, at 9:11 PM, Mark Adams wrote:

Yes, the eigen estimates are converging slowly.

BTW, have you tried hypre? It is a good solver (lots lots more woman years)
These eigen estimates are conceptually simple, but they can lead to problems like this (hypre and an eigen estimate free smoother).

I just moved from petsc 3.3 to main, so my experience with an old version of hyper has not been very convincing. Strangely enough, ML has always been the most efficient PC for me. Maybe it’s time to revisit.
That said, I would really like to get decent performances out of gamg. One day, I’d like to be able to account for the special structure of phase-field fracture in the construction of the coarse space.

But try this (good to have options anyway):

-pc_gamg_esteig_ksp_max_it 20

Chevy will scale the estimate that we give by, I think, 5% by default. Maybe 10.
You can set that with:

-mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.05

0.2 is the scaling of the high eigen estimate for the low eigen value in Chebyshev.

Jed’s suggestion of using -pc_gamg_reuse_interpolation 0 worked. I am testing your options at the moment.

Thanks a lot,

Blaise

—

Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] GAMG failure

2023-03-27 Thread Blaise Bourdin







On Mar 24, 2023, at 3:21 PM, Mark Adams  wrote:



* Do you set:


    PetscCall(MatSetOption(Amat, MAT_SPD, PETSC_TRUE));








    PetscCall(MatSetOption(Amat, MAT_SPD_ETERNAL, PETSC_TRUE));







Yes







Do that to get CG Eigen estimates. Outright failure is usually caused by a bad Eigen estimate.
-pc_gamg_esteig_ksp_monitor_singular_value

Will print out the estimates as its iterating. You can look at that to check that the max has converged.






I just did, and something is off:
I do multiple calls to SNESSolve (staggered scheme for phase-field fracture), but only get informations on the first solve (which is not the one failing, of course)
Here is what I get:


Residual norms for Displacement_pc_gamg_esteig_ solve.
  0 KSP Residual norm 7.636421712860e+01 % max 1.e+00 min 1.e+00 max/min 1.e+00
  1 KSP Residual norm 3.402024867977e+01 % max 1.114319928921e+00 min 1.114319928921e+00 max/min 1.e+00
  2 KSP Residual norm 2.124815079671e+01 % max 1.501143586520e+00 min 5.739351119078e-01 max/min 2.615528402732e+00
  3 KSP Residual norm 1.581785698912e+01 % max 1.644351137983e+00 min 3.263683482596e-01 max/min 5.038329074347e+00
  4 KSP Residual norm 1.254871990315e+01 % max 1.714668863819e+00 min 2.044075812142e-01 max/min 8.388479789416e+00
  5 KSP Residual norm 1.051198229090e+01 % max 1.760078533063e+00 min 1.409327403114e-01 max/min 1.248878386367e+01
  6 KSP Residual norm 9.061658306086e+00 % max 1.792995287686e+00 min 1.023484740555e-01 max/min 1.751853463603e+01
  7 KSP Residual norm 8.015529297567e+00 % max 1.821497535985e+00 min 7.818018001928e-02 max/min 2.329871248104e+01
  8 KSP Residual norm 7.201063258957e+00 % max 1.855140071935e+00 min 6.178572472468e-02 max/min 3.002538337458e+01
  9 KSP Residual norm 6.548491711695e+00 % max 1.903578294573e+00 min 5.008612895206e-02 max/min 3.800609738466e+01
 10 KSP Residual norm 6.002109992255e+00 % max 1.961356890125e+00 min 4.130572033722e-02 max/min 4.748390475004e+01
  Residual norms for Displacement_pc_gamg_esteig_ solve.
  0 KSP Residual norm 2.373573910237e+02 % max 1.e+00 min 1.e+00 max/min 1.e+00
  1 KSP Residual norm 8.845061415709e+01 % max 1.081192207576e+00 min 1.081192207576e+00 max/min 1.e+00
  2 KSP Residual norm 5.607525485152e+01 % max 1.345947059840e+00 min 5.768825326129e-01 max/min 2.333138869267e+00
  3 KSP Residual norm 4.123522550864e+01 % max 1.481153523075e+00 min 3.070603564913e-01 max/min 4.823655974348e+00
  4 KSP Residual norm 3.345765664017e+01 % max 1.551374710727e+00 min 1.953487694959e-01 max/min 7.941563771968e+00
  5 KSP Residual norm 2.859712984893e+01 % max 1.604588395452e+00 min 1.313871480574e-01 max/min 1.221267391199e+01
  6 KSP Residual norm 2.525636054248e+01 % max 1.650487481750e+00 min 9.322735730688e-02 max/min 1.770389646804e+01
  7 KSP Residual norm 2.270711391451e+01 % max 1.697243639599e+00 min 6.945419058256e-02 max/min 2.443687883140e+01
  8 KSP Residual norm 2.074739485241e+01 % max 1.737293728907e+00 min 5.319942519758e-02 max/min 3.265624999621e+01
  9 KSP Residual norm 1.912808268870e+01 % max 1.771708608618e+00 min 4.229776586667e-02 max/min 4.188657656771e+01
 10 KSP Residual norm 1.787394414641e+01 % max 1.802834420843e+00 min 3.460455235448e-02 max/min 5.209818645753e+01
  Residual norms for Displacement_pc_gamg_esteig_ solve.
  0 KSP Residual norm 1.361990679391e+03 % max 1.e+00 min 1.e+00 max/min 1.e+00
  1 KSP Residual norm 5.377188333825e+02 % max 1.086812916769e+00 min 1.086812916769e+00 max/min 1.e+00
  2 KSP Residual norm 2.819790765047e+02 % max 1.474233179517e+00 min 6.475176340551e-01 max/min 2.276745994212e+00
  3 KSP Residual norm 1.856720658591e+02 % max 1.646049713883e+00 min 4.391851040105e-01 max/min 3.747963441500e+00
  4 KSP Residual norm 1.446507859917e+02 % max 1.760403013135e+00 min 2.972886103795e-01 max/min 5.921528614526e+00
  5 KSP Residual norm 1.212491636433e+02 % max 1.839250080524e+00 min 1.921591413785e-01 max/min 9.571494061277e+00
  6 KSP Residual norm 1.052783637696e+02 % max 1.887062042760e+00 min 1.275920366984e-01 max/min 1.478981048966e+01
  7 KSP Residual norm 9.230292625762e+01 % max 1.917891358356e+00 min 8.853577120467e-02 max/min 2.166233300122e+01
  8 KSP Residual norm 8.262607594297e+01 % max 1.935857204308e+00 min 6.706949937710e-02 max/min 2.886345093206e+01
  9 KSP Residual norm 7.616474911000e+01 % max 1.946323901431e+00 min 5.354310733090e-02 max/min 3.635059671458e+01
 10 KSP Residual norm 7.138356892221e+01 % max 1.954382723686e+00 min 4.367661484659e-02 max/min 4.474666204216e+01
  Residual norms for Displacement_pc_gamg_esteig_ solve.
  0 KSP Residual norm 3.702300162209e+03 % max 1.e+00 min 1.e+00 max/min 1.e+00
  1 KSP Residual norm 1.255008322497e+03 % max

[petsc-users] GAMG failure

2023-03-24 Thread Blaise Bourdin

Hi,

I am having issue with GAMG for some very ill-conditioned 2D linearized 
elasticity problems (sharp variation of elastic moduli with thin  regions of 
nearly incompressible material). I use snes_type newtonls, linesearch_type cp, 
and pc_type gamg without any further options. pc_type Jacobi converges fine 
(although slowly of course).


I am not really surprised that gamg would not converge out of the box, but 
don’t know where to start to investigate the convergence failure. Can anybody 
help?

Blaise

— 
Canada Research Chair in Mathematical and Computational Aspects of Solid 
Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] dmplex overlap questions

2023-02-16 Thread Blaise Bourdin







On Feb 16, 2023, at 10:54 AM, Lawrence Mitchell  wrote:


Hi Blaise,

On Thu, 16 Feb 2023 at 15:17, Blaise Bourdin  wrote:

Hi,

I am trying to implement a non-local finite elements reconstruction operator in parallel.

Given a dmplex distributed with an overlap, is there a way to figure out which cells are in the overlap and which are not?


Yes. Get the pointSF of the DM, and the cell chart

DMPlexGetPointSF(dm, );
DMPlexGetHeightStratum(dm, 0, , );

Now get the graph (specifically ilocal of the sf):

PetscSFGetGraph(sf, NULL, , ,  NULL);

Now any value of ilocal that lies in [cstart, cend) is a cell which is
not owned by this process (i.e. in the overlap). Note that ilocal can
be NULL which just means it is the identity map [0, ..., nleaves), so
you just intersect [cstart, cend) with [0, nleaves) in that case to
find the overlap cells.

But that is very unlikely to be true, so:

for (PetscInt i = 0; i < nleaves; i++) {
   if (cstart <= ilocal[i] && ilocal[i] < cend) {
  // i is an overlap cell
   }
}
Alternatively, suppose that I distribute the same DM with and without an overlap. Is there any warranty that the distributions are compatible (i.e. coincide when the overlap is ignored)? If this is the case, can I assume that the non-overlap
 cells are numbered first in the overlapped dm?


If you do:

DMPlexDistribute(dm, 0, , );
DMPlexDistributeOverlap(paralleldm, 1, , );

Then paralleldm and overlapdm will be compatible, and I think it is
still the case that the overlap cells are numbered last (and
contiguously).

If you just do DMPlexDistribute(dm, 1, , ) then
you obviously don't have the non-overlapped one to compare, but it is
in this case still true that the overlap cells are numbered last.







Excellent. That is what I was hoping.


Regards,
Blaise






Thanks,

Lawrence





















— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

[petsc-users] dmplex overlap questions

2023-02-16 Thread Blaise Bourdin

Hi,

I am trying to implement a non-local finite elements reconstruction operator in 
parallel.

Given a dmplex distributed with an overlap, is there a way to figure out which 
cells are in the overlap and which are not?
Alternatively, suppose that I distribute the same DM with and without an 
overlap. Is there any warranty that the distributions are compatible (i.e. 
coincide when the overlap is ignored)? If this is the case, can I assume that 
the non-overlap cells are numbered first in the overlapped dm?


Regards,
Blaise


— 
Canada Research Chair in Mathematical and Computational Aspects of Solid 
Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] PETSc Fortran 64-bit

2023-02-01 Thread Blaise Bourdin







On Feb 1, 2023, at 3:50 PM, Satish Balay via petsc-users  wrote:




call DMPlexGetDepthStratum(dm, 0, vst, vend, ierr);CHKERRA(ierr)  gives an
error regarding the datatype of ierr. The error basically leads:
Error: Type mismatch in argument ‘b’ at (1); passed INTEGER(4) to
INTEGER(8)



Ok - I think the error here is with '0' - not ierr.

Use:







PetscInt zero
zero = 0
call DMPlexGetDepthStratum(dm, zero, vst, vend, ierr);CHKERRA(ierr)
<<<

Satish







Hence my suggestion to systematically use kind when passing constants!


Blaise






On Wed, 1 Feb 2023, Satish Balay via petsc-users wrote:

On Wed, 1 Feb 2023, Mike Michell wrote:

@Satish and Blaise - Thank you for the notes.

@Satish - When you say: "Some routines probably don't have Interface
definitions - hence compiler can't cat this issue with all functions."


sorry - meant to say 'catch' [not cat]

Does it mean that I cannot use 64-bit option for those routines, which do
not have the interface definitions?


All routines are usable. Its just that compiler treats routines
without interface definitions as F77 - and it won't verify the
data-types of arguments passed in. [i.e F77 mode..]

But I do see interface defs for both DMPlexGetDepthStratum() and
DMPlexGetChart() so don't know why they behave differently for you.

src/dm/f90-mod/ftn-auto-interfaces/petscdmplex.h90











 subroutine DMPlexGetDepthStratum(a,b,c,d,z)
  import tDM
  DM a ! DM
  PetscInt b ! PetscInt
  PetscInt c ! PetscInt
  PetscInt d ! PetscInt
  PetscErrorCode z
  end subroutine DMPlexGetDepthStratum

 subroutine DMPlexGetChart(a,b,c,z)
  import tDM
  DM a ! DM
  PetscInt b ! PetscInt
  PetscInt c ! PetscInt
  PetscErrorCode z
  end subroutine DMPlexGetChart

Satish


Thanks,
Mike


I use the following I all my fortran codes (inspired by a post from
Michael Metcalf on comp.lang.fortran many many moons ago):

PetscReal,Parameter :: PReal = 1.0
Integer,Parameter,Public :: Kr = Selected_Real_Kind(Precision(PReal))
PetscInt,Parameter :: PInt = 1
Integer,Parameter,Public :: Ki = kind(PInt)

You will need to pass constant with their kind (i.e. 1_Ki instead of 1)


The advantage of this approach over explicitly trying to figure out the
proper type of integer ois that the same code will automatically work with
32 and 64 bit indices.

I’ve been wondering if petsc should include these definitions (perhaps
under another name).

Blaise


On Feb 1, 2023, at 2:58 PM, Mike Michell  wrote:

Hi all,

I want to use PETSc with 64-bit indices with a Fortran90 code. It seems
some PETSc functions have no problem, but some of the others do not accept
32-bit error code integer (e.g., "ierr" declared as PetscErrorCode type).

For example,

call DMPlexGetChart(dm, pst, pend, ierr);CHKERRA(ierr)    works okay,

but

call DMPlexGetDepthStratum(dm, 0, vst, vend, ierr);CHKERRA(ierr)  gives an
error regarding the datatype of ierr. The error basically leads:
Error: Type mismatch in argument ‘b’ at (1); passed INTEGER(4) to
INTEGER(8)

I tried to declare ierr as integer(kind=8) type, but there are still some
problems. If PETSc is configured with 32-bit indices, the Fortran code
works without problem.

What surprising to me is that as mentioned above, DMPlexGetChart() works
okay, but  DMPlexGetDepthStratum() does not work with "ierr (PetscErrorCode
type)" variable with 64-bit indices.

Can I get any comments on it?

Thanks,
Mike


—
Canada Research Chair in Mathematical and Computational Aspects of Solid
Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243




























— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] PETSc Fortran 64-bit

2023-02-01 Thread Blaise Bourdin




I use the following I all my fortran codes (inspired by a post from Michael Metcalf on comp.lang.fortran many many moons ago):




PetscReal,Parameter :: PReal = 
1.0
Integer,Parameter,Public :: Kr = Selected_Real_Kind(Precision(PReal))

PetscInt,Parameter :: PInt = 
1
Integer,Parameter,Public :: Ki = kind(PInt)




You will need to pass constant with their kind (i.e. 1_Ki instead of 1)




The advantage of this approach over explicitly trying to figure out the proper type of integer ois that the same code will automatically work with 32 and 64 bit indices.


I’ve been wondering if petsc should include these definitions (perhaps under another name).


Blaise
 

On Feb 1, 2023, at 2:58 PM, Mike Michell  wrote:


Hi all,


I want to use PETSc with 64-bit indices with a Fortran90 code. It seems some PETSc functions have no problem, but some of the others do not accept 32-bit error code integer (e.g., "ierr" declared as PetscErrorCode type). 


For example, 


call DMPlexGetChart(dm, pst, pend, ierr);CHKERRA(ierr)    works okay, 


but



call DMPlexGetDepthStratum(dm, 0, vst, vend, ierr);CHKERRA(ierr)  gives an error regarding the datatype of ierr. The error basically leads: 

Error: Type mismatch in argument ‘b’ at (1); passed INTEGER(4) to INTEGER(8)



I tried to declare ierr as integer(kind=8) type, but there are still some problems. If PETSc is configured with 32-bit indices, the Fortran code works without problem. 


What surprising to me is that as mentioned above, DMPlexGetChart() works okay, but  DMPlexGetDepthStratum() does not work with "ierr (PetscErrorCode type)" variable with 64-bit indices. 


Can I get any comments on it? 


Thanks,
Mike





















— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] about repeat of expensive functions using VecScatterCreateToAll

2023-01-17 Thread Blaise Bourdin

Got it. Can you partition your mesh with only one processor in the z-direction? (Trivial if using DMDA)
Blaise

On Jan 17, 2023, at 4:49 PM, Venugopal, Vysakh (venugovh) wrote:

This is the support structure minimization filter. So I need to go layer-by-layer from the bottommost slice of the array and update it
as I move up. Every slice needs the updated values below that slice.
Vysakh

From: Blaise Bourdin <bour...@mcmaster.ca>
Sent: Tuesday, January 17, 2023 4:47 PM
To: Venugopal, Vysakh (venugovh) <venug...@mail.uc.edu>
Cc: Barry Smith <bsm...@petsc.dev>; petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] about repeat of expensive functions using VecScatterCreateToAll

External Email: Use Caution

What type of filter are you implementing?

Convolution filters are expensive to parallelize since you need an overlap of the size of the support of the filter, but it may still not be worst than doing it sequentially (typically
the filter size is only one or 2 element diameters). Or you may be able to apply the filter in Fourier space.

PDE-filters are typically elliptic and can be parallelized.

Blaise

On Jan 17, 2023, at 4:38 PM, Venugopal,
Vysakh (venugovh) via petsc-users <petsc-users@mcs.anl.gov> wrote:

Thank you! I am doing a structural optimization filter that inherently cannot be parallelized.

Vysakh

From: Barry Smith <bsm...@petsc.dev>
Sent: Tuesday, January 17, 2023 3:28 PM
To: Venugopal, Vysakh (venugovh) <venug...@mail.uc.edu>
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] about repeat of expensive functions using VecScatterCreateToAll

External Email: Use Caution

On Jan 17, 2023, at 3:12 PM, Venugopal, Vysakh (venugovh) via petsc-users <petsc-users@mcs.anl.gov>
wrote:

Hi,

I am doing the following thing.

Step 1. Create DM object and get global vector ‘V’ using DMGetGlobalVector.

Step 2. Doing some parallel operations on V.

Step 3. I am using VecScatterCreateToAll on V to create a sequential vector ‘V_SEQ’ using VecScatterBegin/End with SCATTER_FORWARD.

Step 4. I am performing an expensive operation on V_SEQ and outputting the updated V_SEQ.

Step 5. I am using VecScatterBegin/End with SCATTER_REVERSE (global and sequential flipped) to get V that is updated with new values from
V_SEQ.

Step 6. I continue using this new V on the rest of the parallelized program.

Question: Suppose I have n MPI processes, is the expensive operation in Step 4 repeated n times? If yes, is there a workaround such that
the operation in Step 4 is performed only once? I would like to follow the same structure as steps 1 to 6 with step 4 only performed once.

Each MPI rank is doing the same operations on its copy of the sequential vector. Since they are running in parallel it probably does not matter much that each is doing the same computation.
Step 5 does not require any MPI since only part of the sequential vector (which everyone has) is needed in the parallel vector.

You could use VecScatterCreateToZero() but then step 3 would require less communication but step 5 would require communication to get parts of the solution from rank 0 to the other
ranks. The time for step 4 would be roughly the same.

You will likely only see a worthwhile improvement in performance if you can parallelize the computation in 4. What are you doing that is computational intense and requires all the
data on a rank?

Barry

Thanks,

Vysakh Venugopal

---

Vysakh Venugopal

Ph.D. Candidate

Department of Mechanical Engineering

University of Cincinnati, Cincinnati, OH 45221-0072

—
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)

Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] about repeat of expensive functions using VecScatterCreateToAll

2023-01-17 Thread Blaise Bourdin

What type of filter are you implementing?
Convolution filters are expensive to parallelize since you need an overlap of the size of the support of the filter, but it may still not be worst than doing it sequentially (typically the filter size is only one or 2 element diameters). Or you may be
able to apply the filter in Fourier space.
PDE-filters are typically elliptic and can be parallelized.

Blaise

On Jan 17, 2023, at 4:38 PM, Venugopal,
Vysakh (venugovh) via petsc-users wrote:

Thank you! I am doing a structural optimization filter that inherently cannot be parallelized.
Vysakh

From: Barry Smith
Sent: Tuesday, January 17, 2023 3:28 PM
To: Venugopal, Vysakh (venugovh)
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] about repeat of expensive functions using VecScatterCreateToAll

External Email: Use Caution

On Jan 17, 2023, at 3:12 PM, Venugopal, Vysakh (venugovh) via petsc-users
wrote:

Hi,

I am doing the following thing.

Step 1. Create DM object and get global vector ‘V’ using DMGetGlobalVector.

Step 2. Doing some parallel operations on V.

Step 3. I am using VecScatterCreateToAll on V to create a sequential vector ‘V_SEQ’ using VecScatterBegin/End with SCATTER_FORWARD.

Step 4. I am performing an expensive operation on V_SEQ and outputting the updated V_SEQ.

Step 5. I am using VecScatterBegin/End with SCATTER_REVERSE (global and sequential flipped) to get V that is updated with new values from
V_SEQ.

Step 6. I continue using this new V on the rest of the parallelized program.

You will likely only see a worthwhile improvement in performance if you can parallelize the computation in 4. What are you doing that is computational intense and requires all the
data on a rank?

Barry

Thanks,

Vysakh Venugopal

---

Vysakh Venugopal

Ph.D. Candidate

Department of Mechanical Engineering

University of Cincinnati, Cincinnati, OH 45221-0072

Re: [petsc-users] coordinate degrees of freedom for 2nd-order gmsh mesh

2023-01-12 Thread Blaise Bourdin




Out of curiosity, what is the rationale for _reading_ high order gmsh meshes?
Is it so that one can write data back in native gmsh format? 


Regards,
Blaise



On Jan 12, 2023, at 7:13 PM, Matthew Knepley  wrote:



On Thu, Jan 12, 2023 at 1:33 PM Jed Brown  wrote:



It's confusing, but this line makes high order simplices always read as discontinuous coordinate spaces. I would love if someone would revisit that, perhaps also using DMPlexSetIsoperiodicFaceSF(),


Perhaps as a switch, but there is no way I am getting rid of the current periodicity. As we have discussed before, breaking the topological relation is a non-starter for me.


It does look like higher order Gmsh does read as DG. We can just project that to CG for non-periodic stuff.


  Thanks,


    Matt



which should simplify the code and avoid the confusing cell coordinates pattern. Sadly, I don't have time to dive in.

https://gitlab.com/petsc/petsc/-/commit/066ea43f7f75752f012be6cd06b6107ebe84cc6d#3616cad8148970af5b97293c49492ff893e25b59_1552_1724

"Daniel R. Shapero"  writes:

> Sorry either your mail system or mine prevented me from attaching the file,
> so I put it on pastebin:
> https://pastebin.com/awFpc1Js
>
> On Wed, Jan 11, 2023 at 4:54 PM Matthew Knepley  wrote:
>
>> Can you send the .msh file? I still have not installed Gmsh :)
>>
>>   Thanks,
>>
>>      Matt
>>
>> On Wed, Jan 11, 2023 at 2:43 PM Daniel R. Shapero  wrote:
>>
>>> Hi all -- I'm trying to read in 2nd-order / piecewise quadratic meshes
>>> that are generated by gmsh and I don't understand how the coordinates are
>>> stored in the plex. I've been discussing this with Matt Knepley here
>>> 
>>> as it pertains to Firedrake but I think this is more an issue at the PETSc
>>> level.
>>>
>>> This code
>>> 
>>> uses gmsh to generate a 2nd-order mesh of the unit disk, read it into a
>>> DMPlex, print out the number of cells in each depth stratum, and finally
>>> print a view of the coordinate DM's section. The resulting mesh has 64
>>> triangles, 104 edges, and 41 vertices. For 2nd-order meshes, I'd expected
>>> there to be 2 degrees of freedom at each node and 2 at each edge. The
>>> output is:
>>>
>>> ```
>>> Depth strata: [(64, 105), (105, 209), (0, 64)]
>>>
>>> PetscSection Object: 1 MPI process
>>>   type not yet set
>>> 1 fields
>>>   field 0 with 2 components
>>> Process 0:
>>>   (   0) dim 12 offset   0
>>>   (   1) dim 12 offset  12
>>>   (   2) dim 12 offset  24
>>> ...
>>>   (  62) dim 12 offset 744
>>>   (  63) dim 12 offset 756
>>>   (  64) dim  0 offset 768
>>>   (  65) dim  0 offset 768
>>> ...
>>>   ( 207) dim  0 offset 768
>>>   ( 208) dim  0 offset 768
>>>   PetscSectionSym Object: 1 MPI process
>>>     type: label
>>>     Label 'depth'
>>>     Symmetry for stratum value 0 (0 dofs per point): no symmetries
>>>     Symmetry for stratum value 1 (0 dofs per point): no symmetries
>>>     Symmetry for stratum value 2 (12 dofs per point):
>>>       Orientation range: [-3, 3)
>>>     Symmetry for stratum value -1 (0 dofs per point): no symmetries
>>> ```
>>>
>>> The output suggests that there are 12 degrees of freedom in each
>>> triangle. That would mean the coordinate field is discontinuous across cell
>>> boundaries. Can someone explain what's going on? I tried reading the .msh
>>> file but it's totally inscrutable to me. I'm happy to RTFSC if someone
>>> points me in the right direction. Matt tells me that the coordinate field
>>> should only be discontinuous if the mesh is periodic, but this mesh
>>> shouldn't be periodic.
>>>
>>
>>
>> --
>> 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/
>> 
>>





-- 






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/



























— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A,

[petsc-users] GAMG and linearized elasticity

2022-12-13 Thread Blaise Bourdin





Hi,

I am getting close to finish porting a code from petsc 3.3 / sieve to main / dmplex, but am now encountering difficulties

I am reasonably sure that the Jacobian and residual are correct. The codes handle boundary conditions differently (MatZeroRowCols vs dmplex constraints) so it is not trivial to compare them. Running with snes_type ksponly pc_type Jacobi or hyper gives me the
 same results in roughly the same number of iterations.

In my old code, gamg would work out of the box. When using petsc-main, -pc_type gamg -pc_gamg_type agg works for _some_ problems using P1-Lagrange elements, but never for P2-Lagrange. The typical error message is in gamg_agg.txt

When using -pc_type classical, a problem where the KSP would converge in 47 iteration in 3.3 now takes 1400.  ksp_view_3.3.txt and ksp_view_main.txt show the output of -ksp_view for both versions. I don’t notice anything obvious.

Strangely, removing the call to PCSetCoordinates does not have any impact on the convergence.

I am sure that I am missing something, or not passing the right options. What’s a good starting point for 3D elasticity?
Regards,
Blaise






— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243





[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Computed maximum singular value as zero
[0]PETSC ERROR: WARNING! There are option(s) set that were not used! Could be 
the program crashed before they were used or a spelling mistake, etc!
[0]PETSC ERROR: Option left: name:-displacement_ksp_converged_reason value: 
ascii source: file
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.18.2-341-g16200351da0  GIT 
Date: 2022-12-12 23:42:20 +
[0]PETSC ERROR: 
/home/bourdinb/Development/mef90/mef90-dmplex/bbserv-gcc11.2.1-mvapich2-2.3.7-O/bin/ThermoElasticity
 on a bbserv-gcc11.2.1-mvapich2-2.3.7-O named bb01 by bourdinb Tue Dec 13 
17:02:19 2022
[0]PETSC ERROR: Configure options --CFLAGS=-Wunused 
--FFLAGS="-ffree-line-length-none -fallow-argument-mismatch -Wunused" 
--COPTFLAGS="-O2 -march=znver2" --CXXOPTFLAGS="-O2 -march=znver2" 
--FOPTFLAGS="-O2 -march=znver2" --download-chaco=1 --download-exodusii=1 
--download-fblaslapack=1 --download-hdf5=1 --download-hypre=1 
--download-metis=1 --download-ml=1 --download-mumps=1 --download-netcdf=1 
--download-p4est=1 --download-parmetis=1 --download-pnetcdf=1 
--download-scalapack=1 --download-sowing=1 
--download-sowing-cc=/opt/rh/devtoolset-9/root/usr/bin/gcc 
--download-sowing-cxx=/opt/rh/devtoolset-9/root/usr/bin/g++ 
--download-sowing-cpp=/opt/rh/devtoolset-9/root/usr/bin/cpp 
--download-sowing-cxxcpp=/opt/rh/devtoolset-9/root/usr/bin/cpp 
--download-superlu=1 --download-triangle=1 --download-yaml=1 --download-zlib=1 
--with-debugging=0 --with-mpi-dir=/opt/HPC/mvapich2/2.3.7-gcc11.2.1 --with-pic 
--with-shared-libraries=1 --with-mpiexec=srun --with-x11=0
[0]PETSC ERROR: #1 PCGAMGOptProlongator_AGG() at 
/1/HPC/petsc/main/src/ksp/pc/impls/gamg/agg.c:779
[0]PETSC ERROR: #2 PCSetUp_GAMG() at 
/1/HPC/petsc/main/src/ksp/pc/impls/gamg/gamg.c:639
[0]PETSC ERROR: #3 PCSetUp() at 
/1/HPC/petsc/main/src/ksp/pc/interface/precon.c:994
[0]PETSC ERROR: #4 KSPSetUp() at 
/1/HPC/petsc/main/src/ksp/ksp/interface/itfunc.c:405
[0]PETSC ERROR: #5 KSPSolve_Private() at 
/1/HPC/petsc/main/src/ksp/ksp/interface/itfunc.c:824
[0]PETSC ERROR: #6 KSPSolve() at 
/1/HPC/petsc/main/src/ksp/ksp/interface/itfunc.c:1070
[0]PETSC ERROR: #7 SNESSolve_KSPONLY() at 
/1/HPC/petsc/main/src/snes/impls/ksponly/ksponly.c:48
[0]PETSC ERROR: #8 SNESSolve() at 
/1/HPC/petsc/main/src/snes/interface/snes.c:4693
[0]PETSC ERROR: #9 
/home/bourdinb/Development/mef90/mef90-dmplex/ThermoElasticity/ThermoElasticity.F90:228
  Linear solve converged due to CONVERGED_RTOL iterations 46
KSP Object:(Disp_) 32 MPI processes
  type: cg
  maximum iterations=1
  tolerances:  relative=1e-05, absolute=1e-08, divergence=1e+10
  left preconditioning
  using nonzero initial guess
  using PRECONDITIONED norm type for convergence test
PC Object:(Disp_) 32 MPI processes
  type: gamg
MG: type is MULTIPLICATIVE, levels=4 cycles=v
  Cycles per PCApply=1
  Using Galerkin computed coarse grid matrices
  Coarse grid solver -- level ---
KSP Object:(Disp_mg_coarse_) 32 MPI processes
  type: gmres
GMRES: restart=30, using Classical (unmodified) Gram-Schmidt 
Orthogonalization with no iterative refinement
GMRES: happy breakdown tolerance 1e-30
  maximum iterations=1, initial

Re: [petsc-users] Union of sequential vecs

2022-12-09 Thread Blaise Bourdin

Hi,

If you use an IS to store your ints, you can also do
ISAllGather followed by ISSortRemoveDups

Blaise

On Dec 9, 2022, at 4:14 PM, Barry Smith wrote:

Ok, so you want the unique list of integers sorted from all the seq vectors on ever MPI rank?

VecScatterCreateToAll() to get all values on all ranks (make the sequential vectors MPI vectors instead).
create an integer array long enough to hold all of them
Use VecGetArray() and a for loop to copy all the values to the integer array,
Use PetscSortRemoveDupsInt on the integer array

Now each rank has all the desired values.

On Dec 9, 2022, at 3:24 PM, Karthikeyan Chockalingam - STFC UKRI wrote:

That is where I am stuck, I don’t know who to combine them
to get Vec = {2,5,7,8,10,11,12}.
I just want them in an MPI vector.
I finally plan to call VecScatterCreateToAll so that all processor gets a copy.
Thank you.
Kind regards,
Karthik.

From: Barry Smith
Date: Friday, 9 December 2022 at 20:04
To: Chockalingam, Karthikeyan (STFC,DL,HC)
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] Union of sequential vecs

How are you combining them to get Vec = {2,5,7,8,10,11,12}?

Do you want the values to remain on the same MPI rank as before, just in an MPI vector?

On Dec 9, 2022, at 2:28 PM, Karthikeyan Chockalingam - STFC UKRI via petsc-users
wrote:

Hi,

I want to take the union of a set of sequential vectors, each living in a different processor.

Say,

Vec_Seq1 = {2,5,7}

Vec_Seq2 = {5,8,10,11}

Vec_Seq3 = {5,2,12}.

Finally, get the union of all them Vec = {2,5,7,8,10,11,12}.

I initially wanted to create a parallel vector and insert the (sequential vector) values but I do not know, to which index to insert the values to. But
I do know the total size of Vec (which in this case is 7).

Any help is much appreciated.

Kind regards,

Karthik.

This email and any attachments are intended solely for the use of the named recipients. If you are not the intended recipient you must not use, disclose, copy or distribute this email or any of its attachments
and should notify the sender immediately and delete this email from your system. UK Research and Innovation (UKRI) has taken every reasonable precaution to minimise risk of this email or any attachments containing viruses or malware but the recipient should
carry out its own virus and malware checks before opening the attachments. UKRI does not accept any liability for any losses or damages which the recipient may sustain due to presence of any viruses.

[petsc-users] Postdoctoral openings at McMaster University

2022-11-28 Thread Blaise Bourdin

Dear PETSc users and developers,

I am looking for several postdocs for various projects on phase-field models 
including but not limited to:
   - Large scale numerical simulation of frack propagation in polycrystals and 
comparison with experiments
   - New algorithms and discretization schemes for phase-field fracture
   - Study of crack nucleation in nominally brittle materials

Candidates should have a background in applied mathematics or solid mechanics 
and familiarity with PETSc. Start dates are flexible but I need to fill at 
least one position immediately. Interested candidates should reach out to me by 
email (bour...@mcmaster.ca <mailto:bour...@mcmaster.ca>) and attach a vitae and 
list of publications.


Feel free to forward this posting.

Blaise Bourdin

— 
Canada Research Chair in Mathematical and Computational Aspects of Solid 
Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

2022-11-14 Thread Blaise Bourdin




Replying to myself so that it gets into the googles:


The reference tetrahedron used by DMPlexComputeCellGeometryAffineFEM has vertices at (-1,1,-1), (-1,-1,-1), (1, -1, -1) and (-1,-1,1).


Blaise
 


On Nov 10, 2022, at 6:42 PM, Matthew Knepley  wrote:



On Thu, Nov 10, 2022 at 3:46 PM Blaise Bourdin <bour...@mcmaster.ca> wrote:




I am not sure I am buying this… If the tet was inverted, detJ would be negative, but it is always 1/8, as expected.


The attached mesh is a perfectly valid tet generated by Cubit, with orientation matching the exodus documentation (ignore the mid-edge dof since this is a tet4).
Here is what I get out of the code I attached in my previous email:





Yes, I use the opposite convention from ExodusII. In my opinion, orienting face (1, 2, 3) to have an inward normal is sacrilegious.


  Thanks,


     Matt
 





SiMini:Tests (dmplex)$ ./TestDMPlexComputeCellGeometryAffineFEM -i ../TestMeshes/TetCubit.gen 

 filename ../TestMeshes/TetCubit.gen

Vec Object: coordinates 1 MPI process

  type: seq

0.

0.

0.

1.

0.

0.

0.

1.

0.

0.

0.

1.

v0

 0:   1.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

 0:   0.e+00   0.e+00   2.e+00

detJ : 0.125



From J, invJ, and v0, I still can’t reconstruct a reasonable reference tet which I was naively assuming was either the unit simplex, or the simplex with vertices at (-1,-1,-1), (-1,0,-1), (0, -1, -1), and (-1,-1,1) not necessarily in this order. In order
 to build my FE basis functions on the reference element, I really need to know what this element is…


Blaise






 







On Nov 9, 2022, at 6:56 PM, Matthew Knepley <knep...@gmail.com> wrote:



On Wed, Nov 9, 2022 at 10:46 AM Blaise Bourdin <bour...@mcmaster.ca> wrote:







On Nov 9, 2022, at 10:04 AM, Matthew Knepley <knep...@gmail.com> wrote:



On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin <bour...@mcmaster.ca> wrote:



Hi,

What reference simplex is DMPlexComputeCellGeometryAffineFEM using in 2 and 3D?
I am used to computing my shape functions on the unit simplex (vertices at the origin and each e_i), but it does not look to be the reference simplex in this function:

In 3D, for the unit simplex with vertices at (0,0,0) (1,0,0) (0,1,0) (0,0,1) (in this order), I get J = 1 / 2 . [[-1,-1,-1],[1,0,0],[0,0,1]] and v0 = [0,0,1]

In 2D, for the unit simplex with vertices at (0,0), (1,0), and (0,1), I get J = 1 / 2. I and v0 = [0,0], which does not make any sense to me (I was assuming that the 2D reference simplex had vertices at (-1,-1), (1, -1) and (-1,1), but if this were the case,
 v0 would not be 0).

I can build a simple example with meshes consisting only of the unit simplex in 2D and 3D if that would help.



I need to rewrite the documentation on geometry, but I was waiting until I rewrite the geometry calculations to fit into libCEED. Toby found a nice
way to express them in BLAS form which I need to push through everything.


I always think of operating on the cell with the first vertex at the origin (I think it is easier), so I have a xi0 that translates the first vertex
of the reference to the origin, and a v0 that translates the first vertex of the real cell to the origin. You can see this here


  https://gitlab.com/petsc/petsc/-/blob/main/include/petsc/private/petscfeimpl.h#L251


This explains the 2D result. I cannot understand your 3D result, unless the vertices are in another order.






That makes two of us, then… I am attaching a small example and test meshes (one cell being the unit simplex starting with the origin and numbered in direct order when looking from (1,1,1)





Oh, it is probably inverted. All faces are oriented for outward normals. It is in the Orientation chapter in the book :)


  Thanks,


      Matt
 





 filename ../TestMeshes/1Tri.gen

Vec Object: coordinates 1 MPI process

  type: seq

0.

0.

1.

0.

0.

1.

v0

 0:   0.e+00   0.e+00

J

 0:   5.e-01   0.e+00

 0:   0.e+00   5.e-01

invJ

 0:   2.e+00  -0.e+00

 0:  -0.e+00   2.e+00

detJ : 0.25


And 


 filename ../TestMeshes/1Tet.gen

Vec Object: coordinates 1 MPI process

  type: seq

0.

0.

0.

1.
0.


0.

0.

1.

0.

0.

0.

1.

v0

 0:   1.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

2022-11-10 Thread Blaise Bourdin




I am not sure I am buying this… If the tet was inverted, detJ would be negative, but it is always 1/8, as expected.


The attached mesh is a perfectly valid tet generated by Cubit, with orientation matching the exodus documentation (ignore the mid-edge dof since this is a tet4).
Here is what I get out of the code I attached in my previous email:




SiMini:Tests (dmplex)$ ./TestDMPlexComputeCellGeometryAffineFEM -i ../TestMeshes/TetCubit.gen 

 filename ../TestMeshes/TetCubit.gen

Vec Object: coordinates 1 MPI process

  type: seq

0.

0.

0.

1.

0.

0.

0.

1.

0.

0.

0.

1.

v0

 0:   1.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

 0:   0.e+00   0.e+00   2.e+00

detJ : 0.125



From J, invJ, and v0, I still can’t reconstruct a reasonable reference tet which I was naively assuming was either the unit simplex, or the simplex with vertices at (-1,-1,-1), (-1,0,-1), (0, -1, -1), and (-1,-1,1) not necessarily in this order. In order
 to build my FE basis functions on the reference element, I really need to know what this element is…


Blaise






 







On Nov 9, 2022, at 6:56 PM, Matthew Knepley  wrote:



On Wed, Nov 9, 2022 at 10:46 AM Blaise Bourdin <bour...@mcmaster.ca> wrote:







On Nov 9, 2022, at 10:04 AM, Matthew Knepley <knep...@gmail.com> wrote:



On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin <bour...@mcmaster.ca> wrote:



Hi,

What reference simplex is DMPlexComputeCellGeometryAffineFEM using in 2 and 3D?
I am used to computing my shape functions on the unit simplex (vertices at the origin and each e_i), but it does not look to be the reference simplex in this function:

In 3D, for the unit simplex with vertices at (0,0,0) (1,0,0) (0,1,0) (0,0,1) (in this order), I get J = 1 / 2 . [[-1,-1,-1],[1,0,0],[0,0,1]] and v0 = [0,0,1]

In 2D, for the unit simplex with vertices at (0,0), (1,0), and (0,1), I get J = 1 / 2. I and v0 = [0,0], which does not make any sense to me (I was assuming that the 2D reference simplex had vertices at (-1,-1), (1, -1) and (-1,1), but if this were the case,
 v0 would not be 0).

I can build a simple example with meshes consisting only of the unit simplex in 2D and 3D if that would help.



I need to rewrite the documentation on geometry, but I was waiting until I rewrite the geometry calculations to fit into libCEED. Toby found a nice
way to express them in BLAS form which I need to push through everything.


I always think of operating on the cell with the first vertex at the origin (I think it is easier), so I have a xi0 that translates the first vertex
of the reference to the origin, and a v0 that translates the first vertex of the real cell to the origin. You can see this here


  https://gitlab.com/petsc/petsc/-/blob/main/include/petsc/private/petscfeimpl.h#L251


This explains the 2D result. I cannot understand your 3D result, unless the vertices are in another order.






That makes two of us, then… I am attaching a small example and test meshes (one cell being the unit simplex starting with the origin and numbered in direct order when looking from (1,1,1)





Oh, it is probably inverted. All faces are oriented for outward normals. It is in the Orientation chapter in the book :)


  Thanks,


      Matt
 





 filename ../TestMeshes/1Tri.gen

Vec Object: coordinates 1 MPI process

  type: seq

0.

0.

1.

0.

0.

1.

v0

 0:   0.e+00   0.e+00

J

 0:   5.e-01   0.e+00

 0:   0.e+00   5.e-01

invJ

 0:   2.e+00  -0.e+00

 0:  -0.e+00   2.e+00

detJ : 0.25


And 


 filename ../TestMeshes/1Tet.gen

Vec Object: coordinates 1 MPI process

  type: seq

0.

0.

0.

1.
0.


0.

0.

1.

0.

0.

0.

1.

v0

 0:   1.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

 0:   0.e+00   0.e+00   2.e+00

detJ : 0.125




I don’t understand why v0=(0,0) in 2D and (1,0,0) in 3D (but don’t really care) since I only want J. J makes no sense to me in 3D. In particular, one does not seem to have X~ = invJ.X + v0 (X = J.(X~-v0) as stated in CoordinatesRefToReal (it works in 2D
 if V0 = (1,1), which is consistent with a reference simplex with vertices at (-1,-1), (1,-1) and (-1,1)).


What am I missing?


Blaise





/














  Thanks,


      Matt
 

Regards,
Blaise



—

Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

2022-11-09 Thread Blaise Bourdin







On Nov 9, 2022, at 10:04 AM, Matthew Knepley  wrote:



On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin <bour...@mcmaster.ca> wrote:



Hi,

What reference simplex is DMPlexComputeCellGeometryAffineFEM using in 2 and 3D?
I am used to computing my shape functions on the unit simplex (vertices at the origin and each e_i), but it does not look to be the reference simplex in this function:

In 3D, for the unit simplex with vertices at (0,0,0) (1,0,0) (0,1,0) (0,0,1) (in this order), I get J = 1 / 2 . [[-1,-1,-1],[1,0,0],[0,0,1]] and v0 = [0,0,1]

In 2D, for the unit simplex with vertices at (0,0), (1,0), and (0,1), I get J = 1 / 2. I and v0 = [0,0], which does not make any sense to me (I was assuming that the 2D reference simplex had vertices at (-1,-1), (1, -1) and (-1,1), but if this were the case,
 v0 would not be 0).

I can build a simple example with meshes consisting only of the unit simplex in 2D and 3D if that would help.



I need to rewrite the documentation on geometry, but I was waiting until I rewrite the geometry calculations to fit into libCEED. Toby found a nice
way to express them in BLAS form which I need to push through everything.


I always think of operating on the cell with the first vertex at the origin (I think it is easier), so I have a xi0 that translates the first vertex
of the reference to the origin, and a v0 that translates the first vertex of the real cell to the origin. You can see this here


  https://gitlab.com/petsc/petsc/-/blob/main/include/petsc/private/petscfeimpl.h#L251


This explains the 2D result. I cannot understand your 3D result, unless the vertices are in another order.






That makes two of us, then… I am attaching a small example and test meshes (one cell being the unit simplex starting with the origin and numbered in direct order when looking from (1,1,1)


 filename ../TestMeshes/1Tri.gen

Vec Object: coordinates 1 MPI process

  type: seq

0.

0.

1.

0.

0.

1.

v0

 0:   0.e+00   0.e+00

J

 0:   5.e-01   0.e+00

 0:   0.e+00   5.e-01

invJ

 0:   2.e+00  -0.e+00

 0:  -0.e+00   2.e+00

detJ : 0.25


And 


 filename ../TestMeshes/1Tet.gen

Vec Object: coordinates 1 MPI process

  type: seq

0.

0.

0.

1.
0.


0.

0.

1.

0.

0.

0.

1.

v0

 0:   1.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

 0:   0.e+00   0.e+00   2.e+00

detJ : 0.125




I don’t understand why v0=(0,0) in 2D and (1,0,0) in 3D (but don’t really care) since I only want J. J makes no sense to me in 3D. In particular, one does not seem to have X~ = invJ.X + v0 (X = J.(X~-v0) as stated in CoordinatesRefToReal (it works in 2D
 if V0 = (1,1), which is consistent with a reference simplex with vertices at (-1,-1), (1,-1) and (-1,1)).


What am I missing?


Blaise





/














  Thanks,


      Matt
 

Regards,
Blaise



— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243






-- 






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/



























— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243






















TestDMPlexComputeCellGeometryAffineFEM.c
Description: TestDMPlexComputeCellGeometryAffineFEM.c


1Tet.gen
Description: 1Tet.gen


1Tri.gen
Description: 1Tri.gen

[petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

2022-11-08 Thread Blaise Bourdin

Hi,

What reference simplex is DMPlexComputeCellGeometryAffineFEM using in 2 and 3D?
I am used to computing my shape functions on the unit simplex (vertices at the 
origin and each e_i), but it does not look to be the reference simplex in this 
function:

In 3D, for the unit simplex with vertices at (0,0,0) (1,0,0) (0,1,0) (0,0,1) 
(in this order), I get J = 1 / 2 . [[-1,-1,-1],[1,0,0],[0,0,1]] and v0 = [0,0,1]

In 2D, for the unit simplex with vertices at (0,0), (1,0), and (0,1), I get J = 
1 / 2. I and v0 = [0,0], which does not make any sense to me (I was assuming 
that the 2D reference simplex had vertices at (-1,-1), (1, -1) and (-1,1), but 
if this were the case, v0 would not be 0).

I can build a simple example with meshes consisting only of the unit simplex in 
2D and 3D if that would help.

Regards,
Blaise



— 
Canada Research Chair in Mathematical and Computational Aspects of Solid 
Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] Dof ordering in DMPlexVecGet/Set/RestoreClosure

2022-11-07 Thread Blaise Bourdin




Thanks.
I had figured out the ordering by stratum but was dreading reconstructing the ordering on test and hexes…


Blaise



On Nov 7, 2022, at 12:14 PM, Matthew Knepley  wrote:



On Mon, Nov 7, 2022 at 10:51 AM Blaise Bourdin <bour...@mcmaster.ca> wrote:



Hi,

How are degree of freedom ordered in calls to DMPlexVec[Set/Get/Restore]Closure?
Given a FE mesh and a section for P2-Lagrange elements (i.e. 1 dof per vertex and edge), I was naively assuming that the “local” vector would contain dof at vertices then edges (in this order), since it matches the section ordering, but it looks like I get
 edges then vertices…

Here is my section (my mesh has 8 cells, 16 edges, and 9 vertices)
PetscSection Object: U 1 MPI process
  type not yet set
1 fields
  field 0 with 1 components
Process 0:
  (   0) dim  0 offset   0
  (   1) dim  0 offset   0
  (   2) dim  0 offset   0
  (   3) dim  0 offset   0
  (   4) dim  0 offset   0
  (   5) dim  0 offset   0
  (   6) dim  0 offset   0
  (   7) dim  0 offset   0
  (   8) dim  1 offset   0
  (   9) dim  1 offset   1
  (  10) dim  1 offset   2
  (  11) dim  1 offset   3
  (  12) dim  1 offset   4
  (  13) dim  1 offset   5
  (  14) dim  1 offset   6
  (  15) dim  1 offset   7
  (  16) dim  1 offset   8
  (  17) dim  1 offset   9
  (  18) dim  1 offset  10
  (  19) dim  1 offset  11
  (  20) dim  1 offset  12
  (  21) dim  1 offset  13
  (  22) dim  1 offset  14
  (  23) dim  1 offset  15
  (  24) dim  1 offset  16
  (  25) dim  1 offset  17
  (  26) dim  1 offset  18
  (  27) dim  1 offset  19
  (  28) dim  1 offset  20
  (  29) dim  1 offset  21
  (  30) dim  1 offset  22
  (  31) dim  1 offset  23
  (  32) dim  1 offset  24

Start from the following local vector:
Vec Object: U 1 MPI process
  type: seq
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.

Call PetscCallA(DMPlexVecGetClosure(dmU,sectionU,U,0_Ki,UArray,ierr)) and get the following for Array:
10.000        11.000        12.000        1.        2.        3.

Is this ordering predictable and documented somewhere? Is it ordered by stratum?



The ordering is determined. It follows the same order that DMPlexGetTransitiveClosure() gives for the points.


Transitive closure orders points by stratum. It is a BFS of the Hasse Diagram, starting from the initial point. For
your triangle, it would be 


  tri0, e0, e1, e2, v0, v1, v2


For each cell type, the order of faces is specified in Table 1.2 of the attached. This gives an order to each level
of the BFS.


  Thanks,


     Matt
 

Regards,
Blaise

— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243






-- 






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/




























— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

[petsc-users] Dof ordering in DMPlexVecGet/Set/RestoreClosure

2022-11-07 Thread Blaise Bourdin

Hi,

How are degree of freedom ordered in calls to DMPlexVec[Set/Get/Restore]Closure?
Given a FE mesh and a section for P2-Lagrange elements (i.e. 1 dof per vertex 
and edge), I was naively assuming that the “local” vector would contain dof at 
vertices then edges (in this order), since it matches the section ordering, but 
it looks like I get edges then vertices…

Here is my section (my mesh has 8 cells, 16 edges, and 9 vertices)
PetscSection Object: U 1 MPI process
  type not yet set
1 fields
  field 0 with 1 components
Process 0:
  (   0) dim  0 offset   0
  (   1) dim  0 offset   0
  (   2) dim  0 offset   0
  (   3) dim  0 offset   0
  (   4) dim  0 offset   0
  (   5) dim  0 offset   0
  (   6) dim  0 offset   0
  (   7) dim  0 offset   0
  (   8) dim  1 offset   0
  (   9) dim  1 offset   1
  (  10) dim  1 offset   2
  (  11) dim  1 offset   3
  (  12) dim  1 offset   4
  (  13) dim  1 offset   5
  (  14) dim  1 offset   6
  (  15) dim  1 offset   7
  (  16) dim  1 offset   8
  (  17) dim  1 offset   9
  (  18) dim  1 offset  10
  (  19) dim  1 offset  11
  (  20) dim  1 offset  12
  (  21) dim  1 offset  13
  (  22) dim  1 offset  14
  (  23) dim  1 offset  15
  (  24) dim  1 offset  16
  (  25) dim  1 offset  17
  (  26) dim  1 offset  18
  (  27) dim  1 offset  19
  (  28) dim  1 offset  20
  (  29) dim  1 offset  21
  (  30) dim  1 offset  22
  (  31) dim  1 offset  23
  (  32) dim  1 offset  24

Start from the following local vector:
Vec Object: U 1 MPI process
  type: seq
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.

Call PetscCallA(DMPlexVecGetClosure(dmU,sectionU,U,0_Ki,UArray,ierr)) and get 
the following for Array:
10.00011.00012.000
1.2.3.

Is this ordering predictable and documented somewhere? Is it ordered by stratum?

Regards,
Blaise

— 
Canada Research Chair in Mathematical and Computational Aspects of Solid 
Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] Problem about optimized version

2022-08-31 Thread Blaise Bourdin




Hi,


Are you sure that you are using the same version of petsc in your batch script? Can you echo $PETSC_DIR and $PETSC_ARCH?


Also, it looks like you are using my hydraulic fracturing code VPFHF. That's great and I am glad to see that somebody is continuing to use it. We have not touched the code since a bit after Chukwudi graduated and you must have had to do a few
 changes to compile it with petsc 3.17. WOuld you mind contributing your changes back to the main repository? A pull request would be greatly appreciated.


Regards,
Blaise Bourdin


On Aug 31, 2022, at 9:39 AM, wangxq2020--- via petsc-users <petsc-users@mcs.anl.gov> wrote:

Hi
I run the program on the debuged version of PETSC successfuly. Then I want using the optimized version. And configure by ./configure --with-debugging=no --with-packages-download-dir=~/packages_dir
 --with-mpich-dir="/public3/soft/mpich/mpich-3.4.2". Then I installed the Hypre package by
./configure --download-hypre=/public3/home/scg6368/packages_dir/hypre-2.24.0.tar.gz --with-mpich-dir="/public3/soft/mpich/mpich-3.4.2"
. The result as below



=
                         Configuring PETSc to compile on your system                         
=
=                                            * WARNING: You have a version of GNU make older than 4.0. It will
 work,                                                        but may not support all the parallel testing options. You can install the                                                          latest GNU make with your package manager, such as brew or macports,
 or use                                                        the --download-make option to get the latest GNU make *                                                                  = 
                                     =                                            Running configure on HYPRE; this may take several minutes                         
                                           =                                      = 
                                           Running make on HYPRE; this may take several minutes                                                                         = 
                                     =                                            Running make install on HYPRE; this may take several minutes                       
                                          =                                      Compilers:                                                                         
                                                
  C Compiler:         /public3/soft/mpich/mpich-3.4.2/bin/mpicc  -fPIC -Wall -Wwrite-strings -Wno-unknown-pragmas -Wno-lto-type-mismatch -fstack-protector -fvisibility=hidden -g3 -O0  
    Version: gcc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-44)
  C++ Compiler:         /public3/soft/mpich/mpich-3.4.2/bin/mpicxx  -Wall -Wwrite-strings -Wno-strict-aliasing -Wno-unknown-pragmas -Wno-lto-type-mismatch -fstack-protector -fvisibility=hidden -g -O0  -std=gnu++11
 -fPIC 
    Version: g++ (GCC) 4.8.5 20150623 (Red Hat 4.8.5-44)
  Fortran Compiler:         /public3/soft/mpich/mpich-3.4.2/bin/mpif90  -fPIC -Wall -ffree-line-length-0 -Wno-lto-type-mismatch -Wno-unused-dummy-argument -g -O0  
    Version: GNU Fortran (GCC) 4.8.5 20150623 (Red Hat 4.8.5-44)
Linkers:
  Shared linker:   /public3/soft/mpich/mpich-3.4.2/bin/mpicc  -shared  -fPIC -Wall -Wwrite-strings -Wno-unknown-pragmas -Wno-lto-type-mismatch -fstack-protector -fvisibility=hidden -g3 -O0
  Dynamic linker:   /public3/soft/mpich/mpich-3.4.2/bin/mpicc  -shared  -fPIC -Wall -Wwrite-strings -Wno-unknown-pragmas -Wno-lto-type-mismatch -fstack-protector -fvisibility=hidden -g3 -O0
  Libraries linked against:   -lquadmath -lstdc++ -ldl 
BlasLapack:
  Library:  -llapack -lblas
  Unknown if this uses OpenMP (try export OMP_NUM_THREADS=<1-4> yourprogram -log_view) 
  uses 4 byte integers
MPI:
  Version:  3
  Includes: -I/public3/soft/mpich/mpich-3.4.2/include

Re: [petsc-users] DMPlex: a Mapping Array between Natural and Distributed Cell Index

2022-07-18 Thread Blaise Bourdin

Alexis noticed a problem with the natural SF when constraints are defined. He has a MR coming.
@Alexis: Could it be that what Bora sees is exactly what you fixed?

Blaise

On Jul 18, 2022, at 8:35 PM, Bora Jeong wrote:

Thank you for the reply.

I am actually putting 1dof on the vertex points.
In the attached sample code, from line 63 to line 99, 1dof on all vertex points is defined for the section. But my problem is, if I print out Natural SF using the viewer, I only have 10 vertex points, instead of 15. Here 15 is the size of the
global vector defined by the aforementioned section. The working example code was attached to the previous email with a small grid.

Best

On Mon, Jul 18, 2022 at 6:50 PM Matthew Knepley wrote:

On Mon, Jul 18, 2022 at 4:30 PM Bora Jeong wrote:

Actually, it seems that my problem is that I do not properly understand what Natural SF means.

What I want to achieve is to extract natural vertex index from Natural SF, and make sure this natural vertex index is the same with what is written in a mesh file from Gmsh. The natural vertex index defined by Gmsh is attached as a snapshot with sample code
and mesh. I want to reproduce this indexing based on the Natural SF that Petsc defines.

I am struggling with that because the number of entries defined in the Natural SF that I can see through PetscSFView() is different from the number of vertex points that each processor actually owns. For example, with two processors and with the attached "square.msh",
each processor might have 15 vertex points if I configure Petsc with ParMetis. However, the defined Natural SF for proc=0 only has 10 roots. I expected, for each processor, the Natural SF would show me 15 entries through PetscSFView(), if I define 1 degree
of freedom on vertex when the Petsc Section is defined.

When Petsc defines the Natural SF, are there any omitted components since they are trivial to save?

The naturalSF is mapping the field defined by the Section, not the mesh. You can get the effect you are asking for by putting one
degree of freedom (dof) on every vertex for the Section. How are you defining your Section now?

Thanks,

Matt

Best

On Mon, Jul 18, 2022 at 3:50 PM Matthew Knepley wrote:

On Mon, Jul 18, 2022 at 12:14 PM Bora Jeong wrote:

Thank you for the corrections. It works.

However, I still have a problem.
I tested zero stratum for DMGetStratumIS() with the default depth label, since the graph of vertex is needed. Then, PetscSFGetGraph() from the Natural SF is crashed. I checked that pBcCompIS contains proper set of integers via ISView().

On the other hand, it works okay if DMGetStratumIS() is with dim stratum, which is for cell stratum.

The problem is that if the dim stratum is used, the number of entry in ileaves(i)% pointer does not match with the actual number of vertex that each processor has. That is why I tried to put the zero stratum on DMGetStratumIS(), instead of the stratum for cell
(i.e., "dim") to see if any changes on the graph. Can I get comments?

I cannot understand your question. I am not sure what you are using the set of vertices or cell for.

Thanks,

Matt

Best

On Sun, Jul 17, 2022 at 9:59 AM Matthew Knepley wrote:

On Fri, Jul 15, 2022 at 7:05 PM Bora Jeong wrote:

I found that iroots() and ileaves() have size of nleaves and tested through the attached sample code with the grid previously shared. Then I get the results written in the attached monitor file.

It seems ileaves(i)%rank and ileaves(i)%index at line 127 of the attached code have garbage values, different from displayed values by PetscSFView() at line 113.

It is tough to get it why the garbage values are returned from PetscSFGetGraph(). Any comments will be very appreciated.

Unfortunately, Fortran is not very good at reporting declaration errors. The problem is that you did not include or use the Vec module. I have done this and your example runs for me. I have included the modified code.

Thanks,

Matt

Best

On Fri, Jul 15, 2022 at 3:53 PM Matthew Knepley wrote:

On Fri, Jul 15, 2022 at 2:25 PM Bora Jeong wrote:

Thank you for the comments. Updated sample code is attached here; DMPlexDistribute() is explicitly called and by defining a section before mesh distribution, a Natural SF was able to be defined as found from PetscSFView().

However, I am still struggling to call PetscSFGetGraph()
https://petsc.org/main/docs/manualpages/PetscSF/PetscSFGetGraph/
due to data type definition of ilocal and iremote. What is the proper size allocation for those two variables? Does this size for the number of

[petsc-users] List of points with dof>0 in a PetscSection

2022-06-10 Thread Blaise Bourdin

Hi,

Given a PetscSection, is there an easy way to get a list of point at which the 
number of dof is >0?
For instance, when projecting over a FE space, I’d rather do a loop over such 
points than do a loop over all points in a DM, get the number of dof, and test 
if it is >0.

Regards,
Blaise

-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] problems with petsc4py on M1 mac when configuring PETSc

2022-04-08 Thread Blaise Bourdin




Veronika,


There is something screwy going on with some hombrew gcc installs. Not sure why. 
In many case, creating a symbolic link libgcc_s.dylib -> libgcc_s.1.1.dylib in /opt/homebrew/opt/gcc/lib/gcc/11 fixes this


Regards,
Blaise


On Apr 8, 2022, at 5:41 AM, Veronika Ulanova <19veronik...@gmail.com> wrote:


I pulled the updated repo but now I get a different error also related to architecture:


ld: file not found: @rpath/libgcc_s.1.1.dylib for architecture arm64

clang: error: linker command failed with exit code 1 (use -v to see invocation)

error: command '/Users/veronikaulanova/Documents/documents/Oxford_work/petsc/arch-darwin-c-debug/bin/mpicc' failed with exit code 1

**ERROR*

Error building petsc4py.



make[2]: *** [petsc4pybuild] Error 1

**ERROR*

 
Error during compile, check arch-darwin-c-debug/lib/petsc/conf/make.log

 
Send it and arch-darwin-c-debug/lib/petsc/conf/configure.log to 
petsc-ma...@mcs.anl.gov






Any idea about this one? Thank you



On Fri, Apr 8, 2022 at 10:09 AM Veronika Ulanova <19veronik...@gmail.com> wrote:


Hello Satish,

By the repo you mean the main PETSc repository?

Thank you, Veronika


On Fri, Apr 8, 2022 at 9:49 AM Satish Balay  wrote:


A fix for M1 was recently pushed to the git repo.

Can you try using it [with either release or main branches] - and see if that works?

Satish



On Fri, 8 Apr 2022, Veronika Ulanova wrote:

> Hi,
> 
> I have been trying to configure the petsc4py package along with my PETSc
> installation for my research. I have searched for solutions on the internet
> but there was no luck, this is the error I get when I am trying to
> configure PETSC and use it to run examples of code that utilise the
> petsc4py package...
> 
> 
> 
>  File
> "/Users/veronikaulanova/Documents/documents/Oxford_work/petsc/arch-darwin-c-debug/lib/petsc4py/PETSc.py",
> line 3, in 
> 
>     PETSc = ImportPETSc(ARCH)
> 
>   File
> "/Users/veronikaulanova/Documents/documents/Oxford_work/petsc/arch-darwin-c-debug/lib/petsc4py/lib/__init__.py",
> line 29, in ImportPETSc
> 
>     return Import('petsc4py', 'PETSc', path, arch)
> 
>   File
> "/Users/veronikaulanova/Documents/documents/Oxford_work/petsc/arch-darwin-c-debug/lib/petsc4py/lib/__init__.py",
> line 73, in Import
> 
>     module = import_module(pkg, name, path, arch)
> 
>   File
> "/Users/veronikaulanova/Documents/documents/Oxford_work/petsc/arch-darwin-c-debug/lib/petsc4py/lib/__init__.py",
> line 58, in import_module
> 
>     with f: return imp.load_module(fullname, f, fn, info)
> 
>   File "/opt/homebrew/Cellar/python@3.9/3.9.10/Frameworks/Python.framework/Versions/3.9/lib/python3.9/imp.py",
> line 242, in load_module
> 
>     return load_dynamic(name, filename, file)
> 
>   File "/opt/homebrew/Cellar/python@3.9/3.9.10/Frameworks/Python.framework/Versions/3.9/lib/python3.9/imp.py",
> line 342, in load_dynamic
> 
>     return _load(spec)
> 
> ImportError:
> dlopen(/Users/veronikaulanova/Documents/documents/Oxford_work/petsc/arch-darwin-c-debug/lib/petsc4py/lib/arch-darwin-c-debug/
> 
PETSc.cpython-39-darwin.so , 2): no
> suitable image found.  Did find:
> 
> /Users/veronikaulanova/Documents/documents/Oxford_work/petsc/arch-darwin-c-debug/lib/petsc4py/lib/arch-darwin-c-debug/
> 
PETSc.cpython-39-darwin.so : mach-o,
> but wrong architecture
> 
> /Users/veronikaulanova/Documents/documents/Oxford_work/petsc/arch-darwin-c-debug/lib/petsc4py/lib/arch-darwin-c-debug/
> 
PETSc.cpython-39-darwin.so : stat()
> failed with errno=25
> 
> 
> As you can see I am using an M1 mac and therefore that causes the petsc4py
> package to give me an error (highlighted above), but do you have any idea
> of how I can solve it or why it is even there? I know that the new M1 has a
> different architecture to intel macs but I have no idea of where to
> start with fixing this or how to make it work. Any help would be much
> appreciated, thank you in advance.
> 
> Kind Regards, Veronika
> 




















-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243

Re: [petsc-users] Message length unit in log_view

2022-02-25 Thread Blaise Bourdin




Thanks,
Blaise



On Feb 25, 2022, at 10:30 AM, Junchao Zhang <junchao.zh...@gmail.com> wrote:






On Fri, Feb 25, 2022 at 9:20 AM Blaise Bourdin <bour...@mcmaster.ca> wrote:


Hi,

What is the unit for messages length in PetscLogView?
MPI Messages:         3.500e+00     1.000   3.500e+00  7.000e+00
MPI Message Lengths:  9.200e+01     1.000   2.629e+01  1.840e+02
MPI Reductions:       4.000e+01     1.000

Are these bytes, words, KB?

bytes. Yes, we need to make it clear.
 


Regards,
Blaise

-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243



















-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243

[petsc-users] Message length unit in log_view

2022-02-25 Thread Blaise Bourdin

Hi,

What is the unit for messages length in PetscLogView?
MPI Messages: 3.500e+00 1.000   3.500e+00  7.000e+00
MPI Message Lengths:  9.200e+01 1.000   2.629e+01  1.840e+02
MPI Reductions:   4.000e+01 1.000

Are these bytes, words, KB?

Regards,
Blaise

-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243

Re: [petsc-users] Output data using ExodusIIViewer

2021-12-01 Thread Blaise Bourdin




David,


Here is a modified example.
Exodus needs some additional work prior to saving fields. See the attached modified example.


Blaise








On Dec 1, 2021, at 1:54 PM, Blaise Bourdin <bour...@mcmaster.ca> wrote:


OK, let me have a look.
Blaise



On Nov 30, 2021, at 7:31 PM, David Andrs <and...@gmail.com> wrote:


I see. I added a "Cell Sets" label like so:

    DMCreateLabel(dm, "Cell Sets");
    DMLabel cs_label;
    DMGetLabel(dm, "Cell Sets", _label);
    DMLabelAddStratum(cs_label, 0);
    PetscInt idxs[] = { 1 };
    IS is;
    ISCreateGeneral(comm, 1, idxs, PETSC_COPY_VALUES, );
    DMLabelSetStratumIS(cs_label, 0, is);
    ISDestroy();

Note, that I have only a single element (Quad4) in the mesh and I am just trying to get this working, so I understand what needs to happen for larger meshes.

This got me past the segfault, but now I see:

[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Number of vertices 1 in dimension 2 has no ExodusII type

So, I assume I need to do something more to make it work.

David

On Nov 30 2021, at 11:39 AM, Blaise Bourdin <bour...@mcmaster.ca> wrote:

It looks like your DM cannot be saved in exodus format as such. The exodus format requires that all cells be part of a single block (defined by ‘Cell Set’ labels), and that the cell sets consists of sequentially numbered cells.
Can you see if that is enough? If not, I will go through your example


Blaise



On Nov 30, 2021, at 9:50 AM, David Andrs <and...@gmail.com> wrote:



Hello!


I am trying to store data into an ExodusII file using the ExodusIIViewer, but running into a segfault inside PETSc. Attached is a minimal example showing the problem. It can very much be that I am missing something obvious. However, if I change
 the code to VTKViewer I get the desired output file.


Machine: MacBook Pro 2019
OS version/type: Darwin notak.local 21.1.0 Darwin Kernel Version 21.1.0: Wed Oct 13 17:33:23 PDT 2021; root:xnu-8019.41.5~1/RELEASE_X86_64 x86_64
PETSc: Petsc Release Version 3.16.1, Nov 01, 2021
MPI: MPICH 3.4.2
Compiler: clang-12


Call stack (not sure how relevant that is since it is from opt version):

* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=1, address=0xc)
    frame #0: 0x000102303ba9 libpetsc.3.16.dylib`DMView_PlexExodusII(dm=, viewer=) at plexexodusii.c:457:45 [opt]
   454        else if (degree == 2) nodes[cs] = nodesHexP2;
   455      }
   456      /* Compute the number of cells not in the connectivity table */
-> 457      cellsNotInConnectivity -= nodes[cs][3]*csSize;
   458
   459      ierr = ISRestoreIndices(stratumIS, );CHKERRQ(ierr);
   460      ierr = ISDestroy();CHKERRQ(ierr);




With regards,



David Andrs


















-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243





























-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243





























-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243


















exo2.c
Description: exo2.c

Re: [petsc-users] Output data using ExodusIIViewer

2021-12-01 Thread Blaise Bourdin




OK, let me have a look.
Blaise



On Nov 30, 2021, at 7:31 PM, David Andrs <and...@gmail.com> wrote:


I see. I added a "Cell Sets" label like so:

    DMCreateLabel(dm, "Cell Sets");
    DMLabel cs_label;
    DMGetLabel(dm, "Cell Sets", _label);
    DMLabelAddStratum(cs_label, 0);
    PetscInt idxs[] = { 1 };
    IS is;
    ISCreateGeneral(comm, 1, idxs, PETSC_COPY_VALUES, );
    DMLabelSetStratumIS(cs_label, 0, is);
    ISDestroy();

Note, that I have only a single element (Quad4) in the mesh and I am just trying to get this working, so I understand what needs to happen for larger meshes.

This got me past the segfault, but now I see:

[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Number of vertices 1 in dimension 2 has no ExodusII type

So, I assume I need to do something more to make it work.

David

On Nov 30 2021, at 11:39 AM, Blaise Bourdin <bour...@mcmaster.ca> wrote:

It looks like your DM cannot be saved in exodus format as such. The exodus format requires that all cells be part of a single block (defined by ‘Cell Set’ labels), and that the cell sets consists of sequentially numbered cells.
Can you see if that is enough? If not, I will go through your example


Blaise



On Nov 30, 2021, at 9:50 AM, David Andrs <and...@gmail.com> wrote:



Hello!


I am trying to store data into an ExodusII file using the ExodusIIViewer, but running into a segfault inside PETSc. Attached is a minimal example showing the problem. It can very much be that I am missing something obvious. However, if I change
 the code to VTKViewer I get the desired output file.


Machine: MacBook Pro 2019
OS version/type: Darwin notak.local 21.1.0 Darwin Kernel Version 21.1.0: Wed Oct 13 17:33:23 PDT 2021; root:xnu-8019.41.5~1/RELEASE_X86_64 x86_64
PETSc: Petsc Release Version 3.16.1, Nov 01, 2021
MPI: MPICH 3.4.2
Compiler: clang-12


Call stack (not sure how relevant that is since it is from opt version):

* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=1, address=0xc)
    frame #0: 0x000102303ba9 libpetsc.3.16.dylib`DMView_PlexExodusII(dm=, viewer=) at plexexodusii.c:457:45 [opt]
   454        else if (degree == 2) nodes[cs] = nodesHexP2;
   455      }
   456      /* Compute the number of cells not in the connectivity table */
-> 457      cellsNotInConnectivity -= nodes[cs][3]*csSize;
   458
   459      ierr = ISRestoreIndices(stratumIS, );CHKERRQ(ierr);
   460      ierr = ISDestroy();CHKERRQ(ierr);




With regards,



David Andrs


















-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243





























-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243

Re: [petsc-users] Output data using ExodusIIViewer

2021-11-30 Thread Blaise Bourdin




It looks like your DM cannot be saved in exodus format as such. The exodus format requires that all cells be part of a single block (defined by ‘Cell Set’ labels), and that the cell sets consists of sequentially numbered cells.
Can you see if that is enough? If not, I will go through your example


Blaise



On Nov 30, 2021, at 9:50 AM, David Andrs  wrote:


Hello!


I am trying to store data into an ExodusII file using the ExodusIIViewer, but running into a segfault inside PETSc. Attached is a minimal example showing the problem. It can very much be that I am missing something obvious. However, if I change
 the code to VTKViewer I get the desired output file.


Machine: MacBook Pro 2019
OS version/type: Darwin notak.local 21.1.0 Darwin Kernel Version 21.1.0: Wed Oct 13 17:33:23 PDT 2021; root:xnu-8019.41.5~1/RELEASE_X86_64 x86_64
PETSc: Petsc Release Version 3.16.1, Nov 01, 2021 
MPI: MPICH 3.4.2
Compiler: clang-12


Call stack (not sure how relevant that is since it is from opt version):

* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=1, address=0xc)
    frame #0: 0x000102303ba9 libpetsc.3.16.dylib`DMView_PlexExodusII(dm=, viewer=) at plexexodusii.c:457:45 [opt]
   454        else if (degree == 2) nodes[cs] = nodesHexP2;
   455      }
   456      /* Compute the number of cells not in the connectivity table */
-> 457      cellsNotInConnectivity -= nodes[cs][3]*csSize;
   458 
   459      ierr = ISRestoreIndices(stratumIS, );CHKERRQ(ierr);
   460      ierr = ISDestroy();CHKERRQ(ierr);



With regards,



David Andrs
















-- 
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
Tel. +1 (905) 525 9140 ext. 27243

Re: [petsc-users] [petsc-announce] Question regarding updating PETSc Fortran examples to embrace post F77 constructs

2016-08-27 Thread Blaise Bourdin

FINALLY! Let's get rid of fortran77 free form in examples. I can't think of any 
reason to self inflict such a suffering. 

Are there ANY compiler around that people use and would not be able to process 
free form examples?I can see a point in keeping compatibility with fortran77 in 
petsc. It would make sense to keep a few old style pure f77 examples in the 
using-fortran section, but for the rest of the examples, using fixed form 
serves no purpose other than unexplicable bugs caused when a macro expands to 
more than 72 cols.

Going farther, but it is a really un gratifying job that nobody wants to do, it 
would make sense of having fortran77 bindings through iso_c_binding. That would 
allow better argument type checking and debugging (I.e. Inspecting the content 
of a petsc object from the debugger in a fortran program). Would that prevent 
F77 interoperability? I am not sure. 

Blaise

Sent from a mobile device

> On Aug 26, 2016, at 9:54 PM, Barry Smith  wrote:
> 
> 
>   PETSc users,
> 
>We've always been very conservative in PETSc to keep almost all our 
> Fortran examples in a format that works with classic FORTRAN 77 constructs: 
> fixed line format, (72 character limit) and no use of ; to separate 
> operations on the same line, etc. 
> 
>   Is it time to forgo these constructs and use more modern Fortran 
> conventions in all our examples?
> 
>Any feedback is appreciated
> 
>Barry
> 
> Note: it would continue to be possible to use PETSc in the FORTRAN 77 style, 
> this is just a question about updating the examples.
>

[petsc-users] LineSearch question

2012-08-15 Thread Blaise Bourdin

Hi,

I figured it out, and it was my fault...
I was not zeroing out the matrix before assembling the Jacobian. Once I did 
this, every linesearch converges and each SNES converge in 1 iteration, as 
expected.

Sorry for the confusion.

Blaise

On Aug 15, 2012, at 10:51 AM, Barry Smith bsmith at mcs.anl.gov wrote:

 
   Blaise,
 
We are confused. Can you run both the bt and l2 with all those options and 
 send ALL the output from each of the two runs.
 
 Thanks
 
Barry
 
 
 On Aug 14, 2012, at 6:24 PM, Blaise Bourdin bourdin at lsu.edu wrote:
 
 Hi,
 
 
 On Aug 14, 2012, at 5:58 PM, Barry Smith bsmith at mcs.anl.gov wrote:
 
 
  Blaise,
 
   You can run with -snes_linesearch_monitor   -info   
 -ksp_monitor_true_residual -ksp_converged_reason
 
 
 to get much more information about what is happening and why the line 
 search is failing.
 Focussing on the relevant part of the output, I get the following for the 
 first SNES step
 
 *** Using l2
Residual norms for temp_ solve.
0 KSP preconditioned resid norm 2.352873068990e+00 true resid norm 
 3.742138023215e+00 ||r(i)||/||b|| 1.e+00
 [0] KSPDefaultConverged(): user has provided nonzero initial guess, 
 computing 2-norm of preconditioned RHS
1 KSP preconditioned resid norm 7.175926524783e-01 true resid norm 
 8.026926174904e-01 ||r(i)||/||b|| 2.145010719836e-01
2 KSP preconditioned resid norm 4.099791012407e-01 true resid norm 
 6.219898727406e-01 ||r(i)||/||b|| 1.662124349455e-01
3 KSP preconditioned resid norm 2.769612150659e-01 true resid norm 
 4.622335508644e-01 ||r(i)||/||b|| 1.235212458752e-01
4 KSP preconditioned resid norm 8.991164116822e-02 true resid norm 
 1.938840972976e-01 ||r(i)||/||b|| 5.181104921701e-02
5 KSP preconditioned resid norm 1.369794733551e-02 true resid norm 
 2.867541652138e-02 ||r(i)||/||b|| 7.662843097578e-03
6 KSP preconditioned resid norm 3.522245138600e-03 true resid norm 
 5.452585588775e-03 ||r(i)||/||b|| 1.457077626466e-03
7 KSP preconditioned resid norm 1.117008651636e-03 true resid norm 
 1.551905826961e-03 ||r(i)||/||b|| 4.147110067382e-04
8 KSP preconditioned resid norm 5.008635361807e-04 true resid norm 
 6.226120116381e-04 ||r(i)||/||b|| 1.663786872038e-04
9 KSP preconditioned resid norm 2.079118910844e-04 true resid norm 
 3.430664466007e-04 ||r(i)||/||b|| 9.167658821571e-05
   10 KSP preconditioned resid norm 4.528126074206e-05 true resid norm 
 9.520866575330e-05 ||r(i)||/||b|| 2.544231804457e-05
   11 KSP preconditioned resid norm 8.441137224072e-06 true resid norm 
 1.519916221879e-05 ||r(i)||/||b|| 4.061625232553e-06
   12 KSP preconditioned resid norm 1.839317974157e-06 true resid norm 
 3.245208227421e-06 ||r(i)||/||b|| 8.672069836252e-07
   13 KSP preconditioned resid norm 4.346353491153e-07 true resid norm 
 7.198101954057e-07 ||r(i)||/||b|| 1.923526580100e-07
   14 KSP preconditioned resid norm 6.321413805477e-08 true resid norm 
 1.280486229700e-07 ||r(i)||/||b|| 3.421803850515e-08
   15 KSP preconditioned resid norm 9.029476674935e-09 true resid norm 
 2.193598397200e-08 ||r(i)||/||b|| 5.861885327562e-09
 [0] KSPDefaultConverged(): Linear solver has converged. Residual norm 
 9.029476674935e-09 is less than relative tolerance 1.e-08 times 
 initial right hand side norm 2.352873068990e+00 at iteration 15
  Linear solve converged due to CONVERGED_RTOL iterations 15
 [0] SNESSolve_LS(): iter=0, linear solve iterations=15
 [0] SNESLSCheckResidual_Private(): ||J^T(F-Ax)||/||F-AX|| 4.179425981164e+00 
 near zero implies inconsistent rhs
  Line search: lambdas = [1, 0.5, 0], fnorms = [2.1936e-08, 1.87107, 
 3.74214]
  Line search terminated: lambda = 1, fnorms = 2.19338e-08
 [0] SNESSolve_LS(): fnorm=3.7421380232151638e+00, 
 gnorm=2.1933796240714327e-08, ynorm=5.1838420564550498e+00, lssucceed=1
  1 SNES Function norm 2.193379624071e-08 
 [0] SNESDefaultConverged(): Converged due to function norm 
 2.193379624071e-08  3.742138023215e-08 (relative tolerance)
 SNESTemp converged in in1 iterations. SNESConvergedReason is3 
 
 *** Using bt
Residual norms for temp_ solve.
0 KSP preconditioned resid norm 1.176436534193e+00 true resid norm 
 3.742138022596e+00 ||r(i)||/||b|| 9.8344e-01
 [0] KSPDefaultConverged(): user has provided nonzero initial guess, 
 computing 2-norm of preconditioned RHS
1 KSP preconditioned resid norm 3.587963259712e-01 true resid norm 
 8.026926179905e-01 ||r(i)||/||b|| 2.145010721173e-01
2 KSP preconditioned resid norm 2.049895509618e-01 true resid norm 
 6.219898720314e-01 ||r(i)||/||b|| 1.662124347560e-01
3 KSP preconditioned resid norm 1.384806072424e-01 true resid norm 
 4.622335514699e-01 ||r(i)||/||b|| 1.235212460370e-01
4 KSP preconditioned resid norm 4.495582078268e-02 true resid norm 
 1.938840967382e-01 ||r(i)||/||b|| 5.181104906751e-02
5 KSP preconditioned resid norm 6.848973644691e-03 true resid norm 
 2.867541656135e-02 ||r(i)||/||b

[petsc-users] LineSearch question

2012-08-14 Thread Blaise Bourdin

HI,

I am trying to understand if the following behavior is normal / expected:

I am solving a quasi-static evolution where at each time step, SNESSolve is 
called. My validation problem is a _static_ problem with 2 time steps (i.e. 2 
successive calls to SNESSolve with the same operator, jacobian, and right hand 
side). Furthermore, the problem is linear, so that the Jacobian is constant. I 
even reset the solution vector to the same value at each time step.

In this setting, I am expecting that at each time step, each SNESSolve should 
converge in exactly one iteration no matter what linesearch / snes type I use. 
Indeed, when setting the linesearch to l2 or cp, this is what I get. However, 
for all other choices, the second call to SNESSolve fails to converge with a 
SNESConvergedReason of -5 or -6.

The relevant code is as follow:
  Call VecSet(solTemp,0.0_Kr,ierr);CHKERRQ(ierr)
  Call FormInitialGuess(snesTemp,solTemp,MEF90Ctx,ierr);CHKERRQ(ierr)
  Call VecCopy(solTemp,tmpvec,ierr)

  Call SNESSolve(snesTemp,PETSC_NULL_OBJECT,solTemp,ierr);CHKERRQ(ierr)

  Call VecCopy(tmpvec,soltemp,ierr)
  Call SNESSolve(snesTemp,PETSC_NULL_OBJECT,solTemp,ierr);CHKERRQ(ierr)


Is this expected? I tried to call SNESLineSearchReset before the second call to 
SNESSolve, but this does not seem to have any effect.

Blaise



Below is the sample output using cp as the linesearch type in which case I get 
the expected convergence:
Solving time step1: t= 1.0E+00 
  0 SNES Function norm 3.742138023215e+00 
  Line search terminated: lambda = 1, fnorms = 2.1936e-08
  1 SNES Function norm 2.193598339906e-08 
SNES Object:(temp_) 1 MPI processes
  type: ls
  maximum iterations=50, maximum function evaluations=1
  tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
  total number of linear solver iterations=15
  total number of function evaluations=3
  KSP Object:  (temp_)   1 MPI processes
type: cg
maximum iterations=1
tolerances:  relative=1e-08, absolute=1e-50, divergence=1
left preconditioning
using nonzero initial guess
using PRECONDITIONED norm type for convergence test
  PC Object:  (temp_)   1 MPI processes
type: icc
  0 levels of fill
  tolerance for zero pivot 2.22045e-14
  using Manteuffel shift
  matrix ordering: natural
  factor fill ratio given 1, needed 1
Factored matrix follows:
  Matrix Object:   1 MPI processes
type: seqsbaij
rows=104, cols=104
package used to perform factorization: petsc
total: nonzeros=381, allocated nonzeros=381
total number of mallocs used during MatSetValues calls =0
block size is 1
linear system matrix = precond matrix:
Matrix Object: 1 MPI processes
  type: seqaij
  rows=104, cols=104
  total: nonzeros=658, allocated nonzeros=658
  total number of mallocs used during MatSetValues calls =0
not using I-node routines
  SNESLineSearch Object:  (temp_)   1 MPI processes
type: cp
maxstep=1.00e+08, minlambda=1.00e-12
tolerances: relative=1.00e-08, absolute=1.00e-15, 
lambda=1.00e-08
maximum iterations=1
SNESTemp converged in in1 iterations. SNESConvergedReason is3 
  0 SNES Function norm 3.742138023215e+00 
  Line search: lambdas = [1, 0], ftys = [2.42597, 4.85193]
  Line search terminated: lambda = 2, fnorms = 2.1936e-08
  1 SNES Function norm 2.193598717772e-08 
SNES Object:(temp_) 1 MPI processes
  type: ls
  maximum iterations=50, maximum function evaluations=1
  tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
  total number of linear solver iterations=15
  total number of function evaluations=3
  KSP Object:  (temp_)   1 MPI processes
type: cg
maximum iterations=1
tolerances:  relative=1e-08, absolute=1e-50, divergence=1
left preconditioning
using nonzero initial guess
using PRECONDITIONED norm type for convergence test
  PC Object:  (temp_)   1 MPI processes
type: icc
  0 levels of fill
  tolerance for zero pivot 2.22045e-14
  using Manteuffel shift
  matrix ordering: natural
  factor fill ratio given 1, needed 1
Factored matrix follows:
  Matrix Object:   1 MPI processes
type: seqsbaij
rows=104, cols=104
package used to perform factorization: petsc
total: nonzeros=381, allocated nonzeros=381
total number of mallocs used during MatSetValues calls =0
block size is 1
linear system matrix = precond matrix:
Matrix Object: 1 MPI processes
  type: seqaij
  rows=104, cols=104
  total: nonzeros=658, allocated nonzeros=658
  total number of mallocs used during MatSetValues calls =0
not using I-node routines
  SNESLineSearch Object:  (temp_)   1 MPI processes
type: cp
maxstep=1.00e+08,

[petsc-users] LineSearch question

2012-08-14 Thread Blaise Bourdin

)||/||b|| 2.544231959867e-05
   11 KSP preconditioned resid norm 4.220568874641e-06 true resid norm 
1.519916304124e-05 ||r(i)||/||b|| 4.061625452335e-06
   12 KSP preconditioned resid norm 9.196589910150e-07 true resid norm 
3.245208482213e-06 ||r(i)||/||b|| 8.672070517123e-07
   13 KSP preconditioned resid norm 2.173176852660e-07 true resid norm 
7.198102176806e-07 ||r(i)||/||b|| 1.923526639624e-07
   14 KSP preconditioned resid norm 3.160707194729e-08 true resid norm 
1.280486368595e-07 ||r(i)||/||b|| 3.421804221680e-08
   15 KSP preconditioned resid norm 4.514738683754e-09 true resid norm 
2.193598711826e-08 ||r(i)||/||b|| 5.861886168328e-09
[0] KSPDefaultConverged(): Linear solver has converged. Residual norm 
4.514738683754e-09 is less than relative tolerance 1.e-08 times 
initial right hand side norm 1.176436534495e+00 at iteration 15
  Linear solve converged due to CONVERGED_RTOL iterations 15
[0] SNESSolve_LS(): iter=0, linear solve iterations=15
[0] SNESLSCheckResidual_Private(): ||J^T(F-Ax)||/||F-AX|| 8.358851664967e+00 
near zero implies inconsistent rhs
[0] VecScatterCreate(): Special case: sequential vector general scatter
[0] SNESLineSearchApply_BT(): Initial fnorm 3.742138023215e+00 gnorm 
1.871069011453e+00
  Line search: Using full step: fnorm 3.742138023215e+00 gnorm 
1.871069011453e+00
[0] SNESSolve_LS(): fnorm=3.7421380232151638e+00, gnorm=1.8710690114527022e+00, 
ynorm=2.5919210284812890e+00, lssucceed=1
  1 SNES Function norm 1.871069011453e+00 


As expected, the KSP residuals are exactly the same. I am not sure what to make 
of
[0] SNESLSCheckResidual_Private(): ||J^T(F-Ax)||/||F-AX|| 8.358851664967e+00 
near zero implies inconsistent rhs. 
Which RHS is this referring to? the RHS for SNESSolve is 0 (second argument of 
SNESSolve is PETSC_NULL_OBJECT). Could this mean that there is an 
incompatibility between my Jacobian and my Function?

In all case that diverge, it looks like the gnorm in the linesearch does not 
decrease.

 
  If that doesn't help you can send the code and we can play with it.

The code is a bit of a pain to build but both Matt and Jed have accounts on my 
systems. I can arrange to give access to others if necessary.

Blaise


 
Barry
 
 
 
 On Aug 14, 2012, at 5:53 PM, Blaise Bourdin bourdin at lsu.edu wrote:
 
 HI,
 
 I am trying to understand if the following behavior is normal / expected:
 
 I am solving a quasi-static evolution where at each time step, SNESSolve is 
 called. My validation problem is a _static_ problem with 2 time steps (i.e. 
 2 successive calls to SNESSolve with the same operator, jacobian, and right 
 hand side). Furthermore, the problem is linear, so that the Jacobian is 
 constant. I even reset the solution vector to the same value at each time 
 step.
 
 In this setting, I am expecting that at each time step, each SNESSolve 
 should converge in exactly one iteration no matter what linesearch / snes 
 type I use. Indeed, when setting the linesearch to l2 or cp, this is what I 
 get. However, for all other choices, the second call to SNESSolve fails to 
 converge with a SNESConvergedReason of -5 or -6.
 
 The relevant code is as follow:
 Call VecSet(solTemp,0.0_Kr,ierr);CHKERRQ(ierr)
 Call FormInitialGuess(snesTemp,solTemp,MEF90Ctx,ierr);CHKERRQ(ierr)
 Call VecCopy(solTemp,tmpvec,ierr)
 
 Call SNESSolve(snesTemp,PETSC_NULL_OBJECT,solTemp,ierr);CHKERRQ(ierr)
 
 Call VecCopy(tmpvec,soltemp,ierr)
 Call SNESSolve(snesTemp,PETSC_NULL_OBJECT,solTemp,ierr);CHKERRQ(ierr)
 
 
 Is this expected? I tried to call SNESLineSearchReset before the second call 
 to SNESSolve, but this does not seem to have any effect.
 
 Blaise
 
 
 
 Below is the sample output using cp as the linesearch type in which case I 
 get the expected convergence:
 Solving time step1: t= 1.0E+00 
 0 SNES Function norm 3.742138023215e+00 
 Line search terminated: lambda = 1, fnorms = 2.1936e-08
 1 SNES Function norm 2.193598339906e-08 
 SNES Object:(temp_) 1 MPI processes
 type: ls
 maximum iterations=50, maximum function evaluations=1
 tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
 total number of linear solver iterations=15
 total number of function evaluations=3
 KSP Object:  (temp_)   1 MPI processes
   type: cg
   maximum iterations=1
   tolerances:  relative=1e-08, absolute=1e-50, divergence=1
   left preconditioning
   using nonzero initial guess
   using PRECONDITIONED norm type for convergence test
 PC Object:  (temp_)   1 MPI processes
   type: icc
 0 levels of fill
 tolerance for zero pivot 2.22045e-14
 using Manteuffel shift
 matrix ordering: natural
 factor fill ratio given 1, needed 1
   Factored matrix follows:
 Matrix Object:   1 MPI processes
   type: seqsbaij
   rows=104, cols=104
   package used to perform factorization: petsc
   total: nonzeros=381, allocated nonzeros=381
   total number of mallocs used during

[petsc-users] pc_gamg_type geo broken in 3.3?

2012-07-25 Thread Blaise Bourdin

Sorry, I should have done that with the initial bug report. I 

Here is stack trace I get when using the intel 11.1 compilers. The segfault 
appears to be in triangle. WHen using gcc, everything seems to work fine.

Blaise

Program received signal EXC_BAD_ACCESS, Could not access memory.
Reason: KERN_INVALID_ADDRESS at address: 0x1402e000
0x00010f1e55aa in poolinit (pool=0x10fc31580, bytecount=32, itemcount=4092, 
wtype=FLOATINGPOINT, alignment=0) at src/triangle.c:3227
3227  *(pool-firstblock) = (VOID *) NULL;
(gdb) where
#0  0x00010f1e55aa in poolinit (pool=0x10fc31580, bytecount=32, 
itemcount=4092, wtype=FLOATINGPOINT, alignment=0) at src/triangle.c:3227
#1  0x00010f1e5da7 in initializepointpool () at src/triangle.c:3547
#2  0x00010f203462 in transfernodes (pointlist=0x7fb408aa2370, 
pointattriblist=0x0, pointmarkerlist=0x0, numberofpoints=2307, 
numberofpointattribs=0) at src/triangle.c:11486
#3  0x00010f205f69 in triangulate (triswitches=0x7fff6d9429c8 npczQ, 
in=0x7fff6d9429d0, out=0x7fff6d942828, vorout=0x0) at src/triangle.c:12965
#4  0x00010e46c325 in triangulateAndFormProl (selected_2=0x7fb408a9ad70, 
data_stride=12100, coords=0x113b0c770, nselected_1=2307, 
clid_lid_1=0x7fb408ac8f70, agg_lists_1=0x7fb408ab9b70, crsGID=0x7fb408b00770, 
bs=1, a_Prol=0x7fb408abb370, a_worst_best=0x7fff6d942dd0) at geo.c:253
#5  0x00010e47420d in PCGAMGProlongator_GEO (pc=0x7fb408a07370, 
Amat=0x7fb40889cb70, Gmat=0x7fb408b4d770, agg_lists=0x7fb408ab9b70, 
a_P_out=0x7fff6d944420) at geo.c:848
#6  0x00010e429526 in PCSetUp_GAMG (pc=0x7fb408a07370) at gamg.c:658
#7  0x00010ecc05d8 in PCSetUp (pc=0x7fb408a07370) at precon.c:832
#8  0x00010e5df6ed in KSPSetUp (ksp=0x7fb4089ccf70) at itfunc.c:278
#9  0x00010e5e12b3 in KSPSolve (ksp=0x7fb4089ccf70, b=0x7fb408968570, 
x=0x7fb40894cd70) at itfunc.c:402
#10 0x00010dd4c632 in main ()



On Jul 25, 2012, at 8:56 AM, Mark F. Adams wrote:

 Humm, I did check petsc-3.3 before it was released and petsc-dev works.
 
 Could you attache the debugger and get a line number?
 
 Mark
 
 On Jul 24, 2012, at 10:59 PM, Blaise Bourdin wrote:
 
 Hi,
 
 is pc_gamg_type geo broken in petsc-3.3? I get the following when I make 
 runex54 in  src/ksp/ksp/examples/tutorials:
 
 iMac:tutorials blaise$ make runex54
 1,28c1,128
0 KSP Residual norm 132.598 
1 KSP Residual norm 39.159 
2 KSP Residual norm 15.7856 
3 KSP Residual norm 8.91321 
4 KSP Residual norm 6.95961 
5 KSP Residual norm 15.1387 
6 KSP Residual norm 12.5547 
7 KSP Residual norm 3.78056 
8 KSP Residual norm 2.38412 
9 KSP Residual norm 3.91952 
   10 KSP Residual norm 1.3192 
   11 KSP Residual norm 1.70681 
   12 KSP Residual norm 1.74052 
   13 KSP Residual norm 0.482779 
   14 KSP Residual norm 0.634571 
   15 KSP Residual norm 0.264686 
   16 KSP Residual norm 0.240607 
   17 KSP Residual norm 0.10998 
   18 KSP Residual norm 0.0853072 
   19 KSP Residual norm 0.0551939 
   20 KSP Residual norm 0.0231032 
   21 KSP Residual norm 0.0286234 
   22 KSP Residual norm 0.0110345 
   23 KSP Residual norm 0.00582651 
   24 KSP Residual norm 0.00256534 
   25 KSP Residual norm 0.00176112 
   26 KSP Residual norm 0.000730267 
  Linear solve converged due to CONVERGED_RTOL iterations 26
 ---
 [2]PETSC ERROR: 
 
 [2]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
 probably memory access out of range
 [2]PETSC ERROR: [0]PETSC ERROR: 
 
 [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
 probably memory access out of range
 [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
 [0]PETSC ERROR: or see 
 http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC ERROR: 
 or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory 
 corruption errors
 [0]PETSC ERROR: [1]PETSC ERROR: 
 
 [1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
 probably memory access out of range
 [1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
 [1]PETSC ERROR: or see 
 http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[1]PETSC ERROR: 
 or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory 
 corruption errors
 [1]PETSC ERROR: likely location of problem given in stack below
 [1]PETSC ERROR: -  Stack Frames 
 
 [1]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,
 [1]PETSC ERROR:   INSTEAD the line number of the start of the function
 [1]PETSC ERROR:   is given.
 [1]PETSC ERROR: [1] triangulateAndFormProl line 180 
 /opt/HPC/petsc-3.3/src/ksp/pc/impls/gamg/geo.c
 [1]PETSC ERROR: [1] PCGAMGProlongator_GEO line 756

[petsc-users] pc_gamg_type geo broken in 3.3?

2012-07-24 Thread Blaise Bourdin

Hi,

is pc_gamg_type geo broken in petsc-3.3? I get the following when I make 
runex54 in  src/ksp/ksp/examples/tutorials:

iMac:tutorials blaise$ make runex54
1,28c1,128
   0 KSP Residual norm 132.598 
   1 KSP Residual norm 39.159 
   2 KSP Residual norm 15.7856 
   3 KSP Residual norm 8.91321 
   4 KSP Residual norm 6.95961 
   5 KSP Residual norm 15.1387 
   6 KSP Residual norm 12.5547 
   7 KSP Residual norm 3.78056 
   8 KSP Residual norm 2.38412 
   9 KSP Residual norm 3.91952 
  10 KSP Residual norm 1.3192 
  11 KSP Residual norm 1.70681 
  12 KSP Residual norm 1.74052 
  13 KSP Residual norm 0.482779 
  14 KSP Residual norm 0.634571 
  15 KSP Residual norm 0.264686 
  16 KSP Residual norm 0.240607 
  17 KSP Residual norm 0.10998 
  18 KSP Residual norm 0.0853072 
  19 KSP Residual norm 0.0551939 
  20 KSP Residual norm 0.0231032 
  21 KSP Residual norm 0.0286234 
  22 KSP Residual norm 0.0110345 
  23 KSP Residual norm 0.00582651 
  24 KSP Residual norm 0.00256534 
  25 KSP Residual norm 0.00176112 
  26 KSP Residual norm 0.000730267 
 Linear solve converged due to CONVERGED_RTOL iterations 26
---
 [2]PETSC ERROR: 
 
 [2]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
 probably memory access out of range
 [2]PETSC ERROR: [0]PETSC ERROR: 
 
 [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
 probably memory access out of range
 [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
 [0]PETSC ERROR: or see 
 http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC ERROR: 
 or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory 
 corruption errors
 [0]PETSC ERROR: [1]PETSC ERROR: 
 
 [1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
 probably memory access out of range
 [1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
 [1]PETSC ERROR: or see 
 http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[1]PETSC ERROR: 
 or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory 
 corruption errors
 [1]PETSC ERROR: likely location of problem given in stack below
 [1]PETSC ERROR: -  Stack Frames 
 
 [1]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,
 [1]PETSC ERROR:   INSTEAD the line number of the start of the function
 [1]PETSC ERROR:   is given.
 [1]PETSC ERROR: [1] triangulateAndFormProl line 180 
 /opt/HPC/petsc-3.3/src/ksp/pc/impls/gamg/geo.c
 [1]PETSC ERROR: [1] PCGAMGProlongator_GEO line 756 
 /opt/HPC/petsc-3.3/src/ksp/pc/impls/gamg/geo.c
 [1]PETSC ERROR: [1] PCSetUp_GTry option -start_in_debugger or 
 -on_error_attach_debugger
 [2]PETSC ERROR: or see 
 http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[2]PETSC ERROR: 
 or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory 
 corruption errors
 [2]PETSC ERROR: likely location of problem given in stack below
 [2]PETSC ERROR: -  Stack Frames 
 
 [2]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,
 [2]PETSC ERROR:   INSTEAD the line number of the start of the function
 [2]PETSC ERROR:   is given.
 [2]PETSC ERROR: [2] triangulateAndFormProl line 180 
 /opt/HPC/petsc-3.3/src/ksp/pc/impls/gamg/geo.c
 [2]PETSC ERROR: [2] PCGAMGProlongator_GEO line 756 
 /opt/HPC/petsc-3.3/src/ksp/pc/impls/gamg/geo.c
 [2]PETSC ERROR: [2] PCSetUp_GAMG line 560 
 /opt/HPC/petsc-3.3/src/ksp/pc/impls/gamg/gamg.c
 [2]PETSC ERROR: [2] PCSetUp line 810 
 /opt/HPC/petsc-3.3/src/ksp/pc/interface/precon.c
 [2]PETSC ERROR: [2] KSPSetUp line 182 /opt/HPC/petsc-3.3/src/ksp[3]PETSC 
 ERROR: 
 
 [3]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
 probably memory access out of range
 [3]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
 [3]PETSC ERROR: or see 
 http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[3]PETSC ERROR: 
 or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory 
 corruption errors
 [3]PETSC ERROR: likely location of problem given in stack below
 [3]PETSC ERROR: -  Stack Frames 
 
 [3]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,
 [3]PETSC ERROR:   INSTEAD the line number of the start of the function
 [3]PETSC ERROR:   is given.
 [3]PETSC ERROR: [3] triangulateAndFormProl line 180 
 /opt/HPC/petsc-3.3/src/ksp/pc/impls/gamg/geo.c
 [3]PETSC ERROR: [3] PCGAMGProlongator_GEO line 756 
 /opt/HPC/petsc-3.3/src/ksp/pc/impls/gamg/geo.c
 [3]PETSC ERROR: [3] PCSetUp_Glikely

[petsc-users] block size

2012-07-18 Thread Blaise Bourdin

search for an older thread on this list

possible fixes are to either delete / not create the .info files or to assign
different prefix to each vec

B
On Jul 18, 2012, at 12:39 PM, Anton Popov wrote:

Comment out one of vectors, your code works fine.
I wouldn't say fine because I need to store both vectors.

So what is the workaround?

- I'm not allowed to mix vectors with different block size in the same file?
- I'm not allowed to write more then ONE vector in a file?

On 7/18/12 4:10 PM, Hong Zhang wrote:
You use same datafile 'vec.dat' for writing two different vectors,
V3 and V1:

ierr = VecView(V1, view_out);
ierr = VecView(V3, view_out);

// here, vec.dat holds V3

Then read it in the order
ierr = VecLoad(V1, view_in);
//crash here because reading V3 into V1

ierr = VecLoad(V3, view_in);

Comment out one of vectors, your code works fine.

Hong

On Wed, Jul 18, 2012 at 7:54 AM, Anton Popov popov at uni-mainz.de wrote:
Dear petsc team,

could you please tell me what's wrong with the attached example file?
I run it on 4 processors with petsc-3.3-p1.

What could error message Local size 1000 not compatible with block size 3!
mean?

I've another question related to this issue. What is the real purpose of
PetscViewerBinarySkipInfo function?
I see no reason to skip creating info file, because the file produced by
the attached example seems to be correct.

Moreover, similar block size error occurs in our code while reading file
with multiple vectors, irrespective whether info file exists or not.

Thank you,

Anton

--
Department of Mathematics and Center for Computation Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

-- next part --
An HTML attachment was scrubbed...
URL:
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120718/070d044c/attachment-0001.html

[petsc-users] [petsc-dev] second patch for DMDAVecGetArrayF90

2012-07-11 Thread Blaise Bourdin

On Jul 11, 2012, at 3:44 PM, TAY wee-beng wrote:

 On 11/7/2012 8:10 PM, Satish Balay wrote:
 On Wed, 11 Jul 2012, TAY wee-beng wrote:
 
 I just tried on a different system and the same error occurs. The 1st 
 command
 :
 /
 /
 
 /patch -Np1 -R  DMDAVecRestoreArrayF90-2.patch
 
 $ patch -Np1 -R  DMDAVecRestoreArrayF90-2.patch
 (Stripping trailing CRs from patch.)
 patching file src/dm/impls/da/f90-custom/zda1f90.c
 (Stripping trailing CRs from patch.)
 patching file src/sys/f90-src/f90_cwrap.c
 (Stripping trailing CRs from patch.)
 patching file src/sys/f90-src/fsrc/f90_fwrap.F/
 
 works but the 2nd one doesn't:
 /
 $ patch -Np1   DMDAVecRestoreArrayF90-2.1.patch
 (Stripping trailing CRs from patch.)
 patching file src/dm/impls/da/f90-custom/zda1f90.c
 Hunk #1 FAILED at 175.
 1 out of 1 hunk FAILED -- saving rejects to file
 src/dm/impls/da/f90-custom/zda1f90.c.rej
 (Stripping trailing CRs from patch.)
 patching file src/sys/f90-src/f90_cwrap.c
 Hunk #1 FAILED at 307.
 Hunk #2 FAILED at 333.
 Hunk #3 FAILED at 436.
 3 out of 3 hunks FAILED -- saving rejects to file
 src/sys/f90-src/f90_cwrap.c.rej
 (Stripping trailing CRs from patch.)
 patching file src/sys/f90-src/fsrc/f90_fwrap.F
 Hunk #1 FAILED at 322.
 Hunk #2 FAILED at 426.
 2 out of 2 hunks FAILED -- saving rejects to file
 src/sys/f90-src/fsrc/f90_fwrap.F.rej/
 
 
 It is possible to send me the 3 patch files? Or is there some other ways?
 This doesn't make any sense. What petsc-dev do you have? Is it hg
 repo? or something else?  What revisions? Are the sources locally
 modified?
 
 You can apply DMDAVecRestoreArrayF90-2.patch and then edit
 src/sys/f90-src/f90_cwrap.c and replace all occurances of 'F90ARRAY4d'
 with 'F90ARRAY4D' in the editor to get the same effect.
 
 Satish
 
 I just tried but the same error occurs:
 
 1D:\Lib\petsc-3.3-dev_win32_vs2008/include/finclude/ftn-custom/petscdmcomposite.h90(8):
  warning #6717: This name has not been given an explicit type.   [D1]
 1Linking...
 1libpetsc.lib(f90_cwrap.o) : error LNK2019: unresolved external symbol 
 F90ARRAY4dCREATESCALAR referenced in function F90Array4dCreate
 1libpetsc.lib(f90_cwrap.o) : error LNK2019: unresolved external symbol 
 F90ARRAY4dACCESSFORTRANADDR referenced in function F90Array4dAccess
 1libpetsc.lib(f90_cwrap.o) : error LNK2019: unresolved external symbol 
 F90ARRAY4dACCESSINT referenced in function F90Array4dAccess
 1libpetsc.lib(f90_cwrap.o) : error LNK2019: unresolved external symbol 
 F90ARRAY4dACCESSREAL referenced in function F90Array4dAccess
 1libpetsc.lib(f90_cwrap.o) : error LNK2019: unresolved external symbol 
 F90ARRAY4dACCESSSCALAR referenced in function F90Array4dAccess
 1libpetsc.lib(f90_cwrap.o) : error LNK2019: unresolved external symbol 
 F90ARRAY4dDESTROYSCALAR referenced in function F90Array4dDestroy
 1c:\obj_tmp\dm_test2d\Debug/dm_test2d.exe : fatal error LNK1120: 6 
 unresolved externals
 1
 1Build log written to  file://c:\obj_tmp\dm_test2d\Debug\BuildLog.htm
 
 Any other solutions?


There is something fishy in your tree or your build: the patch that Satish sent 
addressed this issue, but it clearly has not been applied to your tree. The 
patch I had originally sent contained an error. For instance, 
F90ARRAY4DACCESSFORTRANADDR was replaced by F90ARRAY4dACCESSFORTRANADDR (note 
the capitalization of the d). These are the symbols that are unresolved in 
your build.



How about starting from a clean state: 
- rename your current petsc-dev tree
- clone petsc-dev from the official tree (BTW, do you really need petsc-dev, 
and not 3.3?)
- apply the patch Satish sent you or, quoting Satich 
 edit src/sys/f90-src/f90_cwrap.c and replace all occurrences of 'F90ARRAY4d' 
 with 'F90ARRAY4D'

- run configure with your favorite options, then make

This should fix your problem.

Blaise

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin







-- next part --
An HTML attachment was scrubbed...
URL: 
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120711/cde280b1/attachment.html

[petsc-users] [petsc-dev] second patch for DMDAVecGetArrayF90

2012-07-11 Thread Blaise Bourdin

Rectification: clone petsc-3.3, not petsc-dev, of course...
B

 There is something fishy in your tree or your build: the patch that Satish 
 sent addressed this issue, but it clearly has not been applied to your tree. 
 The patch I had originally sent contained an error. For instance, 
 F90ARRAY4DACCESSFORTRANADDR was replaced by F90ARRAY4dACCESSFORTRANADDR (note 
 the capitalization of the d). These are the symbols that are unresolved in 
 your build.
 
 
 
 How about starting from a clean state: 
 - rename your current petsc-dev tree
 - clone petsc-dev from the official tree (BTW, do you really need petsc-dev, 
 and not 3.3?)
 - apply the patch Satish sent you or, quoting Satich 
 edit src/sys/f90-src/f90_cwrap.c and replace all occurrences of 
 'F90ARRAY4d' with 'F90ARRAY4D'
 
 - run configure with your favorite options, then make
 
 This should fix your problem.

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin







-- next part --
An HTML attachment was scrubbed...
URL: 
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120711/aa046b61/attachment.html

[petsc-users] [petsc-dev] Getting 2d array with updated ghost values from DM global vector

2012-07-06 Thread Blaise Bourdin


On Jul 6, 2012, at 6:02 PM, TAY wee-beng wrote:

 On 3/7/2012 1:23 PM, Satish Balay wrote:
 On Tue, 3 Jul 2012, Barry Smith wrote:
 
 On Jul 3, 2012, at 3:08 AM, Blaise Bourdin wrote:
 
 On Jul 3, 2012, at 4:10 AM, Barry Smith wrote:
 
 Blaise,
 
  I don't understand why the patch does anything:
 
 -  *ierr = VecRestoreArray(*v,0);if (*ierr) return;
 +  PetscScalar *fa;
 +  *ierr = F90Array1dAccess(a,PETSC_SCALAR,(void**)fa 
 PETSC_F90_2PTR_PARAM(ptrd));
 +  *ierr = VecRestoreArray(*v,fa);if (*ierr) return;
 *ierr = F90Array1dDestroy(a,PETSC_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
 
 All that passing fa into VecRestoreArray() does is cause fa to be 
 zeroed. Why would that have any affect on anything?
 
 Not sure either, I quite don't understand this code, but I noticed that 
 the logic of VecRestoreArrayF90 was different from that of 
 DMDAVecRestoreArrayF90
 
 src/vec/vec/interface/f90-custom/zvectorf90.c:33
 PetscScalar *fa;
 *__ierr = F90Array1dAccess(ptr,PETSC_SCALAR,(void**)fa 
 PETSC_F90_2PTR_PARAM(ptrd));if (*__ierr) return;
 *__ierr = F90Array1dDestroy(ptr,PETSC_SCALAR 
 PETSC_F90_2PTR_PARAM(ptrd));if (*__ierr) return;
  It could be the above line is important; but the Accesser and restore 
 array are not.
 
  I'll have Satish apply the patch.
 pushed to petsc-3.3 [petsc-dev will get this update]
 
 Satish
 Hi,
 
 I just tested with the latest petsc-dev but it doesn't work in intel linux 
 for ex11f90. Has the patch been applied?
Try to clone petsc-3.3  from the mercurial repository 
http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/ It looks like the first patch 
has not made its way to the tarball or petsc-dev yet. 
You can also apply the patch manually: 
cd $PETSC_DIR
patch -p1  DMDAVecGetArrayF90.patch

 Also, is there any chance of it working under 3d with multiple dof since 
 that's what I'm using and I have other problems with gfortran. Lastly, if the 
 patch is applied, it works with 3d da with 1 dof? Is that right?

I am sending another patch that should take care of the 3d case with 1 dof to 
the developers list.

Blaise

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] [petsc-dev] Getting 2d array with updated ghost values from DM global vector

2012-07-03 Thread Blaise Bourdin


On Jul 3, 2012, at 4:10 AM, Barry Smith wrote:

 
   Blaise,
 
I don't understand why the patch does anything:
 
 -  *ierr = VecRestoreArray(*v,0);if (*ierr) return;
 +  PetscScalar *fa;
 +  *ierr = F90Array1dAccess(a,PETSC_SCALAR,(void**)fa 
 PETSC_F90_2PTR_PARAM(ptrd));
 +  *ierr = VecRestoreArray(*v,fa);if (*ierr) return;
   *ierr = F90Array1dDestroy(a,PETSC_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
 
 All that passing fa into VecRestoreArray() does is cause fa to be zeroed. 
 Why would that have any affect on anything?


Not sure either, I quite don't understand this code, but I noticed that the 
logic of VecRestoreArrayF90 was different from that of DMDAVecRestoreArrayF90

src/vec/vec/interface/f90-custom/zvectorf90.c:33
  PetscScalar *fa;
  *__ierr = F90Array1dAccess(ptr,PETSC_SCALAR,(void**)fa 
PETSC_F90_2PTR_PARAM(ptrd));if (*__ierr) return;
  *__ierr = F90Array1dDestroy(ptr,PETSC_SCALAR PETSC_F90_2PTR_PARAM(ptrd));if 
(*__ierr) return;
  *__ierr = VecRestoreArray(*x,fa);

Why aren't the calls to F90Array1dAccess and F90Array1dDestroy necessary in the 
context of DMDAVecGetArrayF90?

Blaise



 
   Thanks
 
 Barry
 
 On Jul 2, 2012, at 10:10 AM, Blaise Bourdin wrote:
 
 Hi,
 
 There appears to be a bug in DMDAVecRestoreArrayF90. It is probably only 
 triggered when the intel compilers. gfortran and intel seem to have very 
 different internal implementations of fortran90 allocatable arrays. 
 
 Developers, can you check if the attached patch makes sense? It will not fix 
 the case of a 3d da with dof1 since F90Array4dAccess is not implemented. 
 Other than that, it seems to fix ex11f90 under linux and mac OS
 
 DMDAVecRestoreArrayF90.patch
 
 Blaise
 
 
 
 On Jul 2, 2012, at 5:41 PM, TAY wee-beng wrote:
 
 On 2/7/2012 2:49 PM, Matthew Knepley wrote:
 On Mon, Jul 2, 2012 at 3:58 AM, TAY wee-beng zonexo at gmail.com wrote:
 Hi,
 
 I have used DMDACreate2d for my code and then use:
 
 call 
 DMLocalToGlobalBegin(da,b_rhs_semi_local,INSERT_VALUES,b_rhs_semi_global,ierr)
 
 call 
 DMLocalToGlobalEnd(da,b_rhs_semi_local,INSERT_VALUES,b_rhs_semi_global,ierr)
 
 to construct the global DM vector b_rhs_semi_global
 
 Now I want to get the values with ghost values in a 2d array locally which 
 is declared as:
 
 real(8), allocatable :: array2d(:,:)
 
 I guess I should use DMDAGetGhostCorners to get the corressponding indices 
 and allocate it.  But what should I do next? How can I use something like 
 VecGetArrayF90 to get to the pointer to access the local vector?
 
 I can't use DMDAVecGetArrayF90/DMDAVecRestoreArrayF90 since I'm using 
 intel fortran and they can't work. I can't use gfortran at the moment 
 since I've problems with HYPRE with   gfortran in 3D.
 
 Are you certain of this? That used to be true, but the current version 
 should work for any F90.
 
   Matt
 
 I just tested 3.3-p1 and it still doesn't work (example ex11f90 in dm). Is 
 there a chance petsc-dev can work?
 
 Thanks
 
 -- 
 Yours sincerely,
 
 TAY wee-beng
 
 
 
 
 -- 
 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
 
 
 
 -- 
 Department of Mathematics and Center for Computation  Technology
 Louisiana State University, Baton Rouge, LA 70803, USA
 Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 
 http://www.math.lsu.edu/~bourdin
 
 
 
 
 
 
 
 

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] [petsc-dev] Getting 2d array with updated ghost values from DM global vector

2012-07-03 Thread Blaise Bourdin


On Jul 3, 2012, at 4:10 AM, Barry Smith wrote:

 
  Blaise,
 
   I don't understand why the patch does anything:
 
 -  *ierr = VecRestoreArray(*v,0);if (*ierr) return;
 +  PetscScalar *fa;
 +  *ierr = F90Array1dAccess(a,PETSC_SCALAR,(void**)fa 
 PETSC_F90_2PTR_PARAM(ptrd));
 +  *ierr = VecRestoreArray(*v,fa);if (*ierr) return;
  *ierr = F90Array1dDestroy(a,PETSC_SCALAR PETSC_F90_2PTR_PARAM(ptrd));
 
 All that passing fa into VecRestoreArray() does is cause fa to be zeroed. 
 Why would that have any affect on anything?


Not sure either, I quite don't understand this code, but I noticed that the 
logic of VecRestoreArrayF90 was different from that of DMDAVecRestoreArrayF90

src/vec/vec/interface/f90-custom/zvectorf90.c:33
 PetscScalar *fa;
 *__ierr = F90Array1dAccess(ptr,PETSC_SCALAR,(void**)fa 
PETSC_F90_2PTR_PARAM(ptrd));if (*__ierr) return;
 *__ierr = F90Array1dDestroy(ptr,PETSC_SCALAR PETSC_F90_2PTR_PARAM(ptrd));if 
(*__ierr) return;
 *__ierr = VecRestoreArray(*x,fa);

Why aren't the calls to F90Array1dAccess and F90Array1dDestroy necessary in the 
context of DMDAVecGetArrayF90?

Blaise



 
  Thanks
 
Barry
 
 On Jul 2, 2012, at 10:10 AM, Blaise Bourdin wrote:
 
 Hi,
 
 There appears to be a bug in DMDAVecRestoreArrayF90. It is probably only 
 triggered when the intel compilers. gfortran and intel seem to have very 
 different internal implementations of fortran90 allocatable arrays. 
 
 Developers, can you check if the attached patch makes sense? It will not fix 
 the case of a 3d da with dof1 since F90Array4dAccess is not implemented. 
 Other than that, it seems to fix ex11f90 under linux and mac OS
 
 DMDAVecRestoreArrayF90.patch
 
 Blaise
 
 
 
 On Jul 2, 2012, at 5:41 PM, TAY wee-beng wrote:
 
 On 2/7/2012 2:49 PM, Matthew Knepley wrote:
 On Mon, Jul 2, 2012 at 3:58 AM, TAY wee-beng zonexo at gmail.com wrote:
 Hi,
 
 I have used DMDACreate2d for my code and then use:
 
 call 
 DMLocalToGlobalBegin(da,b_rhs_semi_local,INSERT_VALUES,b_rhs_semi_global,ierr)
 
 call 
 DMLocalToGlobalEnd(da,b_rhs_semi_local,INSERT_VALUES,b_rhs_semi_global,ierr)
 
 to construct the global DM vector b_rhs_semi_global
 
 Now I want to get the values with ghost values in a 2d array locally which 
 is declared as:
 
 real(8), allocatable :: array2d(:,:)
 
 I guess I should use DMDAGetGhostCorners to get the corressponding indices 
 and allocate it.  But what should I do next? How can I use something like 
 VecGetArrayF90 to get to the pointer to access the local vector?
 
 I can't use DMDAVecGetArrayF90/DMDAVecRestoreArrayF90 since I'm using 
 intel fortran and they can't work. I can't use gfortran at the moment 
 since I've problems with HYPRE with   gfortran in 3D.
 
 Are you certain of this? That used to be true, but the current version 
 should work for any F90.
 
  Matt
 
 I just tested 3.3-p1 and it still doesn't work (example ex11f90 in dm). Is 
 there a chance petsc-dev can work?
 
 Thanks
 
 -- 
 Yours sincerely,
 
 TAY wee-beng
 
 
 
 
 -- 
 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
 
 
 
 -- 
 Department of Mathematics and Center for Computation  Technology
 Louisiana State University, Baton Rouge, LA 70803, USA
 Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 
 http://www.math.lsu.edu/~bourdin
 
 
 
 
 
 
 
 

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] Getting 2d array with updated ghost values from DM global vector

2012-07-02 Thread Blaise Bourdin

Hi,

There appears to be a bug in DMDAVecRestoreArrayF90. It is probably only
triggered when the intel compilers. gfortran and intel seem to have very
different internal implementations of fortran90 allocatable arrays.

Developers, can you check if the attached patch makes sense? It will not fix
the case of a 3d da with dof1 since F90Array4dAccess is not implemented. Other
than that, it seems to fix ex11f90 under linux and mac OS

Blaise

On Jul 2, 2012, at 5:41 PM, TAY wee-beng wrote:

On 2/7/2012 2:49 PM, Matthew Knepley wrote:
On Mon, Jul 2, 2012 at 3:58 AM, TAY wee-beng zonexo at gmail.com wrote:
Hi,

I have used DMDACreate2d for my code and then use:

call
DMLocalToGlobalBegin(da,b_rhs_semi_local,INSERT_VALUES,b_rhs_semi_global,ierr)

call
DMLocalToGlobalEnd(da,b_rhs_semi_local,INSERT_VALUES,b_rhs_semi_global,ierr)

to construct the global DM vector b_rhs_semi_global

Now I want to get the values with ghost values in a 2d array locally which
is declared as:

real(8), allocatable :: array2d(:,:)

I guess I should use DMDAGetGhostCorners to get the corressponding indices
and allocate it. But what should I do next? How can I use something like
VecGetArrayF90 to get to the pointer to access the local vector?

I can't use DMDAVecGetArrayF90/DMDAVecRestoreArrayF90 since I'm using intel
fortran and they can't work. I can't use gfortran at the moment since I've
problems with HYPRE with gfortran in 3D.

Are you certain of this? That used to be true, but the current version
should work for any F90.

Matt

I just tested 3.3-p1 and it still doesn't work (example ex11f90 in dm). Is
there a chance petsc-dev can work?

Thanks

--
Yours sincerely,

TAY wee-beng

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

-- next part --
An HTML attachment was scrubbed...
URL:
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120702/c28a29ed/attachment.html
-- next part --
A non-text attachment was scrubbed...
Name: DMDAVecRestoreArrayF90.patch
Type: application/octet-stream
Size: 2211 bytes
Desc: not available
URL:
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120702/c28a29ed/attachment.obj
-- next part --
An HTML attachment was scrubbed...
URL:
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120702/c28a29ed/attachment-0001.html

[petsc-users] VecValid

2012-05-03 Thread Blaise Bourdin

Hi,

As a side note: this is actually a problem I have met in one of the TS 
examples. When using fortran datatypes, the statements A=0 and if (A .eq. 
0) are not valid. They would require overloading the assignment and comparison 
operators.
Can anybody see an alternative?

Blaise

 
  VecValid was unsafe because it looks inside a data structure that may or not 
 be valid. 
 
   In Fortran the safe way to do this is to set the value to 0 right after 
 declaration and then check if not zero. So
 
   MatA 
 
   A  = 0
 
...
 
if (A .eq. 0) then 
 call VecCreate(   A,ierr)
endif
 
   Barry
 
 On May 3, 2012, at 9:17 PM, Sanjay Govindjee wrote:
 
 The entry points for our code can vary.  Thus a vector for example can get 
 created in a handful
 of routines.  But it is possible that the vector has already been created in 
 another routine.  So
 we perform a VecValid check on it.  If it is false we create it, else we go 
 on and use it.  Also at
 the end when we exit the program we clean up with a series of VecDestroy 
 calls but we check first if the vector is valid before calling VecDestroy on 
 it.  Same for MatValid/MatDestroy.
 
 Obviously we can store a logical flag in a common block or the like upon 
 creation; but VecValid and
 MatValid were handy utility functions for us.
 
 Is is feasible to call PetscValidHeaderSpecific( ) out of fortran?  Or 
 should should I just set up my own
 collection of allocate/deallocate flags.
 
 -sanjay
 
 On 5/3/12 7:08 PM, Matthew Knepley wrote:
 On Thu, May 3, 2012 at 9:20 PM, Sanjay Govindjee s_g at berkeley.edu 
 wrote:
 We have some code that has worked up through petsc-3.1 and we are in the 
 process of updating it to petsc-3.2.
 
 One thing that I can not find mentioned in the change logs (looking all the 
 way back to petsc-2.1.0) or the FAQs is the elimination of the VecValid and 
 MatValid tests.
 
 In C, this is now the macro PetscValidHeaderSpecific(). This check now 
 happens in every function call. What
 logic would need to look for a corrupt Vec beyond CHKERRQ?
 
  Thanks,
 
  Matt
 
 We used to perform operations like:
 
 call VecValid(xvec, chk, ierr)
 
 if(chk .eqv. PETSC_FALSE) then
   call VecCreate(PETSC_COMM_WORLD, xvec, ierr)
   call VecSetSizes  (xvec, numpeq, PETSC_DECIDE, ierr)
   call VecSetFromOptions(xvec, ierr)
 else
   
 endif
 
 But it seems VecValid and MatValid have been eliminated.  Is there a reason 
 for this?  or have
 I made a mistake?  Is there a replacement test for invalid Vec and Mat 
 pointers?
 
 -sanjay
 
 
 
 
 -- 
 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
 
 -- 
 ---
 Sanjay Govindjee, PhD, PE
 Professor of Civil Engineering
 
 779 Davis Hall
 Structural Engineering, Mechanics and Materials
 Department of Civil Engineering
 University of California
 Berkeley, CA 94720-1710
 
 Voice:  +1 510 642 6060
 FAX:+1 510 643 5264
 
 s_g at berkeley.edu
 http://www.ce.berkeley.edu/~sanjay
 
 ---
 
 

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] ex52_integrateElement.cu

2012-03-27 Thread Blaise Bourdin

On Mar 27, 2012, at 1:23 PM, Matthew Knepley wrote:

On Tue, Mar 27, 2012 at 12:58 PM, David Fuentes fuentesdt at gmail.com
wrote:
Hi,

I had a question about the status of example 52.

http://petsc.cs.iit.edu/petsc/petsc-dev/file/a8e2f2c19319/src/snes/examples/tutorials/ex52.c
http://petsc.cs.iit.edu/petsc/petsc-dev/file/a8e2f2c19319/src/snes/examples/tutorials/ex52_integrateElement.cu

Can this example be used with a DM object created from an unstructured
exodusII mesh, DMMeshCreateExodus, And the FEM assembly done on GPU ?

1) I have pushed many more tests for it now. They can be run using the Python
build system

./config/builder2.py check src/snes/examples/tutorials/ex52.c

in fact, you can build any set of files this way.

2) The Exodus creation has to be converted to DMComplex from DMMesh. That
should not take me very long. Blaise maintains that
so maybe there will be help :) You will just replace
DMComplexCreateBoxMesh() with DMComplexCreateExodus(). If you request
it, I will bump it up the list.

DMMeshCreateExodusNG is much more flexible than DMMeshCreateExodus in that it
can read meshes with multiple element types and should have a much lower memory
footprint. The code should be fairly easy to read. you can email me directly if
you have specific questions. I had looked at creating a DMComplex and it did
not look too difficult, as long as interpolation is not needed. I have plans to
write DMComplexCreateExodus, but haven't had time too so far. Updating the Vec
viewers and readers may be a bit more involved. In perfect world, one would
write an EXODUS viewer following the lines of the VTK and HDF5 ones.

Blaise

Let me know if you can run the tests.

Thanks

Matt

Thanks,
David
--
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

-- next part --
An HTML attachment was scrubbed...
URL:
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120327/0252007f/attachment.htm

[petsc-users] Binary VTK viewer

2012-03-08 Thread Blaise Bourdin

No. 

Imtiaz had to leave LSU before we could accomplish anything meaningful.

Blaise

On Mar 8, 2012, at 11:46 AM, Max Rudolph wrote:

 I was looking into the various PetscViewerFormats available and it appears 
 that only ASCII vtk format is supported. I found an old thread on this 
 mailing list with talk between Blaise Bourdin and Matt Knepley of adding 
 support for output of binary vtk files. Was this ever done?
 
 Thanks for your help,
 Max

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] ISAllGather withoutduplicates

2012-02-17 Thread Blaise Bourdin

Hi,

Is there an easy way to gather all values of an IS across all processes in a 
communicator while removing duplicates?

Basically, I want to go from 
[0] Number of indices in set 2
[0] 0 1
[0] 1 2
[1] Number of indices in set 2
[1] 0 2
[1] 1 3

to
[0] Number of indices in set 3
[0] 0 1
[0] 1 2
[0] 2 3
[1] Number of indices in set 3
[1] 0 1
[1] 1 2
[1] 2 3

The way I do it right now is 
  ierr = ISGetTotalIndices(csIS,labels);CHKERRQ(ierr);
  ierr = ISGetSize(csIS,num_cs_global);CHKERRQ(ierr);
  ierr = PetscMalloc(num_cs_global * sizeof(PetscInt),labels2);
  for (i = 0; i  num_cs_global; i++) {
labels2[i] = labels[i];
  }
  ierr = PetscSortRemoveDupsInt(num_cs_global,labels2);CHKERRQ(ierr);
  ierr = 
ISCreateGeneral(comm,num_cs_global,labels2,PETSC_COPY_VALUES,csIS_global);CHKERRQ(ierr);
  ierr = ISRestoreTotalIndices(csIS,labels);CHKERRQ(ierr);
  ierr = PetscFree(labels2);CHKERRQ(ierr);

but there has to be a better way (or at least one that does not involve copying 
from const PetscInt *labels to PetscInt *labels2, and then uses again 
PETSC_COPY_VALUES).


Blaise
-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] ISAllGather withoutduplicates

2012-02-17 Thread Blaise Bourdin

 The way I do it right now is
  ierr = ISGetTotalIndices(csIS,labels);CHKERRQ(ierr);
  ierr = ISGetSize(csIS,num_cs_global);CHKERRQ(ierr);
 
 I would violate PETSc semantics here since you are going to destroy csIS 
 anyway:
 
   PetscSortRemoveDupsInt(num_cs_global, labels);

That's what I was thinking too, but PetscSortRemoveDupsInt expect a PetscInt* 
not a const PetscInt*

It's not a big deal, I can live with the 2 copies, considering that the local 
size of the IS is going to be quite small.

Blaise




   ISCreateGeneral(comm, num_cs_global, labels, PETSC_COPY_VALUES, 
 csIS_global);
  
  ierr = ISRestoreTotalIndices(csIS,labels);CHKERRQ(ierr);
 
Matt
  
 
 Blaise
 --
 Department of Mathematics and Center for Computation  Technology
 Louisiana State University, Baton Rouge, LA 70803, USA
 Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 
 http://www.math.lsu.edu/~bourdin
 
 
 
 
 
 
 
 
 
 
 -- 
 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

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin







-- next part --
An HTML attachment was scrubbed...
URL: 
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120217/d402f129/attachment.htm

[petsc-users] ISAllGather withoutduplicates

2012-02-17 Thread Blaise Bourdin


On Feb 17, 2012, at 2:45 PM, Matthew Knepley wrote:

 On Fri, Feb 17, 2012 at 2:39 PM, Blaise Bourdin bourdin at lsu.edu wrote:
 The way I do it right now is
  ierr = ISGetTotalIndices(csIS,labels);CHKERRQ(ierr);
  ierr = ISGetSize(csIS,num_cs_global);CHKERRQ(ierr);
 
 I would violate PETSc semantics here since you are going to destroy csIS 
 anyway:
 
   PetscSortRemoveDupsInt(num_cs_global, labels);
 
  PetscSortRemoveDupsInt(num_cs_global, (PetscInt *) labels);

Facepalm...

Blaise

 
   Matt
  
 That's what I was thinking too, but PetscSortRemoveDupsInt expect a PetscInt* 
 not a const PetscInt*
 
 It's not a big deal, I can live with the 2 copies, considering that the local 
 size of the IS is going to be quite small.
 
 Blaise
 
 
 
 
   ISCreateGeneral(comm, num_cs_global, labels, PETSC_COPY_VALUES, 
 csIS_global);
  
  ierr = ISRestoreTotalIndices(csIS,labels);CHKERRQ(ierr);
 
Matt
  
 
 Blaise
 --
 Department of Mathematics and Center for Computation  Technology
 Louisiana State University, Baton Rouge, LA 70803, USA
 Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 
 http://www.math.lsu.edu/~bourdin
 
 
 
 
 
 
 
 
 
 
 -- 
 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
 
 -- 
 Department of Mathematics and Center for Computation  Technology
 Louisiana State University, Baton Rouge, LA 70803, USA
 Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 
 http://www.math.lsu.edu/~bourdin
 
 
 
 
 
 
 
 
 
 
 -- 
 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

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin







-- next part --
An HTML attachment was scrubbed...
URL: 
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120217/c35cfbdd/attachment-0001.htm

[petsc-users] SNESVI convergence spped

2012-01-16 Thread Blaise Bourdin

Hi,

Ata and I are working together on this. The problem he describes is 1/2 of the
iteration of our variational fracture code.
In our application, E is position dependant, and typically becomes very large
along very thin bands with width of the order of epsilon in the domain.
Essentially, we expect that V will remain exactly equal to 1 almost everywhere,
and will transition to 0 on these bands. Of course, we are interested in the
limit as epsilon goes to 0.

If the problem indeed is that it takes many steps to add the degrees of
freedom. Is there any way to initialize manually the list of active
constraints? To give you an idea, here is a link to a picture of the type of
solution we expect. blue=1
https://www.math.lsu.edu/~bourdin/377451-.png

Blaise

It seems to me that the problem is that ultimately ALL of the degrees of
freedom are in the active set,
but they get added to it a few at a time -- and there may even be some
chatter there -- necessitating many SNESVI steps.
Could it be that the regularization makes things worse? When \epsilon \ll 1,
the unconstrained solution is highly oscillatory, possibly further
exacerbating the problem. It's possible that it would be better if V just
diverged uniformly. Then nearly all of the degrees of freedom would bump up
against the upper obstacle all at once.

Dmitry.

On Mon, Jan 16, 2012 at 8:05 PM, Barry Smith bsmith at mcs.anl.gov wrote:

What do you get with -snes_vi_monitor it could be it is taking a while to
get the right active set.

Barry

On Jan 16, 2012, at 6:20 PM, Ataollah Mesgarnejad wrote:

Dear all,

I'm trying to use SNESVI to solve a quadratic problem with box constraints.
My problem in FE context reads:

(\int_{Omega} E phi_i phi_j + \alpha \epsilon dphi_i dphi_j dx) V_i -
(\int_{Omega} \alpha \frac{phi_j}{\epsilon} dx) = 0 , 0= V = 1

or:

[A]{V}-{b}={0}

here phi is the basis function, E and \alpha are positive constants, and
\epsilon is a positive regularization parameter in order of mesh
resolution. In this problem we expect V =1 a.e. and go to zero very fast
at some places.
I'm running this on a rather small problem (50 DOFS) on small number
of processors (72). I expected SNESVI to converge in couple of iterations
(10) since my A matrix doesn't change, however I'm experiencing a slow
convergence (~50-70 iterations). I checked KSP solver for SNES and it
converges with a few iterations.

I would appreciate any suggestions or observations to increase the
convergence speed?

Best,
Ata

-- next part --
An HTML attachment was scrubbed...
URL:
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120116/86bd6064/attachment-0001.htm

[petsc-users] query about PETSc usage in Fortran applications

2011-10-27 Thread Blaise Bourdin

Lois,

I use petsc through the fortran interface for my variational fracture mechanics 
code. (unstructured finite elements, 2d-3d). petsc has been instrumental in 
getting a parallel version. Sieve really helped me getting to large problems 
(largest to date being a 24M elements, 2,400 cores simulation, after which I 
cannot find tools to post-process my results)


The application area is somewhere between engineering, computational mechanics, 
and applied mathematics

I currently use Sieve, the KSP solvers and am in the (slow) process of adding 
TS. In the future, I will try to use VI in order to replace the optimization 
routines from TAO.

I have acknowledged petsc in the following publications (I don't think that any 
of them is listed on the petsc web page).
[Bourdin et al., 2011] Bourdin, B., Larsen, C., and Richardson, C. (2011). A 
time-discrete model for dynamic fracture based on crack regularization. 
International Journal of Fracture, 168:133?143. 10.1007/s10704- 010-9562-x.
[Bourdin et al., 2010] Bourdin, B., Bucur, D., and Oudet, E. (2009/2010). 
Optimal partitions. SIAM J. Sci. Comput., 31(6):4100?4114.
[Bourdin et al., 2008] Bourdin, B., Francfort, G., and Marigo, J.-J. (2008). 
The Variational Approach to Fracture. (reprinted from J. Elasticity 
91(1-3):1?148, 2008). Springer.
[Bourdin et al., 2008] Bourdin, B., Francfort, G., and Marigo, J.-J. (2008). 
The variational approach to fracture. J. Elasticity, 91(1-3):1?148.
[Kimn and Bourdin, 2007] Kimn, J.-H. and Bourdin, B. (2007). Numerical 
implementation of overlapping balancing domain decomposition methods on 
unstructured meshes. In Widlund, O. B. and Keyes, D. E., editors, Domain 
Decomposition Methods in Science and Engineering XVI, volume 55 of Lecture 
Notes in Computational Science and Engineering, pages 309?315. Springer-Verlag.

as well as in several conference proceedings, and talks.

Blaise

On Oct 25, 2011, at 7:46 AM, Lois Curfman McInnes wrote:

 
 I am collecting information about PETSc use in Fortran applications.  If you 
 are a using PETSc via the Fortran interface, please send email (to me only, 
 curfman at mcs.anl.gov) to indicate:
- application area
- what parts of PETSc are used 
- pointer to any publications or other references 
 
 Of particular interest are applications in which PETSc facilitated a 
 transition to parallelism for an existing application that had previously 
 been only sequential.  
 
 Thanks,
 Lois
 

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin







-- next part --
An HTML attachment was scrubbed...
URL: 
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20111027/f9c2d3bd/attachment-0001.htm

[petsc-users] Questions about TS

2011-10-19 Thread Blaise Bourdin

Hi,

I am trying to use TS to solve a simple transient problem in an unstructured 
finite element f90 code.

1. Section 6.1.1 of the manual refers to a TSSetMatrices function that can be 
used to set the RHS and LHS matrices, but I can;t find it. Is this section 
outdated?
2. Since we are using unstructured finite elements, the LHS matrix is not the 
identity. As far as I understand, we have two possible choices:
   - Use a mass lumping approximation of the variational identity matrix (mass 
matrix), M, and use M^{-1}K for the RHS matrix instead of K.
   - Use an IMEX method where the implicit matrix is the variational identity M.
 Is this right? What is the recommended way?
3. Is src/ts/examples/tests/ex1f.F supposed to work? Is there a fortran or 
fortran90 example I could start from? I get the following:
galerkin:tests bourdin$ ./ex1f -linear_constant_matrix
   1  Procs Avg. error 2 normNaN
   NaN  euler  
galerkin:tests bourdin$ ./ex1f
   1  Procs Avg. error 2 normNaN
   NaN  euler  

Thanks,

Blaise 

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] Questions about TS

2011-10-19 Thread Blaise Bourdin

Hi,


 On Wed, Oct 19, 2011 at 10:54, Blaise Bourdin bourdin at math.lsu.edu wrote:
 Hi,
 
 I am trying to use TS to solve a simple transient problem in an unstructured 
 finite element f90 code.
 
 1. Section 6.1.1 of the manual refers to a TSSetMatrices function that can be 
 used to set the RHS and LHS matrices, but I can;t find it. Is this section 
 outdated?
 
 Yes, I must have missed this section when updating the documentation.

OK, that makes more sense now.
  
 2. Since we are using unstructured finite elements, the LHS matrix is not the 
 identity. As far as I understand, we have two possible choices:
   - Use a mass lumping approximation of the variational identity matrix (mass 
 matrix), M, and use M^{-1}K for the RHS matrix instead of K.
 
 You can also write this as a special case of the choice below where you use 
 -ksp_type preonly -pc_type jacobi.
I am not sure I am following you. Can you elaborate?

   - Use an IMEX method where the implicit matrix is the variational identity 
 M.
  Is this right? What is the recommended way?
 
 I would just do this because it's the most flexible. See the user's manual 
 section on IMEX methods. If you are interested in adaptive error control, 
 then you should also check out -ts_type rosw in petsc-dev. In any case, you 
 can write your mass matrix as well as any stiff terms that you want to treat 
 implicitly into TSSetIFunction(), provide an (approximate) Jacobian with 
 TSSetIJacobian(), and put the rest in TSSetRHSFunction(). You can be even 
 more sloppy about it with TSROSW.
 
 Look at src/ts/examples/tutorials/ex22.c (has a Fortran twin, ex22f.F) or 
 ex25.c in petsc-dev.
 

How about recasting the problem as a DAE? The documentation seems to imply that 
this is feasible. For ODE with nontrivial mass matrices such as arise in FEM, 
the implicit/DAE interface significantly reduces overhead to prepare the system 
for algebraic solvers (SNES/KSP) by having the user assemble the correctly 
shifted matrix.

Following ex15, solving 
u_t = \Delta u
can be recast as solving 
F(t,u,\dot u) = 0
with
F(t,u,v) = v-\Delta u.

In this case, the IJacobian would be 
K+aM,
where K is the stiffness matrix (K_{i,j} = \int_\Omega \nabla \phi_i \cdot 
\nabla \phi_j\,dx) and M the mass matrix (M_{i,j} = \int_\Omega \phi_i  
\phi_j\,dx) .

At the continuous level, the IFunction would be v-\Delta u, which cannot be 
evaluated directly in the finite element basis by either solving
M F = M\dot u + Ku
or using mass lumping.

Am I expected to do this in my IFunction or is there a way to pass the mass 
matrix to the TS?

As far as I understand, using IMPEX will lead to the same issue regardless of 
the way I write my problem,  i.e. wether I write G(t,u,v) =  v-\Delta u and 
F(t,u) = 0, or G(t,u,v) =  v and F(t,u) = \Delta u.

As far as I can see, all examples and tests use structured meshes where the 
mass matrix is the identity matrix.

Thanks for the clarification.
Blaise 
-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] Questions about TS

2011-10-19 Thread Blaise Bourdin

Dear Jed,

I get it now... My confusion was due to being used to derived all scheme on the 
pde, then discretizing, whereas the documentation assumes that the equation is 
already discretized. I should have figured it out. 

Is it right to think of the division between ODE, DAE and IMEX in the 
documentation as Fully Explicit vs. Fully implicit vs. Semi-implicit?

Blaise

On Oct 19, 2011, at 12:55 PM, Jed Brown wrote:

 On Wed, Oct 19, 2011 at 12:32, Blaise Bourdin bourdin at lsu.edu wrote:
 Hi,
 
 
  On Wed, Oct 19, 2011 at 10:54, Blaise Bourdin bourdin at math.lsu.edu 
  wrote:
  Hi,
 
  I am trying to use TS to solve a simple transient problem in an 
  unstructured finite element f90 code.
 
  1. Section 6.1.1 of the manual refers to a TSSetMatrices function that can 
  be used to set the RHS and LHS matrices, but I can;t find it. Is this 
  section outdated?
 
  Yes, I must have missed this section when updating the documentation.
 
 OK, that makes more sense now.
 
  2. Since we are using unstructured finite elements, the LHS matrix is not 
  the identity. As far as I understand, we have two possible choices:
- Use a mass lumping approximation of the variational identity matrix 
  (mass matrix), M, and use M^{-1}K for the RHS matrix instead of K.
 
  You can also write this as a special case of the choice below where you use 
  -ksp_type preonly -pc_type jacobi.
 I am not sure I am following you. Can you elaborate?
 
 To solve
 
 M*Xdot = F(X)
 
 IFunction(X,Xdot) = M*Xdot
 IJacobian(X,Xdot) = M
 RHSFunction(X) = F(X)
 
 Now if you lump M, then you can solve it exactly with one iteration of 
 Jacobi. My options were just turn off any extra norms that would normally be 
 computed by an implicit solve. If M is the consistent mass matrix, then you 
 will need some iterations.
  
 
- Use an IMEX method where the implicit matrix is the variational 
  identity M.
   Is this right? What is the recommended way?
 
  I would just do this because it's the most flexible. See the user's manual 
  section on IMEX methods. If you are interested in adaptive error control, 
  then you should also check out -ts_type rosw in petsc-dev. In any case, you 
  can write your mass matrix as well as any stiff terms that you want to 
  treat implicitly into TSSetIFunction(), provide an (approximate) Jacobian 
  with TSSetIJacobian(), and put the rest in TSSetRHSFunction(). You can be 
  even more sloppy about it with TSROSW.
 
  Look at src/ts/examples/tutorials/ex22.c (has a Fortran twin, ex22f.F) or 
  ex25.c in petsc-dev.
 
 
 How about recasting the problem as a DAE? The documentation seems to imply 
 that this is feasible. For ODE with nontrivial mass matrices such as arise 
 in FEM, the implicit/DAE interface significantly reduces overhead to prepare 
 the system for algebraic solvers (SNES/KSP) by having the user assemble the 
 correctly shifted matrix.
 
 Following ex15, solving
u_t = \Delta u
 can be recast as solving
F(t,u,\dot u) = 0
 with
F(t,u,v) = v-\Delta u.
 
 In this case, the IJacobian would be
K+aM,
 where K is the stiffness matrix (K_{i,j} = \int_\Omega \nabla \phi_i \cdot 
 \nabla \phi_j\,dx) and M the mass matrix (M_{i,j} = \int_\Omega \phi_i  
 \phi_j\,dx) .
 
 At the continuous level, the IFunction would be v-\Delta u, which cannot be 
 evaluated directly in the finite element basis by either solving
M F = M\dot u + Ku
 or using mass lumping.
 
 Am I expected to do this in my IFunction or is there a way to pass the mass 
 matrix to the TS?
 
 You should have some sense for what K is and whether it makes sense to 
 integrate your system implicitly or not. If you know that you want to use 
 fully implicit methods, then just dump everything into the IFunction. If part 
 of your system is non-stiff (and especially if it has less than C^1 
 regularity, or you need special properties from a certain explicit method), 
 then you can use the IMEX form G(X,Xdot) = F(X) with
 
 G(X,Xdot) := M*Xdot
 
 You can pass TSComputeIFunctionLinear and/or TSComputeIJacobianConstant to 
 TSSetIFunction() and TSSetIJacobian respectively if you have a linear 
 constant mass matrix.

Got it.


 As far as I understand, using IMPEX will lead to the same issue regardless of 
 the way I write my problem,  i.e. wether I write G(t,u,v) =  v-\Delta u and 
 F(t,u) = 0, or G(t,u,v) =  v and F(t,u) = \Delta u.
 
 As far as I can see, all examples and tests use structured meshes where the 
 mass matrix is the identity matrix.
 
 Identity as a mass matrix has nothing to do with structured vs. unstructured, 
 it's a matter of continuous finite elements versus finite difference/finite 
 volume/certain Petrov-Galerkin methods.

yep. Sorry, I meant finite differences, not structured mesh.


 
 

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] Questions about TS

2011-10-19 Thread Blaise Bourdin

On Oct 19, 2011, at 4:16 PM, Jed Brown wrote:

On Wed, Oct 19, 2011 at 16:08, Blaise Bourdin bourdin at lsu.edu wrote:
I get it now... My confusion was due to being used to derived all scheme on
the pde, then discretizing, whereas the documentation assumes that the
equation is already discretized. I should have figured it out.

For software purposes and sometimes also for analysis, the method of lines
approach is often useful. That's how the TS interfaces are set up.
That make sense. I work mostly on problems where whose time-continuous
formulation is obtained by writing a time discrete problem, then letting the
time discretization interval go to 0. I have too much of a tendency to think
this way.

Is it right to think of the division between ODE, DAE and IMEX in the
documentation as Fully Explicit vs. Fully implicit vs. Semi-implicit?

Sure, but it's a matter of the interface more than the method. You can write

Xdot = F(X)

and use -ts_type beuler to solve it fully implicitly. I would consider

G(X,Xdot) = F(X)

to be the most general interface. When an IMEX method is used, this has the
clear semantics that G is implicit and F is explicit. Explicit methods
usually assume G(X,Xdot) = Xdot which is the default if you never call
TSSetIFunction. I think we will eventually have support for using standard
explicit methods where you just put the mass matrix into G, but that isn't
done yet.

OK, this is clear. It will be really nice indeed when everything is reorganized
this way.

Thanks again,

Blaise

-- next part --
An HTML attachment was scrubbed...
URL:
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20111019/cab5bd66/attachment-0001.htm

[petsc-users] Questions about TS

2011-10-19 Thread Blaise Bourdin


On Oct 19, 2011, at 4:13 PM, Sean Farley wrote:

 Is it right to think of the division between ODE, DAE and IMEX in the 
 documentation as Fully Explicit vs. Fully implicit vs. Semi-implicit?
 
 Mostly, yes. Indeed, we made 'shortcuts' such as if the user provides no 
 IFunction, then treat it as udot = F(u,t) and if the user provides no RHS, 
 then treat it as G(udot, u, t) = 0.
That make sense, yes.

Thanks,
Blaise

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin







-- next part --
An HTML attachment was scrubbed...
URL: 
http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20111019/bfa08100/attachment.htm

[petsc-users] hdf5 directive

2011-07-26 Thread Blaise Bourdin

PETSC_HAVE_HDF5
On Jul 26, 2011, at 5:25 PM, Ataollah Mesgarnejad wrote:

 Dear all,
 
 Is there a preprocessor directive to check if PETSc has been compiled with 
 hdf5 or not ? (something like PETSC_HAS_HDF5?!!)
 
 Best,
 Ata

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] tecplot viewer

2011-05-25 Thread Blaise Bourdin

Hi,

Has anybody in this list ever written a tecplot viewer for DA based data? I 
have been fighting with tecplot hdf5 reader and it seems like it is not capable 
of reading vector valued fields.

Blaise

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] Storing an array in a PetscBag

2010-09-10 Thread Blaise Bourdin

Hi,

Is there an easy way to use PetscBag to manage arrays of user data?
optimally, I would like to something like
typedef struct {
PetscInt  n;
PetscReal *values;
} MyParameters;
MyParameters   *params

and be able to register values in the bag so that I can specify its values as
-n 3 -values 1.,2.,3. or -n 2 -values 1.,2.

Right now, I can see how I can get n prior to creating the bag and register 
params-values[0] to params-values[n-1] so that the command line could become 
-n 3 -values0 1 -values1 2 -values2 3 or -n 2 -values0 1 -values1 2

Would that require adding PetscBagRegisterXXXArray function or is it feasible 
in the current implementation?

Blaise 

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] [petsc4py] PETSc.Viewer().createHDF5 wipes out existing files

2010-09-09 Thread Blaise Bourdin

Hi,

Hopefully this does not get double posted. I sent the original email from the 
wrong account to the wrong list?

I am trying to do HDF5 IO in petsc4py. I noticed that when I create a viewer 
using 
 PETSc.Viewer().createHDF5(inputfile,comm= PETSc.COMM_WORLD)
the file outputfile is wiped out and replaced with a empty hdf5 container, 
which is bad since I am trying to read forrm it... 

Unsurprisingly, I have no problems using the longer approach:
 h5in = PETSc.Viewer().create(PETSc.COMM_WORLD)
 h5in.setType(PETSc.Viewer.Type.HDF5)
 h5in.setFileMode(PETSc.Viewer.Mode.READ)
 h5in.setFileName(inputfile)

Is this expected?

Regards,
Blaise


-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] VTKViewer with petsc4py

2010-07-31 Thread Blaise Bourdin


 On 30 July 2010 23:50, Blaise Bourdin bourdin at lsu.edu wrote:
 Great! that works.
 
 Now my problem is with DASetUniformCoordinates or DASetCoordinates. Do they 
 exist in petsc4py? I can't see them when I do a dir(da).
 
 
 No, they are not available. However, adding them should be really
 easy. Not sure iff I'll be able to do that this weekend. In any case,
 I'll have to ask you to hg clone http://path/to/petsc4py-dev; and hg
 update release-1.1. If you can beat me and add the missing calls,
 patches are going to be very welcome.
Here is the patch for DASetUniformCoordinates I don't really know what I am 
doing, but it seems reasonable. I'll work on the others this weekend.

diff -r 493264e77f7c src/PETSc/DA.pyx
--- a/src/PETSc/DA.pyx  Thu Jul 08 12:51:39 2010 -0300
+++ b/src/PETSc/DA.pyx  Fri Jul 30 22:58:33 2010 -0500
@@ -287,6 +287,8 @@ cdef class DA(Object):
 CHKERR( DALocalToLocalBegin(self.da, vl.vec, im, vlg.vec) )
 CHKERR( DALocalToLocalEnd  (self.da, vl.vec, im, vlg.vec) )
 
+def setUniformCoordinates(self, xmin, xmax, ymin, ymax, zmin, zmax):
+CHKERR( DASetUniformCoordinates(self.da, xmin, xmax, ymin, ymax, zmin, 
zmax) )
 #
 
 def getAO(self):
diff -r 493264e77f7c src/PETSc/petscda.pxi
--- a/src/PETSc/petscda.pxi Thu Jul 08 12:51:39 2010 -0300
+++ b/src/PETSc/petscda.pxi Fri Jul 30 22:58:33 2010 -0500
@@ -56,10 +56,10 @@ cdef extern from petscda.h nogil:
   PetscInt*,PetscInt*,PetscInt*,
   PetscInt*,PetscInt*,PetscInt*)
 
-#int DASetUniformCoordinates(PetscDA,
-#PetscReal,PetscReal,
-#PetscReal,PetscReal,
-#PetscReal,PetscReal)
+int DASetUniformCoordinates(PetscDA,
+PetscReal,PetscReal,
+PetscReal,PetscReal,
+PetscReal,PetscReal)
 #int DASetCoordinates(PetscDA,PetscVec)
 #int DAGetCoordinates(PetscDA,PetscVec*)
 #int DAGetGhostedCoordinates(PetscDA,PetscVec*)


 
 PS: The DA support is somewhat unfinished, I've never used DA's too much.

I don't need much. I just need to read a few tecplot files, do some basic 
post-processing and save them. What you have should be enough, once the 
coordinates stuff is added.


Blaise 



-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] VTKViewer with petsc4py

2010-07-30 Thread Blaise Bourdin

Hi,

Is PetscViewerSetFormat implemented in petsc4py (I am using petsc-3.1-p3 and 
petsc4py-1.1)? I am trying to save a DA and Vec in vtk format, but can't figure 
out the python equivalent to 
PetscViewerSetFormat(VTKViewer,PETSC_VIEWER_ASCII_VTK).

My C code would be 
  ierr = PetscViewerCreate(PETSC_COMM_WORLD,VTKViewer);CHKERRQ(ierr);
  ierr = PetscViewerSetType(VTKViewer,PETSC_VIEWER_ASCII);CHKERRQ(ierr);
  ierr = PetscViewerFileSetName(VTKViewer,filename);CHKERRQ(ierr);
  ierr = PetscViewerSetFormat(VTKViewer,PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);

  ierr = DAView(da,VTKViewer);CHKERRQ(ierr);
  ierr = VecView(U,VTKViewer);CHKERRQ(ierr);

What would be the python equivalent?

Thanks
Blaise

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-users] VTKViewer with petsc4py

2010-07-30 Thread Blaise Bourdin

Great! that works.

Now my problem is with DASetUniformCoordinates or DASetCoordinates. Do they 
exist in petsc4py? I can't see them when I do a dir(da).

Thanks,
Blaise 

On Jul 30, 2010, at 9:10 PM, Lisandro Dalcin wrote:

 On 30 July 2010 21:13, Blaise Bourdin bourdin at lsu.edu wrote:
 Hi,
 
 Is PetscViewerSetFormat implemented in petsc4py (I am using petsc-3.1-p3 and 
 petsc4py-1.1)?
 
 Viewer.setFormat() is certainly available.
 
 I am trying to save a DA and Vec in vtk format, but can't figure out the 
 python equivalent to
 PetscViewerSetFormat(VTKViewer,PETSC_VIEWER_ASCII_VTK).
 
 
 However, it seems I missed to add Viewer.Format.ASCII_VTK. As a
 workaround, you could do:
 
 viewer.setFormat(11) # 11 is the corresponding enum value from petscviewer.h
 
 Alternatively, you could monkeypatch like below at the beginning of your code:
 
 from petsc4py import PETSc
 PETSc.Viewer.Format.ASCII_VTK = 11
 
 
 My C code would be
  ierr = PetscViewerCreate(PETSC_COMM_WORLD,VTKViewer);CHKERRQ(ierr);
  ierr = PetscViewerSetType(VTKViewer,PETSC_VIEWER_ASCII);CHKERRQ(ierr);
  ierr = PetscViewerFileSetName(VTKViewer,filename);CHKERRQ(ierr);
  ierr = PetscViewerSetFormat(VTKViewer,PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
 
  ierr = DAView(da,VTKViewer);CHKERRQ(ierr);
  ierr = VecView(U,VTKViewer);CHKERRQ(ierr);
 
 What would be the python equivalent?
 
 
 vtkviewer = PETSc.Viewer().create(PETSc.COMM_WORLD)
 vtkviewer.setType('ascii') # or pass PETSc.Viewer.Type.ASCII
 vtkviewer.setFileName(filename)
 vtkviewer.setFromat(PETSc.Viewer.Format.ASCII_VTK) # of course,
 requires previous hack.
 
 
 However, this way is simpler:
 
 vtkviewer = PETSc.Viewer().createASCII(
filename, format=PETSc.Viewer.Format.ASCII_VTK,
comm= PETSc.COMM_WORLD)
 
 
 PS:  I'll add the missing format enumerate ASAP to petsc4py-dev,
 branch release-1-1. I'm near to make a micro release with a bunch of
 fixes, I?m just waiting for Cython 0.13 to be released.
 
 -- 
 Lisandro Dalcin
 ---
 CIMEC (INTEC/CONICET-UNL)
 Predio CONICET-Santa Fe
 Colectora RN 168 Km 472, Paraje El Pozo
 Tel: +54-342-4511594 (ext 1011)
 Tel/Fax: +54-342-4511169

-- 
Department of Mathematics and Center for Computation  Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

67 matches

Mail list logo