Dear all,
I am using TS to solve the PDE
u_t + c * u_c = eps * u_xx in 1D i.e. smthg. similar to
https://urldefense.us/v3/__https://petsc.org/release/src/ts/tutorials/ex4.c.html__;!!G_uCfscf7eWS!eeErx_FNg50N1eRlD9eyfoEp8I0llTQmVFpPWoz5_LQhL8K6uOId5zbgle0ghE71NJmon8T_vwwfN3LnEbf4cTUXtlkWJRflikU$
except that I am using DMPlex, instead of a structured grid.
Since the problem is linear, I tried smthg. like:
PetscCall(RHSMatrixHeat(ts,0.0,u,A,A,&appctx));
PetscCall(TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx));
PetscCall(TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx));
My problem is that the A matrix (as well as the global vector) is
missing the rows/cols corresponding to the endpoints where Dirichlet
boundary conditions are prescribed;
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.
therefore, A*u = f(u) except in the two vertices adjacent to the
endpoints, i.e. first and last entries of the global vector.
The code is attached; when the flag -use_jacobian is set, f(u) is
computed as A*u, if it is not set, f(u) is computed in FormRHSFunction.
The latter approach works, the former does not.
Regards,
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!eeErx_FNg50N1eRlD9eyfoEp8I0llTQmVFpPWoz5_LQhL8K6uOId5zbgle0ghE71NJmon8T_vwwfN3LnEbf4cTUXtlkWzhmq6Yg$
program main
! solve the PDE: u_t + c * u_x = eps * u_xx on a 1D plex
! using TS
!
! if L = 1 is the length of the box, Pe = c*L/eps
!
!
! 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_simplex true -dm_plex_interpolate -dm_plex_check_all -dm_petscsection_view
! -ts_max_steps 500 -options_left -Peclet xxxx
!
! if the option -Peclet is omitted, then solves u_tt = u_xx
! if the option -use_jacobian is used, compute f(u) as A*u
!
! Author: aldo.bonfigli...@unibas.it
!
#include <petsc/finclude/petsc.h>
use petsc
implicit none
DM :: dm
TS :: ts
Vec :: u,y,rhs
Mat :: jac,Pmat
PetscBool :: lflag
external FormRHSFunction,dtFunction,MyOwnTSMonitoringRoutine,TSSetRHSJacobian
PetscErrorCode :: ierr
PetscReal:: user(5)
PetscReal:: dt
PetscReal:: Peclet,rnorm_2,rnorm_max,s
PetscInt:: nitems
character(len=PETSC_MAX_PATH_LEN) :: IOBuffer
PetscViewer :: viewer
logical, parameter :: verbose = .FALSE.
PetscCallA(PetscInitialize(ierr))
!
! Make the mesh (the mesh length is returned in user(3))
!
call Createmesh(dm,user(3))
!
PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER , "-Peclet", Peclet, lflag, ierr))
!
IF(.NOT.lflag)then
Peclet = 0.d0
user(1) = 1.d0 ! diffusion coefficient
user(2) = 0.d0 ! convective speed
else
user(2) = 1.d0 ! convective speed
user(1) = user(2) * user(3) / Peclet ! diffusion coefficient
endif
write(IOBuffer,FMT=*)"\n Peclet number is: ",Peclet,"\n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
write(IOBuffer,FMT=*)"\n diffusion coefficient is: ",user(1),"\n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
!
! Is the following call really needed ?!?!?!
!
PetscCallA(DMPlexCreateClosureIndex(dm, PETSC_NULL_SECTION, ierr))
!
! Attach a section to the DM
!
! call SetupSection(dm)
call CreateSectionAlternate(dm)
!
! u is the vector where the solution is stored
!
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(DMCreateMatrix(dm, jac, ierr))
PetscCallA(DMCreateMatrix(dm, Pmat, 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)
! Build the matrix
PetscCallA(TSRHSJacobianFn(ts, dt, u, jac, Pmat, user, ierr))
!
! This is a linear problem, i.e. u_t = F(u) = A*u
!
! https://petsc.org/release/src/ts/tutorials/ex3.c.html
!
PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-use_jacobian',lflag, ierr))
if(lflag)then
PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,TSComputeRHSFunctionLinear,user,ierr))
PetscCallA(TSSetRHSJacobian(ts,jac,jac,TSComputeRHSJacobianConstant,user,ierr))
else
! Set the user-defined function for computing the residual
PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr))
! PetscCallA(TSSetRHSFunction(ts,r,FormRHSFunction,user,ierr))
!
! Sets the function to compute the Jacobian of G, where U_t = G(U,t), as well as the location to store the matrix.
PetscCallA(TSSetRHSJacobian(ts, jac, Pmat, TSRHSJacobianFn, user, ierr))
endif
!
! looks like TSPseudoSetTimeStep is NOT available in fortran
! PetscCallA(TSPseudoSetTimeStep(ts, dtFunction, PETSC_NULL_VEC,ierr))
!
! PetscCallA(TSMonitorSet(ts, MyOwnTSMonitoringRoutine, PETSC_NULL_REAL, PETSC_NULL_FUNCTION, ierr))
PetscCallA(TSMonitorSet(ts, MyOwnTSMonitoringRoutine, user, PETSC_NULL_FUNCTION, 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))
!
! Here we check if rhs = jac * u
!
! +++++++++++++++ start checking tha f(u) = A*u ++++++++++++++++
!
if(verbose)then
write(IOBuffer,FMT="(A)")"This is u \n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
PetscCallA(VecView(u,PETSC_VIEWER_STDOUT_SELF,ierr))
endif
!
PetscCallA(VecDuplicate(u,rhs,ierr))
PetscCallA(VecSet(rhs,0.d0,ierr))
call FormRHSFunction(ts,dt,u,rhs,user,ierr)
if(verbose)then
write(IOBuffer,FMT="(A)")"This is f(u) \n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
PetscCallA(VecView(rhs,PETSC_VIEWER_STDOUT_SELF,ierr))
endif
!
PetscCallA(VecDuplicate(u,y,ierr))
PetscCallA(MatMult(jac,u,y,ierr))
if(verbose)then
write(IOBuffer,FMT="(A)")"This is A*u \n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
PetscCallA(VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr))
endif
!
PetscCallA(MatView(jac,PETSC_VIEWER_STDOUT_SELF,ierr))
!
PetscCallA(VecAXPY(y, -1.d0, rhs, ierr))
write(IOBuffer,FMT=*)"This is A*u-f(u): should be 0.d0 \n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
PetscCallA(VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr))
!
PetscCallA(VecDestroy(y,ierr))
PetscCallA(VecDestroy(rhs,ierr))
!
! +++++++++++++++ finished checking ++++++++++++++++
!
! 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))
! Need X support
! PetscCallA(PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,PETSC_NULL_CHARACTER,0,0,300,300,viewer,ierr))
!
PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD,viewer,ierr))
PetscCallA(PetscViewerSetType(viewer,PETSCVIEWERASCII,ierr))
PetscCallA(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE,ierr))
PetscCallA(PetscViewerFileSetName(viewer,"sol.dat",ierr))
! View the vector
PetscCallA(VecView(u,viewer,ierr))
! Deallocate memory.
PetscCallA(PetscViewerDestroy(viewer,ierr))
!
!
PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD,viewer,ierr))
PetscCallA(PetscViewerSetType(viewer,PETSCVIEWERASCII,ierr))
PetscCallA(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE,ierr))
PetscCallA(PetscViewerFileSetName(viewer,"exact.dat",ierr))
!
PetscCallA(VecDuplicate(u,y,ierr))
call FormInitialSolution(ts,y,user,ierr) ! returns the exact solution except when Peclet = 0
! View the vector
PetscCallA(VecView(y,viewer,ierr))
! Deallocate memory.
PetscCallA(PetscViewerDestroy(viewer,ierr))
!
! Steady Error calculation
!
PetscCallA(PetscViewerCreate(PETSC_COMM_WORLD,viewer,ierr))
PetscCallA(PetscViewerSetType(viewer,PETSCVIEWERASCII,ierr))
PetscCallA(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE,ierr))
PetscCallA(PetscViewerFileSetName(viewer,"err.dat",ierr))
!
PetscCallA(VecGetSize(y, nitems, ierr))
s = 1.d0/dble(nitems)
s = user(4) ! average mesh spacing
PetscCallA(VecAXPY(y, -1.d0, u, ierr))
PetscCallA(VecView(y,viewer,ierr))
! Deallocate memory.
PetscCallA(PetscViewerDestroy(viewer,ierr))
!
PetscCallA(VecNorm(y, NORM_2, rnorm_2, ierr))
rnorm_2 = dsqrt(s) * rnorm_2
PetscCallA(VecNorm(y, NORM_MAX, rnorm_max, ierr))
write(IOBuffer,FMT=*)"h;L2;Lmax error norms are : ",user(4),rnorm_2,rnorm_max,"\n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
!
!
!
PetscCallA(VecDestroy(y,ierr))
!
! Free work space. All PETSc objects should be destroyed when they
! are no longer needed.
!
PetscCallA(MatDestroy(jac,ierr))
PetscCallA(MatDestroy(Pmat,ierr))
PetscCallA(VecDestroy(u,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,length)
#include <petsc/finclude/petsc.h>
use petsc
implicit none
DM :: dm
PetscReal :: length
PetscErrorCode :: ierr
PetscInt :: depth,ndim,ibgn,iend
character(len=PETSC_MAX_PATH_LEN) :: IOBuffer
!
! 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.
!
!
PetscCall(DMPlexCreate(PETSC_COMM_WORLD, dm, ierr))
PetscCall(PetscObjectSetName(dm, '1D plex', ierr))
PetscCall(DMSetFromOptions(dm,ierr))
PetscCall(DMGetDimension(dm,ndim,ierr))
if( ndim .NE. 1 )call PetscError(PETSC_COMM_WORLD,1,PETSC_ERROR_INITIAL,'Topological dimension must equal 1')
PetscCall(DMViewFromOptions(dm, PETSC_NULL_VEC, '-dm_view', ierr))
PetscCall(DMView(dm,PETSC_VIEWER_STDOUT_WORLD,ierr))
PetscCall(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
PetscCall(DMPlexGetDepthStratum(dm, depth, ibgn , iend,ierr))
write(IOBuffer,FMT=*)"there are ",iend-ibgn," entries in depth: ",depth," \n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
enddo
! can I get the length from Petsc?
length = 1.d0 ! should be set to the length
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(*)
PetscErrorCode :: ierr
DM :: dm
PetscInt :: depth,cell,cbgn,cend
PetscReal, pointer :: pcentroid(:)
PetscReal, pointer :: pnormal(:)
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 :: eps,c,length,cplus,cmins
PetscReal :: du,ux,vol,h
pcentroid => centroid
pnormal => normal
p_res => res
eps = user(1) ! viscosity coefficient
c = user(2) ! convective speed
length = user(3)
cplus = 0.5d0*(c+dabs(c))
cmins = c - cplus
h = user(4)
PetscCall(TSGetDM(ts, dm, ierr));
PetscCall(DMGetLocalVector(dm, localu, ierr ))
PetscCall(DMGetLocalVector(dm, localf, ierr ))
PetscCall(DMGlobaltoLocal(dm, u, INSERT_VALUES, localu, ierr ))
PetscCall(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
PetscCall(DMPlexGetDepthStratum(dm, depth, cbgn, cend, ierr))
do cell = cbgn,cend-1 ! loop over 1D cells (edges)
PetscCall(DMPlexVecGetClosure(dm, PETSC_NULL_SECTION, localu, cell, p_values, ierr))
PetscCall(DMPlexComputeCellGeometryFVM(dm, cell, vol, pcentroid, pnormal, ierr))
du = (p_values(2)-p_values(1))
ux = du/vol ! cell-wise gradient
p_res(1) = (-cmins*du+eps*ux)!/h ! send it to the "left"
p_res(2) = (-cplus*du-eps*ux)!/h ! send it to the "right"
! write(6,*)cell,(p_values(j),j=1,2),du
PetscCall(DMPlexVecRestoreClosure(dm, PETSC_NULL_SECTION, localu, cell, p_values, ierr))
PetscCall(DMPlexVecSetClosure(dm, PETSC_NULL_SECTION, localf, cell, p_res, ADD_VALUES, ierr ))
enddo
! PetscCall(VecScale(localu,1.d0/h, ierr ))
PetscCall(DMRestoreLocalVector(dm, localu, ierr ))
! Insert (do NOT ADD!!) values into the "global" solution vector f
PetscCall(DMLocalToGlobal(dm, localf, INSERT_VALUES, f, ierr))
PetscCall(DMRestoreLocalVector(dm, localf, ierr ))
! do I need to VecAssembly or is it done inside DMLocalToGlobalEnd ?
! PetscCall(VecAssemblyBegin(f,ierr))
! PetscCall(VecAssemblyEnd(f,ierr))
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(*)
PetscErrorCode :: ierr
DM :: dm
PetscBool :: simplex
PetscInt :: size,depth,vbgn,vend
PetscReal :: x
integer :: i,nitems,ifail
PetscScalar, pointer :: xx_g(:),xx_u(:)
PetscReal :: eps,c,length,Pe,rx,fi0,fin
Vec :: coordVec
character(len=PETSC_MAX_PATH_LEN) :: IOBuffer
logical, parameter :: verbose = .FALSE.
eps = user(1)
c = user(2)
length = user(3)
Pe = c*length/eps
rx = 1.d0/length
fi0 = 1.d0
fin = 0.d0
PetscCall(TSGetDM(ts, dm, ierr));
PetscCall(DMPlexIsSimplex(dm, simplex, ierr));
if( simplex .NEQV. PETSC_TRUE )then
PetscCall(PetscError(PETSC_COMM_WORLD,1,PETSC_ERROR_INITIAL,"I can only work with simplices"))
endif
! use DMGetCoordinatesLocal or it will NOT work in the parallel case
PetscCall(DMGetCoordinatesLocal(dm,coordVec,ierr))
PetscCall(VecGetSize(coordVec,size,ierr))
! call VecView(coordVec,PETSC_VIEWER_STDOUT_WORLD, ierr)
!
! do I need the CoordinateSection ?
!
! PetscCall(DMGetCoordinateSection(dm, coordSection,ierr))
!
! I need to retrieve the LocalVector or Dirichlet nodes will be missing
PetscCall(DMGetLocalVector(dm, localu, ierr ))
PetscCall(DMGlobalToLocal(dm, u, INSERT_VALUES, localu, 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 = 0 ! pick up the vertices
PetscCall(DMPlexGetDepthStratum(dm, depth, vbgn, vend, ierr))
nitems = vend-vbgn
if(nitems.NE.size)then
write(IOBuffer,FMT=*)"There is a mismatch ",nitems,size,"\n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
ifail = nitems-size
else
ifail = 0
endif
PetscCall(VecGetArrayF90(localu,xx_u,ierr))
PetscCall(VecGetArrayReadF90(coordVec,xx_g,ierr))
do i = 1,nitems ! loop over the vertices
x = xx_g(i)
if(verbose)then
write(IOBuffer,FMT=*)"x-Coord of gridpoint ",i," is ",x," \n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
endif
if(pe .LE. 1.d-8)then ! diffusive limit
xx_u(i) = fi0 + (x*x*x) * (fin-fi0)
else
xx_u(i) = fi0 + (exp(pe*x*rx)-1.d0)/(exp(pe)-1.d0) * (fin-fi0)
endif
enddo
PetscCall(VecRestoreArrayReadF90(coordVec,xx_g,ierr))
PetscCall(VecRestoreArrayF90(localu,xx_u,ierr))
! call VecView(localu,PETSC_VIEWER_STDOUT_WORLD, ierr)
PetscCall(DMLocalToGlobal(dm, localu, INSERT_VALUES, u, ierr))
PetscCall(DMRestoreLocalVector(dm, localu, ierr ))
if(ifail.NE.0)call PetscError(PETSC_COMM_WORLD,1,PETSC_ERROR_INITIAL,'There is a mismatch')
return
end
subroutine SetupSection(dm)
#include <petsc/finclude/petsc.h>
use petsc
implicit none
DM :: dm
PetscInt :: ndofs
PetscErrorCode :: ierr
PetscSection :: s
PetscInt :: depth, vStart, vEnd, v, numDof, numFixed, nitems , numPoints, i
DMLabel :: label
PetscBool :: hasLabel
! PetscInt, dimension(:),pointer :: p_idx
PetscInt,dimension(:),pointer :: constraints
PetscInt, pointer :: xx_v(:)
IS, target, dimension(1) :: bcPointIS
ndofs = 1
numDof = 0
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
!
PetscCall(DMHasLabel(dm, "marker", hasLabel, ierr))
if(.NOT. hasLabel) call CreateBCLabel(dm, "marker")
PetscCall(DMGetLabel(dm, "marker", label, ierr))
PetscCall(DMGetStratumSize(dm,'marker',1,numPoints,ierr))
PetscCall(DMLabelGetStratumIS(label, 1, bcPointIS(1), ierr))
PetscCall(ISGetSize(bcPointIS(1), nitems, ierr))
write(6,*)"IS bcPoint has ",nitems," items "
PetscCall(ISGetIndicesF90(bcPointIS(1),xx_v,ierr))
do i = 1,numPoints
v = xx_v(i)
write(6,*)"Point ",v,"is a boundary point"
PetscCall(PetscSectionSetConstraintDof(s, v, ndofs, ierr))
! PetscCall(PetscSectionSetConstraintIndicesF90(s,pointID(1),constraints,ierr))
enddo
PetscCall(PetscSectionSetUp(s,ierr))
do v = vStart, vEnd-1
PetscCall(PetscSectionGetDof(s, v, numDof, ierr))
PetscCall(PetscSectionGetConstraintDof(s, v, numFixed, ierr))
write(6,*)"Point ",v," has ",numDof, " dofs; ",numFixed," of these are fixed"
enddo
allocate(constraints(1))
! allocate(pointID(1))
constraints(1) = 0
do i = 1,numPoints
v = xx_v(i)
write(6,*)"Point ",v,"is a boundary point",constraints
! PetscCall(PetscSectionSetConstraintIndicesF90(s,xx_v(i),constraints(1),ierr))
PetscCall(PetscSectionSetConstraintIndicesF90(s,v,constraints,ierr))
enddo
deallocate(constraints)
! deallocate(pointID)
! look @ ./src/dm/impls/plex/tests/ex98f90.F90
! for the use of PetscSectionSetConstraintIndicesF90 and PetscSectionSetFieldConstraintIndicesF90
! allocate(p_idx(1))
!
! am I doing smthg wrong when calling PetscSectionGetConstraintIndicesF90 ?
!
! do v = vStart, vEnd-1
! PetscCall(PetscSectionGetConstraintDof(s, v, numFixed, ierr))
! if(numFixed.NE.0)then
! PetscCall(PetscSectionGetConstraintIndicesF90(s, v, p_idx, ierr))
! write(6,*)v,(p_idx(j),j=1,numFixed)
! PetscCall(PetscSectionRestoreConstraintIndicesF90(s, v, p_idx, ierr))
! endif
! enddo
10 continue
! deallocate(p_idx)
! PetscCall(PetscSectionSetup(s, ierr))
! Get the nof points in depth=1 that are marked; should be 2
! do i = 1,numPoints
! PetscCall(ISGetIndicesF90(bcpointIS,pointID,ierr))
! PetscCall(PetscSectionSetConstraintIndicesF90(s,pointID(1),constraints,ierr))
! enddo
PetscCall(DMSetLocalSection(dm, s,ierr))
PetscCall(PetscSectionGetStorageSize(s, nitems, ierr))
write(6,*)"Constrained section has size = ",nitems
PetscCall(PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD, ierr))
PetscCall(PetscSectionDestroy(s,ierr))
return
end
!
subroutine dtFunction(ts,newdt,user,ierr)
!
! -------------------- Evaluate time-step ---------------------
!
#include <petsc/finclude/petsc.h>
use petsc
implicit none
TS :: ts
PetscReal::newdt
PetscReal::user(*)
PetscErrorCode :: ierr
DM :: dm
PetscInt :: depth,cell,cbgn,cend,nitems
PetscReal :: eps,c
PetscReal :: vol,h,dt(2)
PetscReal, pointer :: pcentroid(:)
PetscReal, pointer :: pnormal(:)
PetscReal, target, dimension(1) :: centroid
PetscReal, target, dimension(1) :: normal
PetscReal, parameter :: cfl = 0.8d0
character(len=PETSC_MAX_PATH_LEN) :: IOBuffer
pcentroid => centroid
pnormal => normal
eps = user(1) ! user(1)
c = user(2) ! user(2)
h = 0.d0 ! average mesh spacing
PetscCall(TSGetDM(ts, dm, ierr))
depth = 1 ! loop over cells (edges), NOT vertices !!
PetscCall(DMPlexGetDepthStratum(dm, depth, cbgn , cend,ierr))
do cell = cbgn,cend-1 ! loop over 1D cells (edges)
PetscCall(DMPlexComputeCellGeometryFVM(dm, cell, vol, pcentroid, pnormal,ierr))
h = h + vol
enddo
nitems = cend-cbgn
h = h/dble(nitems) ! average cell volume
dt(1) = h/c
dt(2) = h*h/(2.d0*eps)
newdt = min(dt(1),dt(2)) ! maximum explicit timestep
newdt = cfl*newdt
write(IOBuffer,FMT=*)"\n Average mesh spacing is ",h," \n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
write(IOBuffer,FMT=*)"\n Time-step is set to ",dt,newdt," \n"
PetscCallA(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
user(4) = h
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
subroutine CreateSectionAlternate(dm)
#include "petsc/finclude/petscdmplex.h"
#include "petsc/finclude/petscdmlabel.h"
use petscdmplex
use petscsys
implicit none
DM :: dm
DMLabel :: label
PetscSection :: section
PetscInt :: ndim,numFields,numBC,size
PetscInt :: i
DMLabel, pointer :: nolabel(:) => NULL()
PetscInt, target, dimension(3) :: numComp
PetscInt, pointer :: pNumComp(:)
PetscInt, target, dimension(12) :: numDof
PetscInt, pointer :: pNumDof(:)
PetscInt, target, dimension(1) :: bcField
PetscInt, pointer :: pBcField(:)
PetscInt, parameter :: zero = 0, one = 1, two = 2
PetscInt :: nitems
! PetscMPIInt ::
IS, target, dimension(1) :: bcCompIS
IS, target, dimension(1) :: bcPointIS
IS, pointer :: pBcCompIS(:)
IS, pointer :: pBcPointIS(:)
PetscErrorCode :: ierr
character(len=PETSC_MAX_PATH_LEN) :: IOBuffer
PetscCall(DMGetDimension(dm, ndim, ierr))
! Create a scalar field u, a vector field v, and a surface vector field w
numFields = 1
numComp(1) = 1
! numComp(2) = ndim
! numComp(3) = ndim-1
pNumComp => numComp
do i = 1, numFields*(ndim+1)
numDof(i) = 0
end do
! Let u be defined on vertices
numDof(0*(ndim+1)+1) = 1
! Let v be defined on cells
! numDof(1*(ndim+1)+ndim+1) = ndim
! Let v be defined on faces
! numDof(2*(ndim+1)+ndim) = ndim-1
do i = 1, numFields*(ndim+1)
write(6,*)"numDof(",i,")= ",numDof(i)
end do
pNumDof => numDof
! Prescribe a Dirichlet condition on u on the boundary
numBC = 1
! Get the nof points in depth=1 that are marked; should be 2
PetscCall(DMGetStratumSize(dm, "marker", one, size, ierr))
write(IOBuffer,FMT=*)"DMGetStratumSize found ",size," points \n"
PetscCall(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
PetscCall(DMGetLabel(dm, "marker", label, ierr))
PetscCall(DMLabelGetStratumIS(label, one, bcPointIS(1), ierr))
PetscCall(ISGetSize(bcPointIS(1), nitems, ierr))
write(IOBuffer,FMT=*)"IS bcPointIS has ",nitems," items \n"
PetscCall(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
! bcField is An array of size numBC giving the field number for each boundary condition
bcField(1) = 0
pBcField => bcField
! PetscErrorCode ISCreateStride(MPI_Comm comm, PetscInt n, PetscInt first, PetscInt step, IS *is)
PetscCall(ISCreateStride(PETSC_COMM_WORLD, one, zero, one, bcCompIS(1), ierr))
PetscCall(ISGetSize(bcCompIS(1), nitems, ierr))
write(IOBuffer,FMT=*)"IS bcCompIS has ",nitems," items \n"
PetscCall(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
! bcCompIS is An array of size numBC giving an IS holding the field components to which each boundary condition applies
pBcCompIS => bcCompIS
pBcPointIS => bcPointIS
! Create a PetscSection with this data layout
PetscCall(DMSetNumFields(dm, numFields,ierr))
PetscCall(DMPlexCreateSection(dm,nolabel,pNumComp,pNumDof,numBC,pBcField,pBcCompIS,pBcPointIS,PETSC_NULL_IS,section,ierr))
PetscCall(ISDestroy(bcCompIS(1), ierr))
PetscCall(ISDestroy(bcPointIS(1), ierr))
! Name the Field variables
PetscCall(PetscSectionSetFieldName(section, zero, 'u', ierr))
PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD, ierr))
! Tell the DM to use this data layout
PetscCall(DMSetLocalSection(dm, section, ierr))
! Cleanup
PetscCall(PetscSectionDestroy(section, ierr))
return
end
!
!
!
subroutine MyOwnTSMonitoringRoutine(ts,iter,time,u,user,ierr)
!
!
!
#include <petsc/finclude/petsc.h>
use petsc
implicit none
TS :: ts
PetscInt :: iter
PetscReal :: time
Vec :: u
PetscReal :: user(*)
PetscErrorCode :: ierr
PetscReal :: s,t
Vec :: f
character(len=PETSC_MAX_PATH_LEN) :: IOBuffer
PetscCall(VecDuplicate(u,f,ierr))
PetscCall(VecNorm(u,NORM_2,t,ierr))
PetscCall(FormRHSFunction(ts,time,u,f,user,ierr))
PetscCall(VecNorm(f,NORM_2,s,ierr))
PetscCall(VecDestroy(f,ierr))
write(IOBuffer,FMT="(A,1X,I6,1X,A,F10.6,2(1X,A,E12.6),A)")"@ iter = ",iter," time = ",time," ||f|| = ",s," ||u|| = ",t,"\n"
PetscCall(PetscPrintf(PETSC_COMM_WORLD, IOBuffer, ierr))
return
end
subroutine TSRHSJacobianFn(ts, time, u, jac, Pmat, user, ierr)
#include <petsc/finclude/petsc.h>
use petsc
implicit none
TS :: ts
PetscReal :: time
Vec :: u
DM :: dm
Mat :: jac,Pmat
PetscReal :: user(*)
!
PetscInt :: ierr,height,cbgn,cend,cell,numDof,point
PetscInt, dimension(:),pointer :: p_points
PetscBool :: useCone = PETSC_TRUE
! PetscBool :: useCone = PETSC_FALSE
integer :: j,ic
PetscReal, target, dimension(2) :: centroid
PetscReal, target, dimension(2) :: normal
PetscReal, pointer :: pcentroid(:)
PetscReal, pointer :: pnormal(:)
PetscInt :: offset
PetscInt :: row(2)
PetscInt :: col(2)
PetscReal :: valA(2,2)
PetscSection :: section
PetscReal :: aplus,amins,eps,vol
pcentroid => centroid
pnormal => normal
aplus = 0.5d0*(user(2)+dabs(user(2)))
amins = user(2) - aplus
eps = user(1)
PetscCall(TSGetDM(ts, dm, ierr))
height = 1
PetscCall(DMGetLocalSection(dm, section, ierr))
PetscCall(DMPlexGetDepthStratum(dm, height, cbgn , cend,ierr))
do cell = cbgn,cend-1
!
PetscCall(DMPlexComputeCellGeometryFVM(dm, cell, vol, pcentroid, pnormal,ierr))
! write(6,*)"Cell # ",cell," ranges btw. : ",(values(j),j=1,2),vol,(pcentroid(j),j=1,1),(pnormal(j),j=1,1)
!
! cell-wise stiffness matrix
! Sorted, row-oriented input will generally assemble the fastest. The default is row-oriented.
!
valA(2,1)=-amins+eps/vol
valA(1,2)=+aplus+eps/vol
valA(1,1)=-valA(1,2)
valA(2,2)=-valA(2,1)
!
PetscCall(DMPlexGetTransitiveClosure(dm, cell, useCone, p_points, ierr))
! write(6,*)"Cell # ",cell," is bounded by vertices: ",(p_points(j),j=1,6,1)
ic = 0
! loop over the (two) vertices of the 1D cell
do j = 3,5,2 ! this is tricky, because vertices are stored in locations 3 and 5 of p_points
ic = ic+1
point = p_points(j)
!
! we know there is only one dof
!
PetscCall(PetscSectionGetDof(section, point, numDof, ierr))
PetscCall(PetscSectionGetOffset(section, point, offset,ierr))
row(ic) = offset
col(ic) = offset
! write(6,*)point,offset,numDof
enddo
PetscCall(MatSetValuesLocal(jac,2,row,2,col,valA(1,1),ADD_VALUES,ierr))
PetscCall(DMPlexRestoreTransitiveClosure(dm, cell,useCone, p_points, ierr))
enddo
PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
return
end
-dm_plex_dim 1
-Peclet 50
-dm_plex_shape box
-dm_plex_box_faces 30
-dm_plex_simplex true
-dm_plex_interpolate
-dm_plex_check_all
-dm_petscsection_view
-ts_max_steps 5000
-options_left
