Dear users,
I have been using petsc's KSP for over 20 years and I am considering
using DMPLEX to replace my own data structure in a mixed FEM/FVM CFD code.
To do so, I am trying to understand DMPLEX by writing a 1D code that
solves u_t = u_xx using PSEUDOTS.
In order to prescribe Dirichlet BCs at the endpoints of the 1D box, I
use PetscSectionSetConstraintDof when building the PetscSection.
A global vector is missing both the shared dofs which are not owned by
this process, as well as /constrained/ dofs. These constraints
represent essential (Dirichlet) boundary conditions. They are dofs
that have a given fixed value, so they are present in local vectors
for assembly purposes, but absent from global vectors since they are
never solved for during algebraic solves.
My global Vec has indeed two entries less than the local one.
When initializing the solution or evaluating the rhs function I transfer
data from the global to local representation, do the calculation, then
transfer back.
I am doing something wrong, though, which shows up in the attached log file.
My simple code is also attached.
Thank you for your advice.
Aldo
--
Dr. Aldo Bonfiglioli
Associate professor of Fluid Machines
Scuola di Ingegneria
Universita' della Basilicata
V.le dell'Ateneo lucano, 10 85100 Potenza ITALY
tel:+39.0971.205203 fax:+39.0971.205215
web:https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?m=002423__;!!G_uCfscf7eWS!cWgDpzKW9KGEkqxu775gDU488CBQ5itLl2l_HUsDKcPstbcuSm0_bZPHkMciXajTi3539fuMPe8TD1YZuirDfmSBMSXXxrQZWns$
DM Object: box 1 MPI process
type: plex
box in 1 dimension:
Number of 0-cells per rank: 11
Number of 1-cells per rank: 10
Labels:
marker: 1 strata with value/size (1 (2))
Face Sets: 2 strata with value/size (1 (1), 2 (1))
depth: 2 strata with value/size (0 (11), 1 (10))
celltype: 2 strata with value/size (0 (11), 1 (10))
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Point 10: Global dof 0 != 1 size - number of constraints
[0]PETSC ERROR: WARNING! There are unused option(s) set! Could be the program
crashed before usage or a spelling mistake, etc!
[0]PETSC ERROR: Option left: name:-ts_max_steps value: 10 source: command line
[0]PETSC ERROR: Option left: name:-ts_monitor (no value) source: command line
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.22.0, Sep 28, 2024
[0]PETSC ERROR: ./ex15f90 with 1 MPI process(es) and PETSC_ARCH linux_gnu on
DESKTOP-04AR6V6.lan by abonfi Mon Nov 4 19:28:00 2024
[0]PETSC ERROR: Configure options: --with-mpi=1
--with-mpi-dir=/usr/lib64/mpi/gcc/mpich --with-shared-libraries=1
--with-debugging=1 --with-blas-lib=/usr/lib64/libblas.a
--with-lapack-lib=/usr/lib64/liblapack.a --with-triangle=1
--download-triangle=yes
[0]PETSC ERROR: #1 PetscSFSetGraphSection() at
/home/abonfi/src/petsc-3.22.0/src/vec/is/sf/utils/sfutils.c:179
[0]PETSC ERROR: #2 DMCreateSectionSF() at
/home/abonfi/src/petsc-3.22.0/src/dm/interface/dm.c:4783
[0]PETSC ERROR: #3 DMGetSectionSF() at
/home/abonfi/src/petsc-3.22.0/src/dm/interface/dm.c:4722
[0]PETSC ERROR: #4 DMGlobalToLocalBegin() at
/home/abonfi/src/petsc-3.22.0/src/dm/interface/dm.c:2851
[0]PETSC ERROR: #5 DMGlobalToLocal() at
/home/abonfi/src/petsc-3.22.0/src/dm/interface/dm.c:2812
[0]PETSC ERROR: #6 ex15f90.F90:280
Abort(62) on node 0 (rank 0 in comm 16): application called
MPI_Abort(MPI_COMM_SELF, 62) - process 0
there are 11 entries in depth: 0
there are 10 entries in depth: 1
Setting up a section on 11 meshpoints with 1 dofs
The DM is marked
Point 10 is a boundary point
Point 20 is a boundary point
Array u has type : seq
(Global) Array u has size : 9
coordVec has size 11
program main
! solve u_t + c u_x = eps * u_xx on a 1D plex
! using pseudo-time stepping
!
! for the time being c = 0
!
! good references are:
!
! https://petsc.org/release/src/ts/tutorials/ex4.c.html
! and
! https://petsc.org/main/src/ts/tutorials/ex17.c.html
!
! run with:
! ./ex15f90 -dm_plex_dim 1 -dm_plex_shape box -dm_plex_box_faces 30 -dm_plex_interpolate -dm_plex_check_all -ts_monitor -ts_max_steps 1
!
! Author:
!
#include <petsc/finclude/petsc.h>
use petsc
implicit none
DM :: dm
TS :: ts
Vec :: u ,r
external FormRHSFunction,dtFunction
PetscInt :: size,ndim
PetscInt :: ndofs = 1
! PetscViewer :: viewer
PetscErrorCode :: ierr
PetscReal::user(3)
PetscReal::dt
PetscReal::s
VecType :: type
! PetscViewer :: myfile
user(1) = 0.01d0 ! viscosity coefficient
user(2) = 0.0d0 ! convective speed
PetscCallA(PetscInitialize(ierr))
!
! Make the mesh
!
call Createmesh(dm)
PetscCallA(DMPlexCreateClosureIndex(dm, PETSC_NULL_SECTION, ierr))
!
! Attach a section to the DM
!
call SetupSection(dm,ndofs)
PetscCallA(DMCreateGlobalVector(dm, u, ierr))
! A global vector is missing both the shared dofs which are not owned by this process, as well as constrained dofs.
PetscCallA(VecSet(u,0.d0,ierr))
PetscCallA(VecGetType(u,type,ierr))
write(6,*)"Array u has type : ",type
PetscCallA(VecGetSize(u,size,ierr))
write(6,*)"(Global) Array u has size : ",size
PetscCallA(VecDuplicate(u, r, ierr))
! look @ https://petsc.org/release/src/ts/utils/dmplexlandau/tutorials/ex1f90.F90.html
PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
PetscCall(TSSetProblemType(ts, TS_LINEAR, ierr));
PetscCallA(TSSetDM(ts, dm, ierr))
! Initialize the solution
call FormInitialSolution(ts,u,user,ierr)
PetscCallA(VecGetSize(u,size,ierr))
PetscCallA(VecNorm(u,NORM_2,s,ierr))
write(6,*)"(Global) Array u has size : ",size," ||u|| = ",s
!
! Set the user-defined function for computing the residual
PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr))
!
! looks like TSPseudoSetTimeStep is NOT available in fortran
! PetscCallA(TSPseudoSetTimeStep(ts, dtFunction, PETSC_NULL_VEC,ierr))
!
! check also the following:
! DMPlexInsertBoundaryValues
! DMPlexMarkBoundaryFaces
!
! PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx));
! PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx));
! PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx));
! PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER,ierr))
PetscCallA(TSSetFromOptions(ts,ierr))
! User-defined function to compute the maximum allowable timestep;
! the sintax complies with TSPseudoSetTimeStep, but there seems to be no fortran interface for TSPseudoSetTimeStep
call dtFunction(ts,dt,user,ierr)
PetscCallA(TSSetTimeStep(ts, dt,ierr))
! call exit(-1)
! DMTSCheckFromOptions does NOT seem to be available in fortran
! PetscCallA(DMTSCheckFromOptions(ts, u, ierr))
PetscCallA(TSSolve(ts, u, ierr))
PetscCallA(VecViewFromOptions(u, PETSC_NULL_VIEWER, "-sol_vec_view", ierr))
!
! Deallocate memory.
!
PetscCallA(VecDestroy(u,ierr))
! PetscCallA(VecDestroy(r,ierr))
PetscCallA(TSDestroy(ts,ierr))
PetscCallA(DMDestroy(dm,ierr))
PetscCallA(PetscFinalize(ierr))
!
! /home/abonfi/src/petsc-3.22.0/src/dm/impls/plex/tutorials ex1F90
!
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Normal end of execution.'
write ( *, '(a)' ) ' '
stop
end
!
!
!
subroutine Createmesh(dm)
#include <petsc/finclude/petsc.h>
use petsc
implicit none
DM :: dm
PetscErrorCode :: ierr
PetscInt :: depth,ndim,ibgn,iend
! The support/closure of a point is defined as the set of all point of the next highest topological dimension that is connected to the current point.
! For example, in a 3 dimensional fully interpolated mesh, the closure of a vertex would be the set of all of the edges which contain this vertex.
!
! The cone of a point is defined as the set of every points of the next lowest topological dimension. For example, in a 3 dimensional fully interpolated mesh,
! the cone of a cell would be the set of every face on belonging to that cell.
!
PetscCallA(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr))
PetscCallA(PetscObjectSetName(dm, '1D plex', ierr))
PetscCallA(DMSetFromOptions(dm,ierr))
PetscCallA(DMGetDimension(dm,ndim,ierr))
if( ndim .NE. 1 )call PetscError(PETSC_COMM_WORLD,1,PETSC_ERROR_INITIAL,'Topological dimension must equal 1')
PetscCallA(DMViewFromOptions(dm, PETSC_NULL_VEC, '-dm_view', ierr))
PetscCallA(DMView(dm,PETSC_VIEWER_STDOUT_WORLD,ierr))
PetscCallA(DMSetup(dm,ierr))
! Depth indexing is related to topological dimension. Depth stratum 0 contains the lowest topological dimension points, often “vertices”.
! If the mesh is “interpolated” (see DMPlexInterpolate()), then depth stratum 1 contains the next higher dimension, e.g., “edges”.
do depth = 0,1
PetscCallA(DMPlexGetDepthStratum(dm, depth, ibgn , iend,ierr))
write(6,*)"there are ",iend-ibgn," entries in depth: ",depth
enddo
return
end
!
! -------------------- Evaluate Function F(x) ---------------------
!
subroutine FormRHSFunction(ts,t,u,F,user,ierr)
#include <petsc/finclude/petsc.h>
use petsc
implicit none
TS :: ts
PetscReal::t
Vec :: u,F
Vec :: localu,localF
PetscReal::user(3)
PetscErrorCode :: ierr
DM :: dm
PetscInt :: depth,cell,cbgn,cend
PetscScalar, dimension(:),pointer :: p_values
PetscScalar, dimension(:),pointer :: p_res
PetscReal, target, dimension(2) :: centroid
PetscReal, target, dimension(2) :: normal
PetscReal, target, dimension(2) :: res
PetscReal, pointer :: pcentroid(:)
PetscReal, pointer :: pnormal(:)
PetscReal :: eps,c
PetscReal :: du,ux,vol,s
integer :: j
pcentroid => centroid
pnormal => normal
p_res => res
eps = user(1) ! viscosity coefficient
c = user(2) ! convective speed
c = 0.d0
PetscCallA(TSGetDM(ts, dm, ierr));
PetscCallA(DMGetLocalVector(dm, localu, ierr ))
PetscCallA(DMGetLocalVector(dm, localf, ierr ))
PetscCallA(DMGlobaltoLocal(dm, u, INSERT_VALUES, localu, ierr ))
PetscCallA(VecSet(localf, 0.d0, ierr ))
!
! Depth indexing is related to topological dimension. Depth stratum 0 contains the lowest topological dimension points, often “vertices”.
! If the mesh is “interpolated” (see DMPlexInterpolate()), then depth stratum 1 contains the next higher dimension, e.g., “edges”.
!
depth = 1 ! pick up cells (edges), NOT vertices
PetscCallA(DMPlexGetDepthStratum(dm, depth, cbgn, cend, ierr))
do cell = cbgn,cend-1 ! loop over 1D cells (edges)
PetscCallA(DMPlexVecGetClosure(dm, PETSC_NULL_SECTION, localu,cell,p_values,ierr))
PetscCallA(DMPlexComputeCellGeometryFVM(dm, cell, vol, pcentroid, pnormal,ierr))
du = (p_values(2)-p_values(1))
ux = du/vol ! cell-wise gradient
p_res(1) = (-(0.5d0*(c+dabs(c))*du)+eps*ux)/vol ! send it to the "left"
p_res(2) = (-(0.5d0*(c-dabs(c))*du)-eps*ux)/vol ! send it to the "right"
PetscCallA(DMPlexVecRestoreClosure(dm, PETSC_NULL_SECTION, localu,cell,p_values,ierr))
PetscCallA(DMPlexVecSetClosure(dm, PETSC_NULL_SECTION, localf, cell, p_res, ADD_VALUES, ierr ))
! write(6,*)(p_res(j),j=1,2)
! write(6,*)"Cell # ",cell," is bounded by vertices: ",(p_points(j),j=3,5,2)
enddo
PetscCallA(DMRestoreLocalVector(dm, localu, ierr ))
PetscCallA(DMLocalToGlobal(dm, localf, ADD_VALUES, f, ierr))
PetscCallA(DMRestoreLocalVector(dm, localf, ierr ))
! Insert values into the "global" solution vector f
! do I need to VecAssembly or is it done inside DMLocalToGlobalEnd ?
! PetscCallA(VecAssemblyBegin(f,ierr))
! PetscCallA(VecAssemblyEnd(f,ierr))
PetscCallA(VecNorm(f,NORM_2,s,ierr))
! call VecView(f,PETSC_VIEWER_STDOUT_WORLD, ierr)
write(6,*)"Cells are ",cend-cbgn," function norm is : ",s
return
end
!
!
!
subroutine FormInitialSolution(ts,u,user,ierr)
!
! build the initial solution: we only need to initialize processor owned vertices, I believe ......
!
#include <petsc/finclude/petsc.h>
use petsc
implicit none
TS :: ts
Vec :: u
Vec :: localu
PetscReal::user(3)
PetscErrorCode :: ierr
DM :: dm
PetscBool :: simplex
PetscInt :: size,depth,vbgn,vend
PetscReal :: x
integer :: i,nitems,ifail
PetscScalar, pointer :: xx_g(:),xx_u(:)
PetscReal, parameter :: c = -30.d0
PetscReal, parameter :: delt = -0.5d0
Vec :: coordVec
PetscCallA(TSGetDM(ts, dm, ierr));
PetscCall(DMPlexIsSimplex(dm, simplex, ierr));
if( simplex .NEQV. PETSC_TRUE )then
call PetscError(PETSC_COMM_WORLD,1,PETSC_ERROR_INITIAL,"I can only work with simplices")
endif
! Should I use DMGetCoordinatesLocal instead ????
PetscCallA(DMGetCoordinates(dm,coordVec,ierr))
PetscCallA(VecGetSize(coordVec,size,ierr))
write(6,*)"coordVec has size ",size
! call VecView(coordVec,PETSC_VIEWER_STDOUT_WORLD, ierr)
!
! do I need the CoordinateSection ?
!
! PetscCallA(DMGetCoordinateSection(dm, coordSection,ierr))
!
! I need to retrieve the LocalVector or Dirichlet nodes will be missing
PetscCallA(DMGetLocalVector(dm, localu, ierr ))
PetscCallA(DMGlobalToLocal(dm, u, INSERT_VALUES, localu, ierr))
PetscCallA(VecGetSize(localu,size,ierr))
write(6,*)"localu has (global) size ",size
PetscCallA(VecGetLocalSize(localu,size,ierr))
write(6,*)"localu has (local) size ",size
!
! Depth indexing is related to topological dimension. Depth stratum 0 contains the lowest topological dimension points, often “vertices”.
! If the mesh is “interpolated” (see DMPlexInterpolate()), then depth stratum 1 contains the next higher dimension, e.g., “edges”.
!
depth = 0 ! pick up the vertices
PetscCallA(DMPlexGetDepthStratum(dm, depth, vbgn, vend, ierr))
nitems = vend-vbgn
if(nitems.NE.size)then
write(6,*)"There is a mismatch ",nitems,size
ifail = nitems-size
else
ifail = 0
endif
call VecGetArrayF90(localu,xx_u,ierr)
call VecGetArrayReadF90(coordVec,xx_g,ierr)
do i = 1,nitems ! loop over the vertices
x = xx_g(i)
! write(6,*)"x @ ",i," is ",x
! r = dsqrt((x - .5) * (x - .5))
! if (r .LT. .125d0) then
! xx_u(i) = exp(c * r * r * r)
xx_u(i) = 1.d0 + delt * x
! else
! xx_u(i) = 0.d0
! endif
enddo
call VecRestoreArrayReadF90(coordVec,xx_g,ierr)
call VecRestoreArrayF90(localu,xx_u,ierr)
call VecView(localu,PETSC_VIEWER_STDOUT_WORLD, ierr)
PetscCallA(DMLocalToGlobal(dm, localu, INSERT_VALUES, u, ierr))
PetscCallA(DMRestoreLocalVector(dm, localu, ierr ))
! write(6,*)ifail
if(ifail.NE.0)call exit(ifail)
return
end
subroutine SetupSection(dm,ndofs)
#include <petsc/finclude/petsc.h>
use petsc
implicit none
DM :: dm
PetscInt :: ndofs
PetscErrorCode :: ierr
PetscSection :: s
PetscInt :: depth, vStart, vEnd, v, numDof,iflag
DMLabel :: label
PetscBool :: hasLabel
numDof = 1
depth = 0 ! loop over vertices
PetscCall(DMPlexGetDepthStratum(dm, depth, vStart, vEnd, ierr))
PetscCall(PetscSectionCreate(PETSC_COMM_WORLD, s, ierr))
! PetscCall(PetscSectionSetNumFields(s, 1, ierr))
! PetscCall(PetscSectionSetFieldComponents(s, 0, 1, ierr))
PetscCall(PetscSectionSetChart(s, vStart, vEnd, ierr))
write(6,*)"Setting up a section on ",vEnd-vStart," meshpoints with ",ndofs," dofs"
do v = vStart, vEnd-1
PetscCall(PetscSectionSetDof(s, v, ndofs, ierr))
! PetscCall(PetscSectionSetFieldDof(s, v, 0, 1));
enddo
! un-comment the following goto to avoid setting BCs
! goto 10
!
! Identify boundary vertices
!
PetscCallA(DMHasLabel(dm, "marker", hasLabel, ierr))
if(.NOT.hasLabel)then
call CreateBCLabel(dm, "marker")
else
PetscCallA(DMGetLabel(dm, "marker", label, ierr))
! PetscCallA(DMLabelGetValue(DMLabel label, v, ilabel, ierr))
write(6,*)"The DM is marked" ! and its label is ",label
endif
do v = vStart, vEnd-1
PetscCallA(DMGetLabelValue(dm, "marker", v, iflag, ierr))
! PetscCallA(DMLabelGetValue(DMLabel label, v, ilabel, ierr))
if( iflag .NE. -1 )then
write(6,*)"Point ",v,"is a boundary point"
PetscCallA(PetscSectionSetConstraintDof(s, v, numDof, ierr))
endif
enddo
10 continue
PetscCall(PetscSectionSetUp(s,ierr))
PetscCall(DMSetLocalSection(dm, s,ierr))
PetscCall(PetscSectionDestroy(s,ierr))
return
end
!
! -------------------- Evaluate time-step ---------------------
!
subroutine dtFunction(ts,newdt,user,ierr)
#include <petsc/finclude/petsc.h>
use petsc
implicit none
TS :: ts
PetscReal::newdt
PetscReal::user(3)
PetscErrorCode :: ierr
DM :: dm
PetscInt :: depth,cell,cbgn,cend,nitems
PetscReal :: eps
PetscReal :: vol,h
PetscReal, pointer :: pcentroid(:)
PetscReal, pointer :: pnormal(:)
PetscReal, parameter :: cfl = 0.8d0
eps = user(1) ! user(1)
h = 0.d0 ! average mesh spacing
PetscCallA(TSGetDM(ts, dm, ierr));
depth = 1 ! loop over cells (edges), NOT vertices !!
PetscCallA(DMPlexGetDepthStratum(dm, depth, cbgn , cend,ierr))
do cell = cbgn,cend-1 ! loop over 1D cells (edges)
PetscCallA(DMPlexComputeCellGeometryFVM(dm, cell, vol, pcentroid, pnormal,ierr))
h = h + vol
enddo
nitems = cend-cbgn
h = h/dble(nitems) ! average cell volume
newdt = h*h/(2.d0*eps) ! maximum explicit timestep
newdt = cfl*newdt
return
end
!
subroutine CreateBCLabel(dm, name)
#include <petsc/finclude/petsc.h>
use petsc
implicit none
DM :: dm
character*(*) name
DMLabel :: label
PetscErrorCode :: ierr
PetscCall(DMCreateLabel(dm, name, ierr))
PetscCall(DMGetLabel(dm, name, label, ierr))
! PetscCall(DMConvert(dm, DMPLEX, plex));
! PetscCall(DMPlexMarkBoundaryFaces(plex, 1, label, ierr))
PetscCall(DMPlexMarkBoundaryFaces(dm, 1, label, ierr))
! PetscCall(DMDestroy(&plex, ierr));
return
end
