From da5ba8229990ceb78ab0bd0a8b8cfdc7d8c41f72 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 9 Mar 2022 20:40:03 +0100 Subject: [PATCH 1/3] initial condition can be specified per point/cell --- src/grid/discretization_grid.f90 | 48 ++++++++++++++++++++++++++++---- 1 file changed, 42 insertions(+), 6 deletions(-) diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 399bb180a..3dcb6446e 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -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) @@ -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) @@ -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) @@ -220,7 +220,7 @@ end function nodes0 !-------------------------------------------------------------------------------------------------- -!> @brief Calculate IP interface areas +!> @brief Calculate IP interface areas. !-------------------------------------------------------------------------------------------------- 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) @@ -257,7 +257,7 @@ end function cellSurfaceNormal !-------------------------------------------------------------------------------------------------- -!> @brief Build IP neighborhood relations +!> @brief Build IP neighborhood relations. !-------------------------------------------------------------------------------------------------- pure function IPneighborhood(cells) @@ -313,4 +313,40 @@ pure function IPneighborhood(cells) 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 From 790ca57ea0eb4e19e196cd329fac6b326c6774ac Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 9 Mar 2022 21:29:40 +0100 Subject: [PATCH 2/3] specify initial temperature per point --- src/grid/DAMASK_grid.f90 | 2 +- src/grid/discretization_grid.f90 | 13 +++++++------ src/grid/grid_thermal_spectral.f90 | 10 ++++------ 3 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index c1e89c88a..1946d31c8 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -308,7 +308,7 @@ program DAMASK_grid case (FIELD_THERMAL_ID) initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict) thermal => initial_conditions%get('thermal',defaultVal=emptyDict) - call grid_thermal_spectral_init(thermal%get_asFloat('T')) + call grid_thermal_spectral_init() case (FIELD_DAMAGE_ID) call grid_damage_spectral_init() diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 3dcb6446e..a32d35059 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -37,7 +37,8 @@ module discretization_grid size3offset !< (local) size offset in 3rd direction public :: & - discretization_grid_init + discretization_grid_init, & + discretization_grid_getInitialCondition contains @@ -316,10 +317,10 @@ end function IPneighborhood !-------------------------------------------------------------------------------------------------- !> @brief Read initial condition from VTI file. !-------------------------------------------------------------------------------------------------- -function getInitialCondition(label) +function discretization_grid_getInitialCondition(label) result(ic) character(len=*), intent(in) :: label - real(pReal), dimension(cells(1),cells(2),cells3) :: getInitialCondition + real(pReal), dimension(cells(1),cells(2),cells3) :: ic real(pReal), dimension(:), allocatable :: ic_global, ic_local integer(MPI_INTEGER_KIND) :: err_MPI @@ -328,7 +329,7 @@ function getInitialCondition(label) displs, sendcounts if (worldrank == 0) then - ic_global = VTI_read_Real(IO_read(interface_geomFile),label) + ic_global = VTI_readDataset_real(IO_read(interface_geomFile),label) else allocate(ic_global(0)) ! needed for IntelMPI endif @@ -345,8 +346,8 @@ function getInitialCondition(label) 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]) + ic = reshape(ic_local,[cells(1),cells(2),cells3]) -end function getInitialCondition +end function discretization_grid_getInitialCondition end module discretization_grid diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 1066a1f61..6bbeabf00 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -65,9 +65,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data !-------------------------------------------------------------------------------------------------- -subroutine grid_thermal_spectral_init(T_0) - - real(pReal), intent(in) :: T_0 +subroutine grid_thermal_spectral_init() PetscInt, dimension(0:worldsize-1) :: localK integer :: i, j, k, ce @@ -105,9 +103,9 @@ subroutine grid_thermal_spectral_init(T_0) !-------------------------------------------------------------------------------------------------- ! init fields - allocate(T_current(cells(1),cells(2),cells3), source=T_0) - allocate(T_lastInc(cells(1),cells(2),cells3), source=T_0) - allocate(T_stagInc(cells(1),cells(2),cells3), source=T_0) + T_current = discretization_grid_getInitialCondition('T') + T_lastInc = T_current + T_stagInc = T_current !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc From 8f9ba9d6e448024f399de3a7cf94138a3e3f351c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 10 Mar 2022 08:52:13 +0100 Subject: [PATCH 3/3] test with initial conditions in grid file --- PRIVATE | 2 +- src/grid/DAMASK_grid.f90 | 4 ---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/PRIVATE b/PRIVATE index 17250a3a2..9c1f83bab 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 17250a3a29f07af3c4d4e4083213d46efd534268 +Subproject commit 9c1f83babb7894bfaa16255d6c15a4a438c7f168 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index beb54127a..39870a8ff 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -107,8 +107,6 @@ program DAMASK_grid load_steps, & load_step, & solver, & - initial_conditions, & - thermal, & step_bc, & step_mech, & step_discretization @@ -303,8 +301,6 @@ program DAMASK_grid select case (ID(field)) case (FIELD_THERMAL_ID) - initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict) - thermal => initial_conditions%get('thermal',defaultVal=emptyDict) call grid_thermal_spectral_init() case (FIELD_DAMAGE_ID)