initial condition can be specified per point/cell

This commit is contained in:
Martin Diehl 2022-03-09 20:40:03 +01:00
parent 2b27388c05
commit da5ba82299
1 changed files with 42 additions and 6 deletions

View File

@ -43,7 +43,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief reads the geometry file to obtain information on discretization !> @brief Read the geometry file to obtain information on discretization.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine discretization_grid_init(restart) subroutine discretization_grid_init(restart)
@ -171,7 +171,7 @@ end subroutine discretization_grid_init
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element) !> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element).
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
function IPcoordinates0(cells,geomSize,cells3Offset) function IPcoordinates0(cells,geomSize,cells3Offset)
@ -196,7 +196,7 @@ end function IPcoordinates0
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief Calculate position of undeformed nodes (pretend to be an element) !> @brief Calculate position of undeformed nodes (pretend to be an element).
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure function nodes0(cells,geomSize,cells3Offset) pure function nodes0(cells,geomSize,cells3Offset)
@ -220,7 +220,7 @@ end function nodes0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate IP interface areas !> @brief Calculate IP interface areas.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function cellSurfaceArea(geomSize,cells) pure function cellSurfaceArea(geomSize,cells)
@ -238,7 +238,7 @@ end function cellSurfaceArea
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate IP interface areas normals !> @brief Calculate IP interface areas normals.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function cellSurfaceNormal(nElems) pure function cellSurfaceNormal(nElems)
@ -257,7 +257,7 @@ end function cellSurfaceNormal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Build IP neighborhood relations !> @brief Build IP neighborhood relations.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function IPneighborhood(cells) pure function IPneighborhood(cells)
@ -313,4 +313,40 @@ pure function IPneighborhood(cells)
end function IPneighborhood end function IPneighborhood
!--------------------------------------------------------------------------------------------------
!> @brief Read initial condition from VTI file.
!--------------------------------------------------------------------------------------------------
function getInitialCondition(label)
character(len=*), intent(in) :: label
real(pReal), dimension(cells(1),cells(2),cells3) :: getInitialCondition
real(pReal), dimension(:), allocatable :: ic_global, ic_local
integer(MPI_INTEGER_KIND) :: err_MPI
integer, dimension(worldsize) :: &
displs, sendcounts
if (worldrank == 0) then
ic_global = VTI_read_Real(IO_read(interface_geomFile),label)
else
allocate(ic_global(0)) ! needed for IntelMPI
endif
call MPI_Gather(product(cells(1:2))*cells3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,&
1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Gather(product(cells(1:2))*cells3, 1_MPI_INTEGER_KIND,MPI_INTEGER,sendcounts,&
1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
allocate(ic_local(product(cells(1:2))*cells3))
call MPI_Scatterv(ic_global,sendcounts,displs,MPI_DOUBLE,ic_local,size(ic_local),&
MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
getInitialCondition = reshape(ic_local,[cells(1),cells(2),cells3])
end function getInitialCondition
end module discretization_grid end module discretization_grid