Commits

Richard Mills committed 0785ea8

Work on getting unstructured grids to use DMShell. This is not yet working (compiles but hits problems at runtime).

Comments (0)

Files changed (2)

src/pflotran/discretization.F90

 #include "finclude/petscdm.h"
 #include "finclude/petscdm.h90"
 #include "finclude/petscdmda.h"
+#include "finclude/petscdmshell.h90"
+       Interface
+        subroutine DMShellDefaultGlobalToLocalBegin(dm, g, mode, l ,ierr&
+     &)
+       DM dm ! DM
+       Vec g ! Vec
+       InsertMode mode ! InsertMode
+       Vec l ! Vec
+       integer ierr
+       end subroutine
+        subroutine DMShellDefaultGlobalToLocalEnd(dm, g, mode, l ,ierr)
+       DM dm ! DM
+       Vec g ! Vec
+       InsertMode mode ! InsertMode
+       Vec l ! Vec
+       integer ierr
+       end subroutine
+       End Interface
 
   type, public :: dm_ptr_type
-    DM :: sgdm  ! structured grid dm (PETSc DM)
-    type(ugdm_type), pointer :: ugdm ! unstructured grid dm
+    DM :: dm  ! PETSc DM
+    type(ugdm_type), pointer :: ugdm
+      ! Unstructured grid "private" dm.  This gets wrapped in a PETSc DM via 
+      ! DMShell routines.
   end type dm_ptr_type
 
   type, public :: discretization_type
   allocate(discretization%dm_1dof)
   allocate(discretization%dm_nflowdof)
   allocate(discretization%dm_ntrandof)
-  discretization%dm_1dof%sgdm = 0
-  discretization%dm_nflowdof%sgdm = 0
-  discretization%dm_ntrandof%sgdm = 0
+  discretization%dm_1dof%dm = 0
+  discretization%dm_nflowdof%dm = 0
+  discretization%dm_ntrandof%dm = 0
   nullify(discretization%dm_1dof%ugdm)
   nullify(discretization%dm_nflowdof%ugdm)
   nullify(discretization%dm_ntrandof%ugdm)
     case(STRUCTURED_GRID, STRUCTURED_GRID_MIMETIC)
       ! this function must be called to set up str_grid%lxs, etc.
       call StructGridComputeLocalBounds(discretization%grid%structured_grid, &
-                                        discretization%dm_1dof%sgdm,option)    
+                                        discretization%dm_1dof%dm,option)    
       discretization%grid%nlmax = discretization%grid%structured_grid%nlmax
       discretization%grid%ngmax = discretization%grid%structured_grid%ngmax
       discretization%grid%global_offset = &
   PetscInt :: ndof
   PetscInt :: stencil_width,stencil_type
   type(option_type) :: option
-  
+  PetscErrorCode :: ierr
+
   select case(discretization%itype)
     case(STRUCTURED_GRID, STRUCTURED_GRID_MIMETIC)
       call StructGridCreateDM(discretization%grid%structured_grid, &
-                                  dm_ptr%sgdm,ndof,stencil_width,stencil_type,option)
+                                  dm_ptr%dm,ndof,stencil_width,stencil_type,option)
     case(UNSTRUCTURED_GRID)
       call UGridCreateUGDM(discretization%grid%unstructured_grid, &
                            dm_ptr%ugdm,ndof,option)
+      call DMShellCreate(option%mycomm,dm_ptr%dm,ierr)
+      call DMShellSetGlobalToLocalVecScatter(dm_ptr%dm,dm_ptr%ugdm%scatter_gtol,ierr)
+      call DMShellSetGlobalToLocal(dm_ptr%dm,DMShellDefaultGlobalToLocalBegin,DMShellDefaultGlobalToLocalEnd,ierr)
   end select
 
 end subroutine DiscretizationCreateDM
     case(STRUCTURED_GRID, STRUCTURED_GRID_MIMETIC)
       select case (vector_type)
         case(GLOBAL)
-          call DMCreateGlobalVector(dm_ptr%sgdm,vector,ierr)
+          call DMCreateGlobalVector(dm_ptr%dm,vector,ierr)
         case(LOCAL)
-          call DMCreateLocalVector(dm_ptr%sgdm,vector,ierr)
+          call DMCreateLocalVector(dm_ptr%dm,vector,ierr)
         case(NATURAL)
-          call DMDACreateNaturalVector(dm_ptr%sgdm,vector,ierr)
+          call DMDACreateNaturalVector(dm_ptr%dm,vector,ierr)
       end select
     case(UNSTRUCTURED_GRID)
       call UGridDMCreateVector(discretization%grid%unstructured_grid, &
   select case(discretization%itype)
     case(STRUCTURED_GRID)
 #ifndef DMGET
-      call DMCreateMatrix(dm_ptr%sgdm,mat_type,Jacobian,ierr)
+      call DMCreateMatrix(dm_ptr%dm,mat_type,Jacobian,ierr)
 #else
-      call DMGetMatrix(dm_ptr%sgdm,mat_type,Jacobian,ierr)
+      call DMGetMatrix(dm_ptr%dm,mat_type,Jacobian,ierr)
 #endif
       call MatSetOption(Jacobian,MAT_KEEP_NONZERO_PATTERN,PETSC_FALSE,ierr)
       call MatSetOption(Jacobian,MAT_ROW_ORIENTED,PETSC_FALSE,ierr)
           call MatSetOption(Jacobian,MAT_ROW_ORIENTED,PETSC_FALSE,ierr)
         case(NTRANDOF)
 #ifndef DMGET
-          call DMCreateMatrix(dm_ptr%sgdm,mat_type,Jacobian,ierr)
+          call DMCreateMatrix(dm_ptr%dm,mat_type,Jacobian,ierr)
 #else
-          call DMGetMatrix(dm_ptr%sgdm,mat_type,Jacobian,ierr)
+          call DMGetMatrix(dm_ptr%dm,mat_type,Jacobian,ierr)
 #endif
           call MatSetOption(Jacobian,MAT_KEEP_NONZERO_PATTERN,PETSC_FALSE,ierr)
           call MatSetOption(Jacobian,MAT_ROW_ORIENTED,PETSC_FALSE,ierr)
     case(NFLOWDOF)
       allocate(discretization%dmc_nflowdof(mg_levels))
       do i=1, mg_levels
-        discretization%dmc_nflowdof(i)%sgdm = 0
+        discretization%dmc_nflowdof(i)%dm = 0
         nullify(discretization%dmc_nflowdof(i)%ugdm)
       enddo
       dmc_ptr => discretization%dmc_nflowdof
     case(NTRANDOF)
       allocate(discretization%dmc_ntrandof(mg_levels))
       do i=1, mg_levels
-        discretization%dmc_ntrandof(i)%sgdm = 0
+        discretization%dmc_ntrandof(i)%dm = 0
         nullify(discretization%dmc_ntrandof(i)%ugdm)
       enddo
       dmc_ptr => discretization%dmc_ntrandof
         if (i <= mg_levels - mg_levels_x ) refine_x = 1
         if (i <= mg_levels - mg_levels_y ) refine_y = 1
         if (i <= mg_levels - mg_levels_z ) refine_z = 1
-        call DMDASetRefinementFactor(dm_fine_ptr%sgdm, refine_x, refine_y, refine_z, &
+        call DMDASetRefinementFactor(dm_fine_ptr%dm, refine_x, refine_y, refine_z, &
                                    ierr)
-        call DMDASetInterpolationType(dm_fine_ptr%sgdm, DMDA_Q0, ierr)
-        call DMCoarsen(dm_fine_ptr%sgdm, option%mycomm, dmc_ptr(i)%sgdm, ierr)
+        call DMDASetInterpolationType(dm_fine_ptr%dm, DMDA_Q0, ierr)
+        call DMCoarsen(dm_fine_ptr%dm, option%mycomm, dmc_ptr(i)%dm, ierr)
 #ifndef DMGET
-        call DMCreateInterpolation(dmc_ptr(i)%sgdm, dm_fine_ptr%sgdm, &
+        call DMCreateInterpolation(dmc_ptr(i)%dm, dm_fine_ptr%dm, &
                                    interpolation(i), PETSC_NULL_OBJECT, ierr)
 #else
-        call DMGetInterpolation(dmc_ptr(i)%sgdm, dm_fine_ptr%sgdm, &
+        call DMGetInterpolation(dmc_ptr(i)%dm, dm_fine_ptr%dm, &
                                 interpolation(i), PETSC_NULL_OBJECT, ierr)
 #endif
         dm_fine_ptr => dmc_ptr(i)
   select case(discretization%itype)
     case(STRUCTURED_GRID)
 #ifndef DMGET
-      call DMCreateColoring(dm_ptr%sgdm,IS_COLORING_GLOBAL,MATBAIJ,coloring,&
+      call DMCreateColoring(dm_ptr%dm,IS_COLORING_GLOBAL,MATBAIJ,coloring,&
                             ierr)
 #else
-      call DMGetColoring(dm_ptr%sgdm,IS_COLORING_GLOBAL,MATBAIJ,coloring,ierr)
+      call DMGetColoring(dm_ptr%dm,IS_COLORING_GLOBAL,MATBAIJ,coloring,ierr)
 #endif
       ! I have set the above to use matrix type MATBAIJ, as that is what we 
       ! usually want (note: for DAs with 1 degree of freedom per grid cell, 
   
   dm_ptr => DiscretizationGetDMPtrFromIndex(discretization,dm_index)
     
-  select case(discretization%itype)
-    case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMGlobalToLocalBegin(dm_ptr%sgdm,global_vec,INSERT_VALUES,local_vec,ierr)
-      call DMGlobalToLocalEnd(dm_ptr%sgdm,global_vec,INSERT_VALUES,local_vec,ierr)
-    case(UNSTRUCTURED_GRID)
-      call VecScatterBegin(dm_ptr%ugdm%scatter_gtol,global_vec,local_vec, &
-                           INSERT_VALUES,SCATTER_FORWARD,ierr)
-      call VecScatterEnd(dm_ptr%ugdm%scatter_gtol,global_vec,local_vec, &
-                         INSERT_VALUES,SCATTER_FORWARD,ierr)
-  end select
+  call DMGlobalToLocalBegin(dm_ptr%dm,global_vec,INSERT_VALUES,local_vec,ierr)
+  call DMGlobalToLocalEnd(dm_ptr%dm,global_vec,INSERT_VALUES,local_vec,ierr)
   
 end subroutine DiscretizationGlobalToLocal
 
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMLocalToGlobalBegin(dm_ptr%sgdm,local_vec,INSERT_VALUES,global_vec,ierr)
-      call DMLocalToGlobalEnd(dm_ptr%sgdm,local_vec,INSERT_VALUES,global_vec,ierr)
+      call DMLocalToGlobalBegin(dm_ptr%dm,local_vec,INSERT_VALUES,global_vec,ierr)
+      call DMLocalToGlobalEnd(dm_ptr%dm,local_vec,INSERT_VALUES,global_vec,ierr)
    case(UNSTRUCTURED_GRID)
       call VecScatterBegin(dm_ptr%ugdm%scatter_ltog,local_vec,global_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMDALocalToLocalBegin(dm_ptr%sgdm,local_vec1,INSERT_VALUES,local_vec2,ierr)
-      call DMDALocalToLocalEnd(dm_ptr%sgdm,local_vec1,INSERT_VALUES,local_vec2,ierr)
+      call DMDALocalToLocalBegin(dm_ptr%dm,local_vec1,INSERT_VALUES,local_vec2,ierr)
+      call DMDALocalToLocalEnd(dm_ptr%dm,local_vec1,INSERT_VALUES,local_vec2,ierr)
     case(UNSTRUCTURED_GRID)
       call VecScatterBegin(dm_ptr%ugdm%scatter_ltol,local_vec1,local_vec2, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMDAGlobalToNaturalBegin(dm_ptr%sgdm,global_vec,INSERT_VALUES,natural_vec,ierr)
-      call DMDAGlobalToNaturalEnd(dm_ptr%sgdm,global_vec,INSERT_VALUES,natural_vec,ierr)
+      call DMDAGlobalToNaturalBegin(dm_ptr%dm,global_vec,INSERT_VALUES,natural_vec,ierr)
+      call DMDAGlobalToNaturalEnd(dm_ptr%dm,global_vec,INSERT_VALUES,natural_vec,ierr)
     case(UNSTRUCTURED_GRID)
       call VecScatterBegin(dm_ptr%ugdm%scatter_gton,global_vec,natural_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMDANaturalToGlobalBegin(dm_ptr%sgdm,natural_vec,INSERT_VALUES,global_vec,ierr)
-      call DMDANaturalToGlobalEnd(dm_ptr%sgdm,natural_vec,INSERT_VALUES,global_vec,ierr)
+      call DMDANaturalToGlobalBegin(dm_ptr%dm,natural_vec,INSERT_VALUES,global_vec,ierr)
+      call DMDANaturalToGlobalEnd(dm_ptr%dm,natural_vec,INSERT_VALUES,global_vec,ierr)
     case(UNSTRUCTURED_GRID)
       call VecScatterBegin(dm_ptr%ugdm%scatter_ntog,natural_vec,global_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMGlobalToLocalBegin(dm_ptr%sgdm,global_vec,INSERT_VALUES,local_vec,ierr)
+      call DMGlobalToLocalBegin(dm_ptr%dm,global_vec,INSERT_VALUES,local_vec,ierr)
     case(UNSTRUCTURED_GRID)
       call VecScatterBegin(dm_ptr%ugdm%scatter_gtol,global_vec,local_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMGlobalToLocalEnd(dm_ptr%sgdm,global_vec,INSERT_VALUES,local_vec,ierr)
+      call DMGlobalToLocalEnd(dm_ptr%dm,global_vec,INSERT_VALUES,local_vec,ierr)
     case(UNSTRUCTURED_GRID)
       call VecScatterEnd(dm_ptr%ugdm%scatter_gtol,global_vec,local_vec, &
                          INSERT_VALUES,SCATTER_FORWARD,ierr)
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMDALocalToLocalBegin(dm_ptr%sgdm,local_vec1,INSERT_VALUES,local_vec2,ierr)
+      call DMDALocalToLocalBegin(dm_ptr%dm,local_vec1,INSERT_VALUES,local_vec2,ierr)
     case(UNSTRUCTURED_GRID)
       call VecScatterBegin(dm_ptr%ugdm%scatter_ltol,local_vec1,local_vec2, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMDALocalToLocalEnd(dm_ptr%sgdm,local_vec1,INSERT_VALUES,local_vec2,ierr)
+      call DMDALocalToLocalEnd(dm_ptr%dm,local_vec1,INSERT_VALUES,local_vec2,ierr)
     case(UNSTRUCTURED_GRID)
       call VecScatterEnd(dm_ptr%ugdm%scatter_ltol,local_vec1,local_vec2, &
                          INSERT_VALUES,SCATTER_FORWARD,ierr)    
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMDAGlobalToNaturalBegin(dm_ptr%sgdm,global_vec,INSERT_VALUES,natural_vec,ierr)
+      call DMDAGlobalToNaturalBegin(dm_ptr%dm,global_vec,INSERT_VALUES,natural_vec,ierr)
     case(UNSTRUCTURED_GRID)
       call VecScatterBegin(dm_ptr%ugdm%scatter_gton,global_vec,natural_vec, &
                            INSERT_VALUES,SCATTER_FORWARD,ierr)
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMDAGlobalToNaturalEnd(dm_ptr%sgdm,global_vec,INSERT_VALUES,natural_vec,ierr)
+      call DMDAGlobalToNaturalEnd(dm_ptr%dm,global_vec,INSERT_VALUES,natural_vec,ierr)
     case(UNSTRUCTURED_GRID)
       call VecScatterEnd(dm_ptr%ugdm%scatter_gton,global_vec,natural_vec, &
                          INSERT_VALUES,SCATTER_FORWARD,ierr)       
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMDANaturalToGlobalBegin(dm_ptr%sgdm,natural_vec,INSERT_VALUES,global_vec,ierr)
+      call DMDANaturalToGlobalBegin(dm_ptr%dm,natural_vec,INSERT_VALUES,global_vec,ierr)
     case(UNSTRUCTURED_GRID)
   end select
   
   
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      call DMDANaturalToGlobalEnd(dm_ptr%sgdm,natural_vec,INSERT_VALUES,global_vec,ierr)
+      call DMDANaturalToGlobalEnd(dm_ptr%dm,natural_vec,INSERT_VALUES,global_vec,ierr)
     case(UNSTRUCTURED_GRID)
   end select
   
       
   select case(discretization%itype)
     case(STRUCTURED_GRID,STRUCTURED_GRID_MIMETIC)
-      if (discretization%dm_1dof%sgdm /= 0) &
-        call DMDestroy(discretization%dm_1dof%sgdm,ierr)
-      discretization%dm_1dof%sgdm = 0
-      if (discretization%dm_nflowdof%sgdm /= 0) &
-        call DMDestroy(discretization%dm_nflowdof%sgdm,ierr)
-      discretization%dm_nflowdof%sgdm = 0
-      if (discretization%dm_ntrandof%sgdm /= 0) &
-        call DMDestroy(discretization%dm_ntrandof%sgdm,ierr)
-      discretization%dm_ntrandof%sgdm = 0
+      if (discretization%dm_1dof%dm /= 0) &
+        call DMDestroy(discretization%dm_1dof%dm,ierr)
+      discretization%dm_1dof%dm = 0
+      if (discretization%dm_nflowdof%dm /= 0) &
+        call DMDestroy(discretization%dm_nflowdof%dm,ierr)
+      discretization%dm_nflowdof%dm = 0
+      if (discretization%dm_ntrandof%dm /= 0) &
+        call DMDestroy(discretization%dm_ntrandof%dm,ierr)
+      discretization%dm_ntrandof%dm = 0
       if (associated(discretization%dmc_nflowdof)) then
         do i=1,size(discretization%dmc_nflowdof)
-          call DMDestroy(discretization%dmc_nflowdof(i)%sgdm,ierr)
+          call DMDestroy(discretization%dmc_nflowdof(i)%dm,ierr)
         enddo
         deallocate(discretization%dmc_nflowdof)
         nullify(discretization%dmc_nflowdof)
       endif
       if (associated(discretization%dmc_ntrandof)) then
         do i=1,size(discretization%dmc_ntrandof)
-          call DMDestroy(discretization%dmc_ntrandof(i)%sgdm,ierr)
+          call DMDestroy(discretization%dmc_ntrandof(i)%dm,ierr)
         enddo
         deallocate(discretization%dmc_ntrandof)
         nullify(discretization%dmc_ntrandof)

src/pflotran/realization.F90

     case(STRUCTURED_GRID, STRUCTURED_GRID_MIMETIC)
       grid => discretization%grid
       ! set up nG2L, nL2G, etc.
-      call GridMapIndices(grid, discretization%dm_1dof%sgdm, &
+      call GridMapIndices(grid, discretization%dm_1dof%dm, &
                           discretization%stencil_type,&
                           discretization%lsm_flux_method, &
                           option)
         call StructGridCreateTVDGhosts(grid%structured_grid, &
                                        realization%reaction%naqcomp, &
                                        field%tran_xx, &
-                                       discretization%dm_1dof%sgdm, &
+                                       discretization%dm_1dof%dm, &
                                        field%tvd_ghosts, &
                                        discretization%tvd_ghost_scatter, &
                                        option)
 
 
      dm_ptr => DiscretizationGetDMPtrFromIndex(discretization, NFLOWDOF)
-     call GridComputeGlobalCell2FaceConnectivity(grid, discretization%MFD, dm_ptr%sgdm, NFLOWDOF, option)
+     call GridComputeGlobalCell2FaceConnectivity(grid, discretization%MFD, dm_ptr%dm, NFLOWDOF, option)
   
    end if
 #endif