diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index f3e74e030..8e3884851 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -36,21 +36,21 @@ variables: # Names of module files to load # =============================================================================================== # ++++++++++++ Compiler +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - COMPILER_GNU: "Compiler/GNU/10" + COMPILER_GNU: "Compiler/GNU/10" COMPILER_INTELLLVM: "Compiler/oneAPI/2022.0.1 Libraries/IMKL/2022.0.1" COMPILER_INTEL: "Compiler/Intel/2022.0.1 Libraries/IMKL/2022.0.1" # ++++++++++++ MPI ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.1" - MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0" - MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0" + MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.2" + MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0" + MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0" # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1" + PETSC_GNU: "Libraries/PETSc/3.16.4/GNU-10-OpenMPI-4.1.2" PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0" - PETSC_INTEL: "Libraries/PETSc/3.16.3/Intel-2022.0.1-IntelMPI-2021.5.0" + PETSC_INTEL: "Libraries/PETSc/3.16.4/Intel-2022.0.1-IntelMPI-2021.5.0" # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - MSC: "FEM/MSC/2021.3.1" - IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020" - HDF5Marc: "HDF5/1.12.1/Intel-19.1.2" + MSC: "FEM/MSC/2021.3.1" + IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020" + HDF5Marc: "HDF5/1.12.1/Intel-19.1.2" ################################################################################################### diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c7b129f0..b4c405319 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,7 +42,7 @@ string(TOUPPER "${CMAKE_BUILD_TYPE}" CMAKE_BUILD_TYPE) if(CMAKE_BUILD_TYPE STREQUAL "DEBUG" OR CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") set(DEBUG_FLAGS "${DEBUG_FLAGS} -DDEBUG") set(PARALLEL "OFF") - set(OPTI "OFF") + set(OPTI "DEBUG") elseif(CMAKE_BUILD_TYPE STREQUAL "RELEASE") set(PARALLEL "ON") set(OPTI "DEFENSIVE") diff --git a/PRIVATE b/PRIVATE index ebb7f0ce7..5774122bf 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit ebb7f0ce78d11275020af0ba60f929f95b446932 +Subproject commit 5774122bf48d637704bb4afb10b87c34a4dbcaba diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index cafcfabca..c204d1af2 100644 --- a/cmake/Compiler-GNU.cmake +++ b/cmake/Compiler-GNU.cmake @@ -9,26 +9,24 @@ if (OPENMP) set (OPENMP_FLAGS "-fopenmp") endif () -if (OPTIMIZATION STREQUAL "OFF") +if (OPTIMIZATION STREQUAL "DEBUG") + set (OPTIMIZATION_FLAGS "-Og") +elseif (OPTIMIZATION STREQUAL "OFF") set (OPTIMIZATION_FLAGS "-O0") elseif (OPTIMIZATION STREQUAL "DEFENSIVE") - set (OPTIMIZATION_FLAGS "-O2 -mtune=generic") + set (OPTIMIZATION_FLAGS "-O2 -mtune=native") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE") - set (OPTIMIZATION_FLAGS "-O3 -march=native -ffast-math -funroll-loops -ftree-vectorize") + set (OPTIMIZATION_FLAGS "-O3 -march=native -funroll-loops -ftree-vectorize -flto") endif () set (STANDARD_CHECK "-std=f2018 -pedantic-errors" ) -set (LINKER_FLAGS "${LINKER_FLAGS} -Wl") -# options parsed directly to the linker -set (LINKER_FLAGS "${LINKER_FLAGS},-undefined,dynamic_lookup" ) -# ensure to link against dynamic libraries #------------------------------------------------------------------------------------------------ # Fine tuning compilation options set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") # preprocessor -set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIC -fPIE") +set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIE") # position independent code set (COMPILE_FLAGS "${COMPILE_FLAGS} -ffree-line-length-132") @@ -123,6 +121,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -ffpe-trap=invalid,zero,overflow") set (DEBUG_FLAGS "${DEBUG_FLAGS} -g") # Generate symbolic debugging information in the object file +set (DEBUG_FLAGS "${DEBUG_FLAGS} -Og") +# Optimize debugging experience + set (DEBUG_FLAGS "${DEBUG_FLAGS} -fbacktrace") set (DEBUG_FLAGS "${DEBUG_FLAGS} -fdump-core") set (DEBUG_FLAGS "${DEBUG_FLAGS} -fcheck=all") diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 3afffd2be..7f34e4a13 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -9,12 +9,12 @@ if (OPENMP) set (OPENMP_FLAGS "-qopenmp -parallel") endif () -if (OPTIMIZATION STREQUAL "OFF") +if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG") set (OPTIMIZATION_FLAGS "-O0 -no-ip") elseif (OPTIMIZATION STREQUAL "DEFENSIVE") set (OPTIMIZATION_FLAGS "-O2") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE") - set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost") + set (OPTIMIZATION_FLAGS "-ipo -O3 -fp-model fast=2 -xHost") # -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost" endif () @@ -110,6 +110,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0") # generate debug information for parameters # Disabled due to ICE when compiling phase_damage.f90 (not understandable, there is no parameter in there) +set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all") +# generate complete debugging information + # Additional options # -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits # -check: Checks at runtime, where diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake index 326cfe319..883873e1c 100644 --- a/cmake/Compiler-IntelLLVM.cmake +++ b/cmake/Compiler-IntelLLVM.cmake @@ -9,7 +9,7 @@ if (OPENMP) set (OPENMP_FLAGS "-qopenmp") endif () -if (OPTIMIZATION STREQUAL "OFF") +if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG") set (OPTIMIZATION_FLAGS "-O0") elseif (OPTIMIZATION STREQUAL "DEFENSIVE") set (OPTIMIZATION_FLAGS "-O2") @@ -109,6 +109,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0") set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all") # generate debug information for parameters +set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all") +# generate complete debugging information + # Additional options # -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits # -check: Checks at runtime, where diff --git a/python/damask/VERSION b/python/damask/VERSION index b4c8c4b3c..00113e2f3 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-556-g97f849c09 +v3.0.0-alpha5-608-g3e8d1a60d diff --git a/python/damask/_grid.py b/python/damask/_grid.py index 8a9a4bbd9..1bf8ff386 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -62,8 +62,8 @@ class Grid: mat_max = np.nanmax(self.material) mat_N = self.N_materials return util.srepr([ - f'cells : {util.srepr(self.cells, " x ")}', - f'size : {util.srepr(self.size, " x ")} / m³', + f'cells: {util.srepr(self.cells, " × ")}', + f'size: {util.srepr(self.size, " × ")} / m³', f'origin: {util.srepr(self.origin," ")} / m', f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else f' (min: {mat_min}, max: {mat_max})') diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 0f7590782..c8b7e3cca 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -90,7 +90,6 @@ subroutine CPFEM_initAll call material_init(.false.) call phase_init call homogenization_init - call crystallite_init call CPFEM_init call config_deallocate diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index ed8fb611b..b24ba5480 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -68,7 +68,6 @@ subroutine CPFEM_initAll call material_init(restart=interface_restartInc>0) call phase_init call homogenization_init - call crystallite_init call CPFEM_init call config_deallocate diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index d2076c2cc..93de4053b 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -10,7 +10,8 @@ module HDF5_utilities #include use PETScSys #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) - use MPI + use MPI_f08 + use MPI, only: MPI_INFO_NULL_F90 => MPI_INFO_NULL #endif #endif @@ -162,9 +163,9 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel) character, intent(in), optional :: mode logical, intent(in), optional :: parallel - character :: m - integer(HID_T) :: plist_id - integer :: hdferr + character :: m + integer(HID_T) :: plist_id + integer :: hdferr if (present(mode)) then @@ -178,9 +179,15 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel) #ifdef PETSC if (present(parallel)) then +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) + if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr) + else + call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr) +#else if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) else call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) +#endif end if if(hdferr < 0) error stop 'HDF5 error' #endif @@ -1860,7 +1867,7 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_ globalShape !< shape of the dataset (all processes) integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer, dimension(worldsize) :: & + integer(MPI_INTEGER_KIND), dimension(worldsize) :: & readSize !< contribution of all processes integer :: hdferr integer(MPI_INTEGER_KIND) :: err_MPI @@ -1871,13 +1878,13 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_ if(hdferr < 0) error stop 'HDF5 error' !-------------------------------------------------------------------------------------------------- - readSize = 0 - readSize(worldrank+1) = int(localShape(ubound(localShape,1))) + readSize = 0_MPI_INTEGER_KIND + readSize(worldrank+1) = int(localShape(ubound(localShape,1)),MPI_INTEGER_KIND) #ifdef PETSC if (parallel) then call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) if(hdferr < 0) error stop 'HDF5 error' - call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process + call MPI_Allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' end if #endif @@ -1954,8 +1961,8 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & totalShape !< shape of the dataset (all processes) integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id - integer, dimension(worldsize) :: writeSize !< contribution of all processes - integer(HID_T) :: dcpl + integer(MPI_INTEGER_KIND), dimension(worldsize) :: writeSize !< contribution of all processes + integer(HID_T) :: dcpl integer :: hdferr integer(MPI_INTEGER_KIND) :: err_MPI integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T @@ -1974,11 +1981,11 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & !-------------------------------------------------------------------------------------------------- ! determine the global data layout among all processes - writeSize = 0 - writeSize(worldrank+1) = int(myShape(ubound(myShape,1))) + writeSize = 0_MPI_INTEGER_KIND + writeSize(worldrank+1) = int(myShape(ubound(myShape,1)),MPI_INTEGER_KIND) #ifdef PETSC if (parallel) then - call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process + call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' end if #endif @@ -2009,7 +2016,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & if (hdferr < 0) error stop 'HDF5 error' end if end if - + !-------------------------------------------------------------------------------------------------- ! create dataspace in memory (local shape) and in file (global shape) call H5Screate_simple_f(size(myShape), myShape, memspace_id, hdferr, myShape) diff --git a/src/IO.f90 b/src/IO.f90 index ec7cf9704..441a0bf3d 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -584,8 +584,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) character(len=pStringLen) :: formatString select case (warning_ID) - case (42) - msg = 'parameter has no effect' case (47) msg = 'no valid parameter for FFTW, using FFTW_PATIENT' case (207) @@ -594,8 +592,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) msg = 'crystallite responds elastically' case (601) msg = 'stiffness close to zero' - case (700) - msg = 'unknown crystal symmetry' case (709) msg = 'read only the first document' case default diff --git a/src/constants.f90 b/src/constants.f90 index 3b99763de..1effb8d95 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -9,8 +9,8 @@ module constants public real(pReal), parameter :: & - T_ROOM = 300.0_pReal, & !< Room temperature in K. ToDo: IUPAC: 298.15 - K_B = 1.38e-23_pReal, & !< Boltzmann constant in J/Kelvin - N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol + T_ROOM = 293.15_pReal, & !< Room temperature in K (20°C) + K_B = 1.380649e-23_pReal, & !< Boltzmann constant in J/Kelvin (https://doi.org/10.1351/goldbook) + N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol (https://doi.org/10.1351/goldbook) end module constants diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 77c5ef6fa..918f477b1 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -32,7 +32,7 @@ program DAMASK_grid implicit none type :: tLoadCase - type(rotation) :: rot !< rotation of BC + type(tRotation) :: rot !< rotation of BC type(tBoundaryCondition) :: stress, & !< stress BC deformation !< deformation BC (dot_F, F, or L) real(pReal) :: t, & !< length of increment diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 8de2de6ef..6e8c154ef 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -27,10 +27,10 @@ module discretization_grid private integer, dimension(3), public, protected :: & - grid !< (global) grid + cells !< (global) cells integer, public, protected :: & - grid3, & !< (local) grid in 3rd direction - grid3Offset !< (local) grid offset in 3rd direction + cells3, & !< (local) cells in 3rd direction + cells3Offset !< (local) cells offset in 3rd direction real(pReal), dimension(3), public, protected :: & geomSize !< (global) physical size real(pReal), public, protected :: & @@ -55,7 +55,7 @@ subroutine discretization_grid_init(restart) mySize, & !< domain size of this process origin !< (global) distance to origin integer, dimension(3) :: & - myGrid !< domain grid of this process + myGrid !< domain cells of this process integer, dimension(:), allocatable :: & materialAt, materialAt_global @@ -77,7 +77,7 @@ subroutine discretization_grid_init(restart) if (worldrank == 0) then fileContent = IO_read(interface_geomFile) - call readVTI(grid,geomSize,origin,materialAt_global,fileContent) + call readVTI(cells,geomSize,origin,materialAt_global,fileContent) fname = interface_geomFile if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) call results_openJobFile(parallel=.false.) @@ -88,37 +88,37 @@ subroutine discretization_grid_init(restart) end if - call MPI_Bcast(grid,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) + call MPI_Bcast(cells,3_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' - if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1') + if (cells(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1') call MPI_Bcast(geomSize,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call MPI_Bcast(origin,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - print'(/,1x,a,3(i12,1x))', 'cells a b c: ', grid - print '(1x,a,3(es12.5,1x))', 'size x y z: ', geomSize - print '(1x,a,3(es12.5,1x))', 'origin x y z: ', origin + print'(/,1x,a,i0,a,i0,a,i0)', 'cells: ', cells(1), ' × ', cells(2), ' × ', cells(3) + print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'size: ', geomSize(1), ' × ', geomSize(2), ' × ', geomSize(3), ' / m³' + print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'origin: ', origin(1), ' ', origin(2), ' ', origin(3), ' / m' - if (worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)') + if (worldsize>cells(3)) call IO_error(894, ext_msg='number of processes exceeds cells(3)') call fftw_mpi_init - devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & - int(grid(2),C_INTPTR_T), & - int(grid(1),C_INTPTR_T)/2+1, & + devNull = fftw_mpi_local_size_3d(int(cells(3),C_INTPTR_T), & + int(cells(2),C_INTPTR_T), & + int(cells(1),C_INTPTR_T)/2+1, & PETSC_COMM_WORLD, & - z, & ! domain grid size along z - z_offset) ! domain grid offset along z + z, & ! domain cells size along z + z_offset) ! domain cells offset along z if (z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes') - grid3 = int(z) - grid3Offset = int(z_offset) - size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) - size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal) - myGrid = [grid(1:2),grid3] + cells3 = int(z) + cells3Offset = int(z_offset) + size3 = geomSize(3)*real(cells3,pReal) /real(cells(3),pReal) + size3Offset = geomSize(3)*real(cells3Offset,pReal)/real(cells(3),pReal) + myGrid = [cells(1:2),cells3] mySize = [geomSize(1:2),size3] - call MPI_Gather(product(grid(1:2))*grid3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,& + 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(myGrid), 1_MPI_INTEGER_KIND,MPI_INTEGER,sendcounts,& @@ -131,10 +131,10 @@ subroutine discretization_grid_init(restart) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call discretization_init(materialAt, & - IPcoordinates0(myGrid,mySize,grid3Offset), & - Nodes0(myGrid,mySize,grid3Offset),& - merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write top layer... - (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process + IPcoordinates0(myGrid,mySize,cells3Offset), & + Nodes0(myGrid,mySize,cells3Offset),& + merge((cells(1)+1) * (cells(2)+1) * (cells3+1),& ! write top layer... + (cells(1)+1) * (cells(2)+1) * cells3,& ! ...unless not last process worldrank+1==worldsize)) !-------------------------------------------------------------------------------------------------- @@ -142,7 +142,7 @@ subroutine discretization_grid_init(restart) if (.not. restart) then call results_openJobFile call results_closeGroup(results_addGroup('geometry')) - call results_addAttribute('cells', grid, '/geometry') + call results_addAttribute('cells', cells, '/geometry') call results_addAttribute('size', geomSize,'/geometry') call results_addAttribute('origin',origin, '/geometry') call results_closeJobFile @@ -170,13 +170,13 @@ end subroutine discretization_grid_init !> @brief Parse vtk image data (.vti) !> @details https://vtk.org/Wiki/VTK_XML_Formats !-------------------------------------------------------------------------------------------------- -subroutine readVTI(grid,geomSize,origin,material, & +subroutine readVTI(cells,geomSize,origin,material, & fileContent) integer, dimension(3), intent(out) :: & - grid ! grid (across all processes!) + cells ! cells (across all processes!) real(pReal), dimension(3), intent(out) :: & - geomSize, & ! size (across all processes!) + geomSize, & ! size (across all processes!) origin ! origin (across all processes!) integer, dimension(:), intent(out), allocatable :: & material @@ -190,7 +190,7 @@ subroutine readVTI(grid,geomSize,origin,material, & s - grid = -1 + cells = -1 geomSize = -1.0_pReal inFile = .false. @@ -215,7 +215,7 @@ subroutine readVTI(grid,geomSize,origin,material, & if (.not. inImage) then if (index(fileContent(startPos:endPos),' @brief Calculate undeformed position of IPs/cell centers (pretend to be an element) !--------------------------------------------------------------------------------------------------- -function IPcoordinates0(grid,geomSize,grid3Offset) +function IPcoordinates0(cells,geomSize,cells3Offset) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) + integer, dimension(3), intent(in) :: cells ! cells (for this process!) real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) - integer, intent(in) :: grid3Offset ! grid(3) offset + integer, intent(in) :: cells3Offset ! cells(3) offset - real(pReal), dimension(3,product(grid)) :: ipCoordinates0 + real(pReal), dimension(3,product(cells)) :: ipCoordinates0 integer :: & a,b,c, & @@ -516,9 +516,9 @@ function IPcoordinates0(grid,geomSize,grid3Offset) i = 0 - do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1) + do c = 1, cells(3); do b = 1, cells(2); do a = 1, cells(1) i = i + 1 - IPcoordinates0(1:3,i) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal) + IPcoordinates0(1:3,i) = geomSize/real(cells,pReal) * (real([a,b,cells3Offset+c],pReal) -0.5_pReal) end do; end do; end do end function IPcoordinates0 @@ -527,22 +527,22 @@ end function IPcoordinates0 !--------------------------------------------------------------------------------------------------- !> @brief Calculate position of undeformed nodes (pretend to be an element) !--------------------------------------------------------------------------------------------------- -pure function nodes0(grid,geomSize,grid3Offset) +pure function nodes0(cells,geomSize,cells3Offset) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) + integer, dimension(3), intent(in) :: cells ! cells (for this process!) real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) - integer, intent(in) :: grid3Offset ! grid(3) offset + integer, intent(in) :: cells3Offset ! cells(3) offset - real(pReal), dimension(3,product(grid+1)) :: nodes0 + real(pReal), dimension(3,product(cells+1)) :: nodes0 integer :: & a,b,c, & n n = 0 - do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1) + do c = 0, cells3; do b = 0, cells(2); do a = 0, cells(1) n = n + 1 - nodes0(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal) + nodes0(1:3,n) = geomSize/real(cells,pReal) * real([a,b,cells3Offset+c],pReal) end do; end do; end do end function nodes0 @@ -551,17 +551,17 @@ end function nodes0 !-------------------------------------------------------------------------------------------------- !> @brief Calculate IP interface areas !-------------------------------------------------------------------------------------------------- -pure function cellSurfaceArea(geomSize,grid) +pure function cellSurfaceArea(geomSize,cells) real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) + integer, dimension(3), intent(in) :: cells ! cells (for this process!) - real(pReal), dimension(6,1,product(grid)) :: cellSurfaceArea + real(pReal), dimension(6,1,product(cells)) :: cellSurfaceArea - cellSurfaceArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3)) - cellSurfaceArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1)) - cellSurfaceArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2)) + cellSurfaceArea(1:2,1,:) = geomSize(2)/real(cells(2)) * geomSize(3)/real(cells(3)) + cellSurfaceArea(3:4,1,:) = geomSize(3)/real(cells(3)) * geomSize(1)/real(cells(1)) + cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1)) * geomSize(2)/real(cells(2)) end function cellSurfaceArea @@ -588,42 +588,42 @@ end function cellSurfaceNormal !-------------------------------------------------------------------------------------------------- !> @brief Build IP neighborhood relations !-------------------------------------------------------------------------------------------------- -pure function IPneighborhood(grid) +pure function IPneighborhood(cells) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) + integer, dimension(3), intent(in) :: cells ! cells (for this process!) - integer, dimension(3,6,1,product(grid)) :: IPneighborhood !< 6 neighboring IPs as [element ID, IP ID, face ID] + integer, dimension(3,6,1,product(cells)) :: IPneighborhood !< 6 neighboring IPs as [element ID, IP ID, face ID] integer :: & x,y,z, & e e = 0 - do z = 0,grid(3)-1; do y = 0,grid(2)-1; do x = 0,grid(1)-1 + do z = 0,cells(3)-1; do y = 0,cells(2)-1; do x = 0,cells(1)-1 e = e + 1 ! element ID - IPneighborhood(1,1,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x+1,grid(1)) & + IPneighborhood(1,1,1,e) = z * cells(1) * cells(2) & + + y * cells(1) & + + modulo(x+1,cells(1)) & + 1 - IPneighborhood(1,2,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x-1,grid(1)) & + IPneighborhood(1,2,1,e) = z * cells(1) * cells(2) & + + y * cells(1) & + + modulo(x-1,cells(1)) & + 1 - IPneighborhood(1,3,1,e) = z * grid(1) * grid(2) & - + modulo(y+1,grid(2)) * grid(1) & + IPneighborhood(1,3,1,e) = z * cells(1) * cells(2) & + + modulo(y+1,cells(2)) * cells(1) & + x & + 1 - IPneighborhood(1,4,1,e) = z * grid(1) * grid(2) & - + modulo(y-1,grid(2)) * grid(1) & + IPneighborhood(1,4,1,e) = z * cells(1) * cells(2) & + + modulo(y-1,cells(2)) * cells(1) & + x & + 1 - IPneighborhood(1,5,1,e) = modulo(z+1,grid(3)) * grid(1) * grid(2) & - + y * grid(1) & + IPneighborhood(1,5,1,e) = modulo(z+1,cells(3)) * cells(1) * cells(2) & + + y * cells(1) & + x & + 1 - IPneighborhood(1,6,1,e) = modulo(z-1,grid(3)) * grid(1) * grid(2) & - + y * grid(1) & + IPneighborhood(1,6,1,e) = modulo(z-1,cells(3)) * cells(1) * cells(2) & + + y * cells(1) & + x & + 1 ! IP ID diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 96b72dbcc..dc21a1a2c 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -106,9 +106,9 @@ subroutine grid_damage_spectral_init() !-------------------------------------------------------------------------------------------------- ! init fields - allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal) - allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) - allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) + allocate(phi_current(cells(1),cells(2),cells3), source=1.0_pReal) + allocate(phi_lastInc(cells(1),cells(2),cells3), source=1.0_pReal) + allocate(phi_stagInc(cells(1),cells(2),cells3), source=1.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -117,23 +117,23 @@ subroutine grid_damage_spectral_init() call SNESSetOptionsPrefix(SNES_damage,'damage_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 1_pPetscInt, 0_pPetscInt, & ! #dof (phi, scalar), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells damage_grid,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(damage_grid,err_PETSc) CHKERRQ(err_PETSc) call DMsetUp(damage_grid,err_PETSc) CHKERRQ(err_PETSc) - call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor) + call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (cells x 1, i.e. every def grad tensor) CHKERRQ(err_PETSc) call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector CHKERRQ(err_PETSc) @@ -213,7 +213,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution) !-------------------------------------------------------------------------------------------------- ! updating damage state ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 call homogenization_set_phi(phi_current(i,j,k),ce) end do; end do; end do @@ -255,7 +255,7 @@ subroutine grid_damage_spectral_forward(cutBack) call DMDAVecRestoreArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc) CHKERRQ(err_PETSc) ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 call homogenization_set_phi(phi_current(i,j,k),ce) end do; end do; end do @@ -289,12 +289,12 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) !-------------------------------------------------------------------------------------------------- ! evaluate polarization field scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current + scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of damage field call utilities_FFTvectorBackward ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k)) end do; end do; end do @@ -302,7 +302,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field call utilities_FFTscalarBackward ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & @@ -315,14 +315,14 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) call utilities_FFTscalarBackward - where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) & - scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc - where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) & - scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness + where(scalarField_real(1:cells(1),1:cells(2),1:cells3) > phi_lastInc) & + scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_lastInc + where(scalarField_real(1:cells(1),1:cells(2),1:cells3) < num%residualStiffness) & + scalarField_real(1:cells(1),1:cells(2),1:cells3) = num%residualStiffness !-------------------------------------------------------------------------------------------------- ! constructing residual - r = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current + r = scalarField_real(1:cells(1),1:cells(2),1:cells3) - phi_current err_PETSc = 0 end subroutine formResidual @@ -339,7 +339,7 @@ subroutine updateReference() K_ref = 0.0_pReal mu_ref = 0.0_pReal - do ce = 1, product(grid(1:2))*grid3 + do ce = 1, product(cells(1:2))*cells3 K_ref = K_ref + homogenization_K_phi(ce) mu_ref = mu_ref + homogenization_mu_phi(ce) end do diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 6a3f1323d..3c42c8c23 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -153,9 +153,9 @@ subroutine grid_mechanical_FEM_init !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate(F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(P_current (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(F_lastInc (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -164,16 +164,16 @@ subroutine grid_mechanical_FEM_init call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if(err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, & DMDA_STENCIL_BOX, & - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 3_pPetscInt, 1_pPetscInt, & ! #dof (u, vector), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells mechanical_grid,err_PETSc) CHKERRQ(err_PETSc) call DMsetFromOptions(mechanical_grid,err_PETSc) @@ -214,7 +214,7 @@ subroutine grid_mechanical_FEM_init call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,err_PETSc) CHKERRQ(err_PETSc) - delta = geomSize/real(grid,pReal) ! grid spacing + delta = geomSize/real(cells,pReal) ! grid spacing detJ = product(delta) ! cell volume BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & @@ -255,11 +255,11 @@ subroutine grid_mechanical_FEM_init call HDF5_read(u_lastInc,groupHandle,'u_lastInc') elseif (interface_restartInc == 0) then restartRead - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) + F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity + F = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response + homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F @@ -339,7 +339,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(rotation), intent(in) :: & + type(tRotation), intent(in) :: & rotation_BC PetscErrorCode :: err_PETSc PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -386,7 +386,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai F_lastInc = F - homogenization_F0 = reshape(F, [3,3,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F, [3,3,product(cells(1:2))*cells3]) endif !-------------------------------------------------------------------------------------------------- @@ -556,13 +556,13 @@ subroutine formResidual(da_local,x_local, & ! get deformation gradient call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc) CHKERRQ(err_PETSc) - do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ctr = 0 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0 ctr = ctr + 1 x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) enddo; enddo; enddo - F(1:3,1:3,i,j,k-grid3offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) + F(1:3,1:3,i,j,k-cells3Offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) enddo; enddo; enddo call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc) CHKERRQ(err_PETSc) @@ -589,14 +589,14 @@ subroutine formResidual(da_local,x_local, & call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc) CHKERRQ(err_PETSc) ele = 0 - do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ctr = 0 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0 ctr = ctr + 1 x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) enddo; enddo; enddo ele = ele + 1 - f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-grid3offset)))*detJ + & + f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-cells3Offset)))*detJ + & matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + & homogenization_dPdF(2,2,2,2,ele) + & homogenization_dPdF(3,3,3,3,ele))/3.0_pReal @@ -615,17 +615,17 @@ subroutine formResidual(da_local,x_local, & ! applying boundary conditions call DMDAVecGetArrayF90(da_local,f_local,r,err_PETSc) CHKERRQ(err_PETSc) - if (grid3offset == 0) then - r(0:2,0, 0, 0) = 0.0_pReal - r(0:2,grid(1),0, 0) = 0.0_pReal - r(0:2,0, grid(2),0) = 0.0_pReal - r(0:2,grid(1),grid(2),0) = 0.0_pReal + if (cells3Offset == 0) then + r(0:2,0, 0, 0) = 0.0_pReal + r(0:2,cells(1),0, 0) = 0.0_pReal + r(0:2,0, cells(2),0) = 0.0_pReal + r(0:2,cells(1),cells(2),0) = 0.0_pReal end if - if (grid3+grid3offset == grid(3)) then - r(0:2,0, 0, grid(3)) = 0.0_pReal - r(0:2,grid(1),0, grid(3)) = 0.0_pReal - r(0:2,0, grid(2),grid(3)) = 0.0_pReal - r(0:2,grid(1),grid(2),grid(3)) = 0.0_pReal + if (cells3+cells3Offset == cells(3)) then + r(0:2,0, 0, cells(3)) = 0.0_pReal + r(0:2,cells(1),0, cells(3)) = 0.0_pReal + r(0:2,0, cells(2),cells(3)) = 0.0_pReal + r(0:2,cells(1),cells(2),cells(3)) = 0.0_pReal end if call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc) CHKERRQ(err_PETSc) @@ -663,7 +663,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc) call MatZeroEntries(Jac,err_PETSc) CHKERRQ(err_PETSc) ce = 0 - do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ctr = 0 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0 ctr = ctr + 1 @@ -719,7 +719,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc) call DMDAVecGetArrayF90(da_local,coordinates,x_scal,err_PETSc) CHKERRQ(err_PETSc) ce = 0 - do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ce = ce + 1 x_scal(0:2,i-1,j-1,k-1) = discretization_IPcoords(1:3,ce) enddo; enddo; enddo diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index c60c542d0..2f2b73f01 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -79,6 +79,12 @@ module grid_mechanical_spectral_basic err_BC, & !< deviation from stress BC err_div !< RMS of div of P +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) + type(MPI_Status) :: status +#else + integer, dimension(MPI_STATUS_SIZE) :: status +#endif + integer :: & totalIter = 0 !< total iteration in current increment @@ -96,7 +102,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine grid_mechanical_spectral_basic_init - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: P PetscErrorCode :: err_PETSc integer(MPI_INTEGER_KIND) :: err_MPI PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -153,8 +159,8 @@ subroutine grid_mechanical_spectral_basic_init !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate(F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_lastInc(3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(Fdot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -163,23 +169,23 @@ subroutine grid_mechanical_spectral_basic_init call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 9_pPetscInt, 0_pPetscInt, & ! #dof (F, tensor), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells da,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(da,err_PETSc) CHKERRQ(err_PETSc) call DMsetUp(da,err_PETSc) CHKERRQ(err_PETSc) - call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor) + call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (cells x 9, i.e. every def grad tensor) CHKERRQ(err_PETSc) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector CHKERRQ(err_PETSc) @@ -217,11 +223,11 @@ subroutine grid_mechanical_spectral_basic_init call HDF5_read(F_lastInc,groupHandle,'F_lastInc') elseif (interface_restartInc == 0) then restartRead - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity + F = reshape(F_lastInc,[9,cells(1),cells(2),cells3]) end if restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response + homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -244,7 +250,7 @@ subroutine grid_mechanical_spectral_basic_init call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', & MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI) + call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call MPI_File_close(fileUnit,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' @@ -310,7 +316,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(rotation), intent(in) :: & + type(tRotation), intent(in) :: & rotation_BC PetscErrorCode :: err_PETSc PetscScalar, pointer, dimension(:,:,:,:) :: F @@ -343,11 +349,11 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ end if Fdot = utilities_calculateRate(guess, & - F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, & + F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, & rotation_BC%rotate(F_aimDot,active=.true.)) - F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) + F_lastInc = reshape(F,[3,3,cells(1),cells(2),cells3]) - homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3]) end if !-------------------------------------------------------------------------------------------------- @@ -359,7 +365,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average - rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) + rotation_BC%rotate(F_aim,active=.true.)),[9,cells(1),cells(2),cells3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc) CHKERRQ(err_PETSc) @@ -530,7 +536,7 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r ! store fPK field for subsequent FFT forward transform + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r ! store fPK field for subsequent FFT forward transform call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier @@ -538,7 +544,7 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! constructing residual - r = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too + r = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too end subroutine formResidual diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 810047403..b72cc4232 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -90,6 +90,12 @@ module grid_mechanical_spectral_polarisation err_curl, & !< RMS of curl of F err_div !< RMS of div of P +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) + type(MPI_Status) :: status +#else + integer, dimension(MPI_STATUS_SIZE) :: status +#endif + integer :: & totalIter = 0 !< total iteration in current increment @@ -107,7 +113,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine grid_mechanical_spectral_polarisation_init - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: P PetscErrorCode :: err_PETSc integer(MPI_INTEGER_KIND) :: err_MPI PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -171,10 +177,10 @@ subroutine grid_mechanical_spectral_polarisation_init !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_lastInc (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(Fdot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(F_tau_lastInc(3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(F_tauDot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -183,23 +189,23 @@ subroutine grid_mechanical_spectral_polarisation_init call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 18_pPetscInt, 0_pPetscInt, & ! #dof (2xtensor), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells da,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(da,err_PETSc) CHKERRQ(err_PETSc) call DMsetUp(da,err_PETSc) CHKERRQ(err_PETSc) - call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor) + call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (cells x 18, i.e. every def grad tensor) CHKERRQ(err_PETSc) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector CHKERRQ(err_PETSc) @@ -241,13 +247,13 @@ subroutine grid_mechanical_spectral_polarisation_init call HDF5_read(F_tau_lastInc,groupHandle,'F_tau_lastInc') elseif (interface_restartInc == 0) then restartRead - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity + F = reshape(F_lastInc,[9,cells(1),cells(2),cells3]) F_tau = 2.0_pReal*F F_tau_lastInc = 2.0_pReal*F_lastInc end if restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response + homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -270,7 +276,7 @@ subroutine grid_mechanical_spectral_polarisation_init call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', & MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI) + call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call MPI_File_close(fileUnit,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' @@ -342,7 +348,7 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(rotation), intent(in) :: & + type(tRotation), intent(in) :: & rotation_BC PetscErrorCode :: err_PETSc PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau @@ -379,15 +385,15 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D end if Fdot = utilities_calculateRate(guess, & - F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, & + F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, & rotation_BC%rotate(F_aimDot,active=.true.)) F_tauDot = utilities_calculateRate(guess, & - F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), Delta_t_old, & + F_tau_lastInc,reshape(F_tau,[3,3,cells(1),cells(2),cells3]), Delta_t_old, & rotation_BC%rotate(F_aimDot,active=.true.)) - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) - F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) + F_lastInc = reshape(F, [3,3,cells(1),cells(2),cells3]) + F_tau_lastInc = reshape(F_tau,[3,3,cells(1),cells(2),cells3]) - homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3]) end if !-------------------------------------------------------------------------------------------------- @@ -400,12 +406,12 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),& - [9,grid(1),grid(2),grid3]) + [9,cells(1),cells(2),cells3]) if (guess) then F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), & - [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + [9,cells(1),cells(2),cells3]) ! does not have any average value as boundary condition else - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) F_lambda33 = math_I3 & + math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, & @@ -597,7 +603,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! tensorField_real = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) tensorField_real(1:3,1:3,i,j,k) = & num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& num%alpha*matmul(F(1:3,1:3,i,j,k), & @@ -612,7 +618,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! constructing residual - r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -629,14 +635,14 @@ subroutine formResidual(in, FandF_tau, & params%stress_mask))) ! calculate divergence tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r_F !< stress field in disguise + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r_F !< stress field in disguise call utilities_FFTtensorForward err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual e = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) e = e + 1 r_F(1:3,1:3,i,j,k) = & math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & @@ -648,7 +654,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! calculating curl tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F call utilities_FFTtensorForward err_curl = utilities_curlRMS() diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 00aa61379..1066a1f61 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -105,9 +105,9 @@ subroutine grid_thermal_spectral_init(T_0) !-------------------------------------------------------------------------------------------------- ! init fields - allocate(T_current(grid(1),grid(2),grid3), source=T_0) - allocate(T_lastInc(grid(1),grid(2),grid3), source=T_0) - allocate(T_stagInc(grid(1),grid(2),grid3), source=T_0) + 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) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -116,23 +116,23 @@ subroutine grid_thermal_spectral_init(T_0) call SNESSetOptionsPrefix(SNES_thermal,'thermal_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 1_pPetscInt, 0_pPetscInt, & ! #dof (T, scalar), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells thermal_grid,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(thermal_grid,err_PETSc) CHKERRQ(err_PETSc) call DMsetUp(thermal_grid,err_PETSc) CHKERRQ(err_PETSc) - call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor) + call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (cells x 1, i.e. every def grad tensor) CHKERRQ(err_PETSc) call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector CHKERRQ(err_PETSc) @@ -150,10 +150,10 @@ subroutine grid_thermal_spectral_init(T_0) call HDF5_read(T_current,groupHandle,'T',.false.) call HDF5_read(T_lastInc,groupHandle,'T_lastInc',.false.) - end if restartRead + end if restartRead ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) ce = ce + 1 call homogenization_thermal_setField(T_current(i,j,k),0.0_pReal,ce) end do; end do; end do @@ -164,7 +164,7 @@ subroutine grid_thermal_spectral_init(T_0) call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc) CHKERRQ(err_PETSc) - call updateReference + call updateReference() end subroutine grid_thermal_spectral_init @@ -214,7 +214,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) !-------------------------------------------------------------------------------------------------- ! updating thermal state ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) end do; end do; end do @@ -257,7 +257,7 @@ subroutine grid_thermal_spectral_forward(cutBack) call DMDAVecRestoreArrayF90(dm_local,solution_vec,T_PETSc,err_PETSc) CHKERRQ(err_PETSc) ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) end do; end do; end do @@ -279,9 +279,9 @@ subroutine grid_thermal_spectral_restartWrite integer(HID_T) :: fileHandle, groupHandle PetscScalar, dimension(:,:,:), pointer :: T - call SNESGetDM(SNES_thermal,dm_local,err_PETSc); + call SNESGetDM(SNES_thermal,dm_local,err_PETSc); CHKERRQ(err_PETSc) - call DMDAVecGetArrayF90(dm_local,solution_vec,T,err_PETSc); + call DMDAVecGetArrayF90(dm_local,solution_vec,T,err_PETSc); CHKERRQ(err_PETSc) print'(1x,a)', 'writing thermal solver data required for restart to file'; flush(IO_STDOUT) @@ -293,7 +293,7 @@ subroutine grid_thermal_spectral_restartWrite call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) - call DMDAVecRestoreArrayF90(dm_local,solution_vec,T,err_PETSc); + call DMDAVecRestoreArrayF90(dm_local,solution_vec,T,err_PETSc); CHKERRQ(err_PETSc) end subroutine grid_thermal_spectral_restartWrite @@ -321,12 +321,12 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) !-------------------------------------------------------------------------------------------------- ! evaluate polarization field scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = T_current + scalarField_real(1:cells(1),1:cells(2),1:cells3) = T_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of temperature field call utilities_FFTvectorBackward ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k)) end do; end do; end do @@ -334,7 +334,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field call utilities_FFTscalarBackward ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & @@ -349,7 +349,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) !-------------------------------------------------------------------------------------------------- ! constructing residual - r = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) + r = T_current - scalarField_real(1:cells(1),1:cells(2),1:cells3) err_PETSc = 0 end subroutine formResidual @@ -366,7 +366,7 @@ subroutine updateReference() K_ref = 0.0_pReal mu_ref = 0.0_pReal - do ce = 1, product(grid(1:2))*grid3 + do ce = 1, product(cells(1:2))*cells3 K_ref = K_ref + homogenization_K_T(ce) mu_ref = mu_ref + homogenization_mu_T(ce) end do diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index ee5fd4c82..d95580145 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -29,9 +29,9 @@ module spectral_utilities include 'fftw3-mpi.f03' !-------------------------------------------------------------------------------------------------- -! grid related information information +! grid related information real(pReal), protected, public :: wgt !< weighting factor 1/Nelems - integer, protected, public :: grid1Red !< grid(1)/2 + integer, protected, public :: grid1Red !< cells(1)/2 real(pReal), protected, public, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence !-------------------------------------------------------------------------------------------------- @@ -86,7 +86,7 @@ module spectral_utilities type, public :: tSolutionParams real(pReal), dimension(3,3) :: stress_BC logical, dimension(3,3) :: stress_mask - type(rotation) :: rotation_BC + type(tRotation) :: rotation_BC real(pReal) :: Delta_t end type tSolutionParams @@ -201,8 +201,8 @@ subroutine spectral_utilities_init num_grid%get_asString('PETSc_options',defaultVal=''),err_PETSc) CHKERRQ(err_PETSc) - grid1Red = grid(1)/2 + 1 - wgt = 1.0/real(product(grid),pReal) + grid1Red = cells(1)/2 + 1 + wgt = 1.0/real(product(cells),pReal) num%memory_efficient = num_grid%get_asInt('memory_efficient', defaultVal=1) > 0 ! ToDo: should be logical in YAML file num%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2) @@ -231,9 +231,9 @@ subroutine spectral_utilities_init enddo elseif (num%divergence_correction == 2) then do j = 1, 3 - if ( j /= int(minloc(geomSize/real(grid,pReal),1)) & - .and. j /= int(maxloc(geomSize/real(grid,pReal),1))) & - scaledGeomSize = geomSize/geomSize(j)*real(grid(j),pReal) + if ( j /= int(minloc(geomSize/real(cells,pReal),1)) & + .and. j /= int(maxloc(geomSize/real(cells,pReal),1))) & + scaledGeomSize = geomSize/geomSize(j)*real(cells(j),pReal) enddo else scaledGeomSize = geomSize @@ -262,11 +262,11 @@ subroutine spectral_utilities_init !-------------------------------------------------------------------------------------------------- ! MPI allocation - gridFFTW = int(grid,C_INTPTR_T) + gridFFTW = int(cells,C_INTPTR_T) alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, & PETSC_COMM_WORLD, local_K, local_K_offset) - allocate (xi1st (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension - allocate (xi2nd (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension + allocate (xi1st (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension + allocate (xi2nd (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension tensorField = fftw_alloc_complex(tensorSize*alloc_local) call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, & @@ -327,27 +327,27 @@ subroutine spectral_utilities_init !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - do k = grid3Offset+1, grid3Offset+grid3 + do k = cells3Offset+1, cells3Offset+cells3 k_s(3) = k - 1 - if (k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 - do j = 1, grid(2) + if (k > cells(3)/2 + 1) k_s(3) = k_s(3) - cells(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + do j = 1, cells(2) k_s(2) = j - 1 - if (j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + if (j > cells(2)/2 + 1) k_s(2) = k_s(2) - cells(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 do i = 1, grid1Red k_s(1) = i - 1 ! symmetry, junst running from 0,1,...,N/2,N/2+1 - xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s) - where(mod(grid,2)==0 .and. [i,j,k] == grid/2+1 .and. & + xi2nd(1:3,i,j,k-cells3Offset) = utilities_getFreqDerivative(k_s) + where(mod(cells,2)==0 .and. [i,j,k] == cells/2+1 .and. & spectral_derivative_ID == DERIVATIVE_CONTINUOUS_ID) ! for even grids, set the Nyquist Freq component to 0.0 - xi1st(1:3,i,j,k-grid3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal) + xi1st(1:3,i,j,k-cells3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal) elsewhere - xi1st(1:3,i,j,k-grid3Offset) = xi2nd(1:3,i,j,k-grid3Offset) + xi1st(1:3,i,j,k-cells3Offset) = xi2nd(1:3,i,j,k-cells3Offset) endwhere enddo; enddo; enddo - if (num%memory_efficient) then ! allocate just single fourth order tensor + if (num%memory_efficient) then ! allocate just single fourth order tensor allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal)) else ! precalculation of gamma_hat field - allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal)) + allocate (gamma_hat(3,3,3,3,grid1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal)) endif end subroutine spectral_utilities_init @@ -373,10 +373,10 @@ subroutine utilities_updateGamma(C) if (.not. num%memory_efficient) then gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A - do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, grid1Red if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 do concurrent (l = 1:3, m = 1:3) - xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset) + xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) end do do concurrent(l = 1:3, m = 1:3) temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) @@ -387,8 +387,8 @@ subroutine utilities_updateGamma(C) call math_invert(A_inv, err, A) temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) do concurrent(l=1:3, m=1:3, n=1:3, o=1:3) - gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* & - conjg(-xi1st(o,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset) + gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_complex(l,n)* & + conjg(-xi1st(o,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) end do end if end if @@ -405,7 +405,7 @@ end subroutine utilities_updateGamma !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTtensorForward - tensorField_real(1:3,1:3,grid(1)+1:grid1Red*2,:,:) = 0.0_pReal + tensorField_real(1:3,1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) end subroutine utilities_FFTtensorForward @@ -429,7 +429,7 @@ end subroutine utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTscalarForward - scalarField_real(grid(1)+1:grid1Red*2,:,:) = 0.0_pReal + scalarField_real(cells(1)+1:grid1Red*2,:,:) = 0.0_pReal call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) end subroutine utilities_FFTscalarForward @@ -454,7 +454,7 @@ end subroutine utilities_FFTscalarBackward !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTvectorForward - vectorField_real(1:3,grid(1)+1:grid1Red*2,:,:) = 0.0_pReal + vectorField_real(1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) end subroutine utilities_FFTvectorForward @@ -493,8 +493,8 @@ subroutine utilities_fourierGammaConvolution(fieldAim) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation (mechanical equilibrium) memoryEfficient: if (num%memory_efficient) then - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red - if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red + if (any([i,j,k+cells3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 do concurrent(l = 1:3, m = 1:3) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) end do @@ -519,7 +519,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim) end if end do; end do; end do else memoryEfficient - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red do concurrent(l = 1:3, m = 1:3) temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k)) end do @@ -527,7 +527,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim) end do; end do; end do end if memoryEfficient - if (grid3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) + if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) end subroutine utilities_fourierGammaConvolution @@ -544,7 +544,7 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation - do k = 1, grid3; do j = 1, grid(2) ;do i = 1, grid1Red + do k = 1, cells3; do j = 1, cells(2) ;do i = 1, grid1Red GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) & / (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal) & * sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k)))) @@ -571,7 +571,7 @@ real(pReal) function utilities_divergenceRMS() !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space utilities_divergenceRMS = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2) + do k = 1, cells3; do j = 1, cells(2) do i = 2, grid1Red -1 ! Has somewhere a conj. complex counterpart. Therefore count it twice. utilities_divergenceRMS = utilities_divergenceRMS & + 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k), & ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2, i.e. do not take square root and square again @@ -579,7 +579,7 @@ real(pReal) function utilities_divergenceRMS() +sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),& conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2)) enddo - utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1) + utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if cells(1) /= 1) + sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), & conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) & + sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), & @@ -589,7 +589,7 @@ real(pReal) function utilities_divergenceRMS() + sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), & conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2) enddo; enddo - if (grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + if (cells(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of cells(1) == 1 call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space @@ -616,7 +616,7 @@ real(pReal) function utilities_curlRMS() ! calculating max curl criterion in Fourier space utilities_curlRMS = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); + do k = 1, cells3; do j = 1, cells(2); do i = 2, grid1Red - 1 do l = 1, 3 curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) & @@ -638,7 +638,7 @@ real(pReal) function utilities_curlRMS() -tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2)) enddo utilities_curlRMS = utilities_curlRMS & - + sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1) + + sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (DC) does not have a conjugate complex counterpart (if cells(1) /= 1) do l = 1, 3 curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) & -tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3)) @@ -648,13 +648,13 @@ real(pReal) function utilities_curlRMS() -tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2)) enddo utilities_curlRMS = utilities_curlRMS & - + sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1) + + sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (Nyquist) does not have a conjugate complex counterpart (if cells(1) /= 1) enddo; enddo call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' utilities_curlRMS = sqrt(utilities_curlRMS) * wgt - if (grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + if (cells(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of cells(1) == 1 end function utilities_curlRMS @@ -666,7 +666,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance real(pReal), intent(in), dimension(3,3,3,3) :: C !< current average stiffness - type(rotation), intent(in) :: rot_BC !< rotation of load frame + type(tRotation), intent(in) :: rot_BC !< rotation of load frame logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC integer :: i, j @@ -736,7 +736,7 @@ subroutine utilities_fourierScalarGradient() integer :: i, j, k - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) ! ToDo: no -conjg? enddo; enddo; enddo @@ -750,7 +750,7 @@ subroutine utilities_fourierVectorDivergence() integer :: i, j, k - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red scalarField_fourier(i,j,k) = sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k))) enddo; enddo; enddo @@ -764,7 +764,7 @@ subroutine utilities_fourierVectorGradient() integer :: i, j, k, m, n - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red do m = 1, 3; do n = 1, 3 tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k) enddo; enddo @@ -780,7 +780,7 @@ subroutine utilities_fourierTensorDivergence() integer :: i, j, k - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red vectorField_fourier(:,i,j,k) = matmul(tensorField_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k))) enddo; enddo; enddo @@ -795,10 +795,10 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress - real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target + real(pReal), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress + real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: F !< deformation gradient target real(pReal), intent(in) :: Delta_t !< loading time - type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame + type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame integer :: i @@ -810,15 +810,15 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& print'(/,1x,a)', '... evaluating constitutive response ......................................' flush(IO_STDOUT) - homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field - call homogenization_mechanical_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field + call homogenization_mechanical_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) ! calculate P field if (.not. terminallyIll) & - call homogenization_thermal_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) + call homogenization_thermal_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) if (.not. terminallyIll) & - call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) + call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) - P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) + P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' @@ -833,7 +833,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& dPdF_norm_max = 0.0_pReal dPdF_min = huge(1.0_pReal) dPdF_norm_min = huge(1.0_pReal) - do i = 1, product(grid(1:2))*grid3 + do i = 1, product(cells(1:2))*cells3 if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i) dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) @@ -878,16 +878,16 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) dt !< Delta_t between field0 and field logical, intent(in) :: & heterogeneous !< calculate field of rates - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & + real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: & field0, & !< data of previous step field !< data of current step - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: & utilities_calculateRate if (heterogeneous) then utilities_calculateRate = (field-field0) / dt else - utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3) + utilities_calculateRate = spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3) endif end function utilities_calculateRate @@ -901,12 +901,12 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim) real(pReal), intent(in) :: & Delta_t !< Delta_t of current step - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & + real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: & field_lastInc, & !< initial field rate !< rate by which to forward real(pReal), intent(in), optional, dimension(3,3) :: & aim !< average field value aim - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: & utilities_forwardField real(pReal), dimension(3,3) :: fieldDiff !< - aim integer(MPI_INTEGER_KIND) :: err_MPI @@ -918,7 +918,7 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' fieldDiff = fieldDiff - aim utilities_forwardField = utilities_forwardField - & - spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid3) + spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3) endif end function utilities_forwardField @@ -936,37 +936,37 @@ pure function utilities_getFreqDerivative(k_s) select case (spectral_derivative_ID) case (DERIVATIVE_CONTINUOUS_ID) - utilities_getFreqDerivative = cmplx(0.0_pReal, 2.0_pReal*PI*real(k_s,pReal)/geomSize,pReal) + utilities_getFreqDerivative = cmplx(0.0_pReal, TAU*real(k_s,pReal)/geomSize,pReal) case (DERIVATIVE_CENTRAL_DIFF_ID) - utilities_getFreqDerivative = cmplx(0.0_pReal, sin(2.0_pReal*PI*real(k_s,pReal)/real(grid,pReal)), pReal)/ & - cmplx(2.0_pReal*geomSize/real(grid,pReal), 0.0_pReal, pReal) + utilities_getFreqDerivative = cmplx(0.0_pReal, sin(TAU*real(k_s,pReal)/real(cells,pReal)), pReal)/ & + cmplx(2.0_pReal*geomSize/real(cells,pReal), 0.0_pReal, pReal) case (DERIVATIVE_FWBW_DIFF_ID) utilities_getFreqDerivative(1) = & - cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) - 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(1)/real(grid(1),pReal), 0.0_pReal, pReal) + cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) - 1.0_pReal, & + sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ & + cmplx(4.0_pReal*geomSize(1)/real(cells(1),pReal), 0.0_pReal, pReal) utilities_getFreqDerivative(2) = & - cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) - 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(2)/real(grid(2),pReal), 0.0_pReal, pReal) + cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) - 1.0_pReal, & + sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ & + cmplx(4.0_pReal*geomSize(2)/real(cells(2),pReal), 0.0_pReal, pReal) utilities_getFreqDerivative(3) = & - cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) - 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(3)/real(grid(3),pReal), 0.0_pReal, pReal) + cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) - 1.0_pReal, & + sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ & + cmplx(4.0_pReal*geomSize(3)/real(cells(3),pReal), 0.0_pReal, pReal) end select end function utilities_getFreqDerivative @@ -979,10 +979,10 @@ end function utilities_getFreqDerivative !-------------------------------------------------------------------------------------------------- subroutine utilities_updateCoords(F) - real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F - real(pReal), dimension(3, grid(1),grid(2),grid3) :: IPcoords - real(pReal), dimension(3, grid(1),grid(2),grid3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI) - real(pReal), dimension(3, grid(1)+1,grid(2)+1,grid3+1) :: nodeCoords + real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: F + real(pReal), dimension(3, cells(1),cells(2),cells3) :: IPcoords + real(pReal), dimension(3, cells(1),cells(2),cells3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI) + real(pReal), dimension(3, cells(1)+1,cells(2)+1,cells3+1) :: nodeCoords integer :: & i,j,k,n, & c @@ -1010,14 +1010,14 @@ subroutine utilities_updateCoords(F) 1, 1, 1, & 0, 1, 1 ], [3,8]) - step = geomSize/real(grid, pReal) + step = geomSize/real(cells, pReal) !-------------------------------------------------------------------------------------------------- ! integration in Fourier space to get fluctuations of cell center discplacements - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F call utilities_FFTtensorForward() - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red - if (any([i,j,k+grid3Offset] /= 1)) then + do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red + if (any([i,j,k+cells3Offset] /= 1)) then vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) & / sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal) else @@ -1029,25 +1029,25 @@ subroutine utilities_updateCoords(F) !-------------------------------------------------------------------------------------------------- ! average F - if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt + if (cells3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' !-------------------------------------------------------------------------------------------------- ! pad cell center fluctuations along z-direction (needed when running MPI simulation) - IPfluct_padded(1:3,1:grid(1),1:grid(2),2:grid3+1) = vectorField_real(1:3,1:grid(1),1:grid(2),1:grid3) - c = product(shape(IPfluct_padded(:,:,:,1))) !< amount of data to transfer + IPfluct_padded(1:3,1:cells(1),1:cells(2),2:cells3+1) = vectorField_real(1:3,1:cells(1),1:cells(2),1:cells3) + c = product(shape(IPfluct_padded(:,:,:,1))) !< amount of data to transfer rank_t = modulo(worldrank+1_MPI_INTEGER_KIND,worldsize) rank_b = modulo(worldrank-1_MPI_INTEGER_KIND,worldsize) ! send bottom layer to process below call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(1),err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI) + call MPI_Irecv(IPfluct_padded(:,:,:,cells3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' ! send top layer to process above - call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI) + call MPI_Isend(IPfluct_padded(:,:,:,cells3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(4),err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' @@ -1063,24 +1063,24 @@ subroutine utilities_updateCoords(F) !-------------------------------------------------------------------------------------------------- ! calculate nodal displacements nodeCoords = 0.0_pReal - do k = 0,grid3; do j = 0,grid(2); do i = 0,grid(1) - nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+grid3Offset],pReal))) + do k = 0,cells3; do j = 0,cells(2); do i = 0,cells(1) + nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+cells3Offset],pReal))) averageFluct: do n = 1,8 me = [i+neighbor(1,n),j+neighbor(2,n),k+neighbor(3,n)] nodeCoords(1:3,i+1,j+1,k+1) = nodeCoords(1:3,i+1,j+1,k+1) & - + IPfluct_padded(1:3,modulo(me(1)-1,grid(1))+1,modulo(me(2)-1,grid(2))+1,me(3)+1)*0.125_pReal + + IPfluct_padded(1:3,modulo(me(1)-1,cells(1))+1,modulo(me(2)-1,cells(2))+1,me(3)+1)*0.125_pReal enddo averageFluct enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! calculate cell center displacements - do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1) + do k = 1,cells3; do j = 1,cells(2); do i = 1,cells(1) IPcoords(1:3,i,j,k) = vectorField_real(1:3,i,j,k) & - + matmul(Favg,step*(real([i,j,k+grid3Offset],pReal)-0.5_pReal)) + + matmul(Favg,step*(real([i,j,k+cells3Offset],pReal)-0.5_pReal)) enddo; enddo; enddo - call discretization_setNodeCoords(reshape(NodeCoords,[3,(grid(1)+1)*(grid(2)+1)*(grid3+1)])) - call discretization_setIPcoords (reshape(IPcoords, [3,grid(1)*grid(2)*grid3])) + call discretization_setNodeCoords(reshape(NodeCoords,[3,(cells(1)+1)*(cells(2)+1)*(cells3+1)])) + call discretization_setIPcoords (reshape(IPcoords, [3,cells(1)*cells(2)*cells3])) end subroutine utilities_updateCoords diff --git a/src/lattice.f90 b/src/lattice.f90 index e18c71edf..2c2840bca 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -498,7 +498,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA) real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin real(pReal), dimension(3,3,sum(Ntwin)):: coordinateSystem - type(rotation) :: R + type(tRotation) :: R integer :: i @@ -537,7 +537,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66 real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S - type(rotation) :: R + type(tRotation) :: R real(pReal) :: a_bcc, a_fcc, cOverA_trans integer :: i @@ -599,7 +599,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem !< coordinate system of slip system real(pReal), dimension(3) :: direction, normal, np - type(rotation) :: R + type(tRotation) :: R integer :: i @@ -1976,7 +1976,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) a_bcc, & !< lattice parameter a for bcc target lattice a_fcc !< lattice parameter a for fcc parent lattice - type(rotation) :: & + type(tRotation) :: & R, & !< Pitsch rotation B !< Rotation of fcc to Bain coordinate system real(pReal), dimension(3,3) :: & diff --git a/src/material.f90 b/src/material.f90 index 31eeef5e1..fd1531964 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -18,7 +18,7 @@ module material private type :: tRotationContainer - type(Rotation), dimension(:), allocatable :: data + type(tRotation), dimension(:), allocatable :: data end type type :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data diff --git a/src/math.f90 b/src/math.f90 index db0666ab0..dd4690672 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -21,10 +21,11 @@ module math config #endif - real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter - real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian to degree - real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree to radian - complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) + real(pReal), parameter :: & + PI = acos(-1.0_pReal), & !< ratio of a circle's circumference to its diameter + TAU = 2.0_pReal*PI, & !< ratio of a circle's circumference to its radius + INDEG = 360.0_pReal/TAU, & !< conversion from radian to degree + INRAD = TAU/360.0_pReal !< conversion from degree to radian real(pReal), dimension(3,3), parameter :: & math_I3 = reshape([& @@ -882,7 +883,7 @@ end function math_Voigt6to33_strain !-------------------------------------------------------------------------------------------------- -!> @brief Convert 3x3 tensor into 6 Voigt stress vector. +!> @brief Convert 3x3 stress tensor into 6 Voigt vector. !-------------------------------------------------------------------------------------------------- pure function math_33toVoigt6_stress(sigma) result(sigma_tilde) @@ -897,7 +898,7 @@ end function math_33toVoigt6_stress !-------------------------------------------------------------------------------------------------- -!> @brief Convert 3x3 tensor into 6 Voigt strain vector. +!> @brief Convert 3x3 strain tensor into 6 Voigt vector. !-------------------------------------------------------------------------------------------------- pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde) @@ -913,48 +914,48 @@ end function math_33toVoigt6_strain !-------------------------------------------------------------------------------------------------- -!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix. +!> @brief Convert 6x6 Voigt stiffness matrix into symmetric 3x3x3x3 tensor. !-------------------------------------------------------------------------------------------------- -pure function math_Voigt66to3333(m66) +pure function math_Voigt66to3333_stiffness(C_tilde) result(C) - real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 - real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix + real(pReal), dimension(3,3,3,3) :: C + real(pReal), dimension(6,6), intent(in) :: C_tilde integer :: i,j do i=1,6; do j=1,6 - math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) - math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) - math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j) - math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j) + C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = C_tilde(i,j) + C(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = C_tilde(i,j) + C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = C_tilde(i,j) + C(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = C_tilde(i,j) end do; end do -end function math_Voigt66to3333 +end function math_Voigt66to3333_stiffness !-------------------------------------------------------------------------------------------------- -!> @brief Convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix. +!> @brief Convert 3x3x3x3 stiffness tensor into 6x6 Voigt matrix. !-------------------------------------------------------------------------------------------------- -pure function math_3333toVoigt66(m3333) +pure function math_3333toVoigt66_stiffness(C) result(C_tilde) - real(pReal), dimension(6,6) :: math_3333toVoigt66 - real(pReal), dimension(3,3,3,3), intent(in) :: m3333 !< symmetric 3x3x3x3 matrix (no internal check) + real(pReal), dimension(6,6) :: C_tilde + real(pReal), dimension(3,3,3,3), intent(in) :: C integer :: i,j #ifndef __INTEL_COMPILER do concurrent(i=1:6, j=1:6) - math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) + C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) end do #else do i=1,6; do j=1,6 - math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) + C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) end do; end do #endif -end function math_3333toVoigt66 +end function math_3333toVoigt66_stiffness !-------------------------------------------------------------------------------------------------- @@ -984,7 +985,7 @@ impure elemental subroutine math_normal(x,mu,sigma) end if call random_number(rnd) - x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(2.0_pReal*PI*(1.0_pReal - rnd(2))) + x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(TAU*(1.0_pReal - rnd(2))) end subroutine math_normal @@ -1088,7 +1089,7 @@ pure function math_rotationalPart(F) result(R) if (dNeq0(x)) then Phi = acos(math_clip((I_C(1)**3 -4.5_pReal*I_C(1)*I_C(2) +13.5_pReal*I_C(3))/x,-1.0_pReal,1.0_pReal)) lambda = I_C(1) +(2.0_pReal * sqrt(math_clip(I_C(1)**2-3.0_pReal*I_C(2),0.0_pReal))) & - *cos((Phi-2.0_pReal * PI*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal) + *cos((Phi-TAU*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal) lambda = sqrt(math_clip(lambda,0.0_pReal)/3.0_pReal) else lambda = sqrt(I_C(1)/3.0_pReal) @@ -1154,8 +1155,8 @@ pure function math_eigvalsh33(m) phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) math_eigvalsh33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & [cos( phi /3.0_pReal), & - cos((phi+2.0_pReal*PI)/3.0_pReal), & - cos((phi+4.0_pReal*PI)/3.0_pReal) & + cos((phi+TAU)/3.0_pReal), & + cos((phi+2.0_pReal*TAU)/3.0_pReal) & ] & + I(1)/3.0_pReal endif @@ -1343,7 +1344,7 @@ subroutine selfTest if (any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & error stop 'math_sym3333to66/math_66toSym3333' - if (any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) & + if (any(dNeq(math_3333toVoigt66_stiffness(math_Voigt66to3333_stiffness(t66)),t66,1.0e-15_pReal))) & error stop 'math_3333toVoigt66/math_Voigt66to3333' call random_number(v6) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index aa156ec2d..065eca1cb 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -307,9 +307,11 @@ program DAMASK_mesh guess = .true. ! start guessing after first converged (sub)inc timeIncOld = timeinc end if - if (.not. cutBack .and. worldrank == 0) & + if (.not. cutBack .and. worldrank == 0) then write(statUnit,*) totalIncsCounter, time, cutBackLevel, & solres%converged, solres%iterationsNeeded ! write statistics about accepted solution + flush(statUnit) + endif end do subStepLooping cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc diff --git a/src/parallelization.f90 b/src/parallelization.f90 index 29deaf724..28ad70d94 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -52,13 +52,13 @@ contains !-------------------------------------------------------------------------------------------------- subroutine parallelization_init - integer(MPI_INTEGER_KIND) :: err_MPI, typeSize + integer(MPI_INTEGER_KIND) :: err_MPI, typeSize, version, subversion, devNull character(len=4) :: rank_str + character(len=MPI_MAX_LIBRARY_VERSION_STRING) :: MPI_library_version !$ integer :: got_env, threadLevel !$ integer(pI32) :: OMP_NUM_THREADS !$ character(len=6) NumThreadsString - PetscErrorCode :: err_PETSc #ifdef _OPENMP ! If openMP is enabled, check if the MPI libary supports it and initialize accordingly. @@ -86,12 +86,22 @@ subroutine parallelization_init if (err_MPI /= 0_MPI_INTEGER_KIND) & error stop 'Could not determine worldrank' - if (worldrank == 0) print'(/,1x,a)', '<<<+- parallelization init -+>>>' + if (worldrank == 0) then + print'(/,1x,a)', '<<<+- parallelization init -+>>>' + + call MPI_Get_library_version(MPI_library_version,devNull,err_MPI) + print'(/,1x,a)', trim(MPI_library_version) + call MPI_Get_version(version,subversion,err_MPI) + print'(1x,a,i0,a,i0)', 'MPI standard: ',version,'.',subversion +#ifdef _OPENMP + print'(1x,a,i0)', 'OpenMP version: ',openmp_version +#endif + end if call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) & error stop 'Could not determine worldsize' - if (worldrank == 0) print'(/,1x,a,i3)', 'MPI processes: ',worldsize + if (worldrank == 0) print'(/,1x,a,i0)', 'MPI processes: ',worldsize call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) & @@ -128,7 +138,7 @@ subroutine parallelization_init !$ OMP_NUM_THREADS = 4_pI32 !$ endif !$ endif -!$ print'(1x,a,1x,i2)', 'OMP_NUM_THREADS:',OMP_NUM_THREADS +!$ print'(1x,a,i0)', 'OMP_NUM_THREADS: ',OMP_NUM_THREADS !$ call omp_set_num_threads(OMP_NUM_THREADS) end subroutine parallelization_init diff --git a/src/phase.f90 b/src/phase.f90 index bea1a5543..3de551131 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -323,7 +323,6 @@ module phase phase_restore, & plastic_nonlocal_updateCompatibility, & converged, & - crystallite_init, & phase_mechanical_constitutive, & phase_thermal_constitutive, & phase_damage_constitutive, & @@ -401,6 +400,8 @@ subroutine phase_init call damage_init call thermal_init(phases) + call crystallite_init() + end subroutine phase_init @@ -408,7 +409,7 @@ end subroutine phase_init !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- subroutine phase_allocateState(state, & - NEntries,sizeState,sizeDotState,sizeDeltaState) + NEntries,sizeState,sizeDotState,sizeDeltaState,offsetDeltaState) class(tState), intent(inout) :: & state @@ -417,12 +418,17 @@ subroutine phase_allocateState(state, & sizeState, & sizeDotState, & sizeDeltaState - + integer, intent(in), optional :: & + offsetDeltaState state%sizeState = sizeState state%sizeDotState = sizeDotState state%sizeDeltaState = sizeDeltaState - state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition + if (present(offsetDeltaState)) then + state%offsetDeltaState = offsetDeltaState ! ToDo: this is a fix for broken nonlocal + else + state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition + end if allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%state0 (sizeState,NEntries), source=0.0_pReal) @@ -431,7 +437,8 @@ subroutine phase_allocateState(state, & allocate(state%dotState (sizeDotState,NEntries), source=0.0_pReal) allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal) - + state%deltaState2 => state%state(state%offsetDeltaState+1: & + state%offsetDeltaState+state%sizeDeltaState,:) end subroutine phase_allocateState @@ -496,22 +503,13 @@ subroutine crystallite_init() ce, & co, & !< counter in integration point component loop ip, & !< counter in integration point loop - el, & !< counter in element loop - cMax, & !< maximum number of integration point components - iMax, & !< maximum number of integration points - eMax !< maximum number of elements + el !< counter in element loop class(tNode), pointer :: & num_crystallite, & phases - print'(/,1x,a)', '<<<+- crystallite init -+>>>' - - cMax = homogenization_maxNconstituents - iMax = discretization_nIPs - eMax = discretization_Nelems - num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal) @@ -545,15 +543,9 @@ subroutine crystallite_init() phases => config_material%get('phase') - print'(/,a42,1x,i10)', ' # of elements: ', eMax - print'( a42,1x,i10)', ' # of integration points/element: ', iMax - print'( a42,1x,i10)', 'max # of constituents/integration point: ', cMax - flush(IO_STDOUT) - - !$OMP PARALLEL DO PRIVATE(ce) - do el = 1, eMax - do ip = 1, iMax + do el = 1, discretization_Nelems + do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) call crystallite_orientations(co,ip,el) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 2252b941d..eb758d8ed 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -40,8 +40,6 @@ submodule(phase) mechanical integer(kind(PLASTIC_undefined_ID)), dimension(:), allocatable :: & phase_plasticity !< plasticity of each phase - integer :: phase_plasticity_maxSizeDotState - interface module subroutine eigen_init(phases) @@ -81,16 +79,17 @@ submodule(phase) mechanical en end subroutine plastic_isotropic_LiAndItsTangent - module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) + module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element + co, & !< constituent + ip, & !< integration point + el, & !< element ph, & en real(pReal), intent(in) :: & - subdt !< timestep - logical :: broken + subdt !< timestep + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState end function plastic_dotState module function plastic_deltaState(ph, en) result(broken) @@ -296,8 +295,6 @@ module subroutine mechanical_init(phases) do ph = 1,phases%length plasticState(ph)%state0 = plasticState(ph)%state enddo - phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) - num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) @@ -588,9 +585,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop logical :: & broken @@ -601,43 +598,43 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul sizeDotState real(pReal) :: & zeta - real(pReal), dimension(phase_plasticity_maxSizeDotState) :: & - r ! state residuum - real(pReal), dimension(phase_plasticity_maxSizeDotState,2) :: & + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + r, & ! state residuum dotState + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: & + dotState_last ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + broken = .true. - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState - plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t + plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState * Delta_t iteration: do NiterationState = 1, num%nState - dotState(1:sizeDotState,2) = merge(dotState(1:sizeDotState,1),0.0, nIterationState > 1) - dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,en) + dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1) + dotState_last(1:sizeDotState,1) = dotState broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) exit iteration - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) exit iteration + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) exit iteration - zeta = damper(plasticState(ph)%dotState(:,en),dotState(1:sizeDotState,1),& - dotState(1:sizeDotState,2)) - plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) * zeta & - + dotState(1:sizeDotState,1) * (1.0_pReal - zeta) - r(1:sizeDotState) = plasticState(ph)%state(1:sizeDotState,en) & - - subState0 & - - plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t - plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) & - - r(1:sizeDotState) - if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then + zeta = damper(dotState,dotState_last(1:sizeDotState,1),dotState_last(1:sizeDotState,2)) + dotState = dotState * zeta & + + dotState_last(1:sizeDotState,1) * (1.0_pReal - zeta) + r = plasticState(ph)%state(1:sizeDotState,en) & + - subState0 & + - dotState * Delta_t + plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) - r + + if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then broken = plastic_deltaState(ph,en) exit iteration endif @@ -652,19 +649,20 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul !-------------------------------------------------------------------------------------------------- real(pReal) pure function damper(omega_0,omega_1,omega_2) - real(pReal), dimension(:), intent(in) :: & - omega_0, omega_1, omega_2 + real(pReal), dimension(:), intent(in) :: & + omega_0, omega_1, omega_2 - real(pReal) :: dot_prod12, dot_prod22 + real(pReal) :: dot_prod12, dot_prod22 - dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2) - dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2) - if (min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif + dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2) + dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2) + + if (min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif end function damper @@ -686,6 +684,8 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res logical :: & broken + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + dotState integer :: & ph, & en, & @@ -694,13 +694,14 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + broken = .true. - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t + + dotState * Delta_t broken = plastic_deltaState(ph,en) if(broken) return @@ -729,20 +730,23 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip ph, & en, & sizeDotState - real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + r, & + dotState ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + broken = .true. - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState - residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,en) * 0.5_pReal * Delta_t + r = - dotState * 0.5_pReal * Delta_t plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotstate(1:sizeDotState,en) * Delta_t + + dotState * Delta_t broken = plastic_deltaState(ph,en) if(broken) return @@ -750,10 +754,10 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) return - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return - broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,en) * Delta_t, & + broken = .not. converged(r + 0.5_pReal * dotState * Delta_t, & plasticState(ph)%state(1:sizeDotState,en), & plasticState(ph)%atol(1:sizeDotState)) @@ -847,44 +851,48 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D ph, & en, & sizeDotState - real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + dotState + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: & + plastic_RKdotState ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + broken = .true. - broken = plastic_dotState(Delta_t,co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState do stage = 1, size(A,1) - plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,en) - plasticState(ph)%dotState(:,en) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) + plastic_RKdotState(1:sizeDotState,stage) = dotState + dotState = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) do n = 2, stage - plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) & - + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) + dotState = dotState & + + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) enddo plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t + + dotState * Delta_t - broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el) + broken = integrateStress(F_0 + (F-F_0) * Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage),co,ip,el) if(broken) exit - broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,en) - if(broken) exit + dotState = plastic_dotState(Delta_t*C(stage), co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) exit enddo if(broken) return - plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,en) - plasticState(ph)%dotState(:,en) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) + plastic_RKdotState(1:sizeDotState,size(B)) = dotState + dotState = matmul(plastic_RKdotState,B) plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t + + dotState * Delta_t if(present(DB)) & broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & @@ -958,7 +966,7 @@ subroutine crystallite_results(group,ph) !-------------------------------------------------------------------------------------------------- function to_quaternion(dataset) - type(rotation), dimension(:), intent(in) :: dataset + type(tRotation), dimension(:), intent(in) :: dataset real(pReal), dimension(4,size(dataset,1)) :: to_quaternion integer :: i diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index a5d9868a8..3c422616b 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -27,7 +27,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics - integer :: Ninstances,p,i,k + integer :: Ninstances, p, k class(tNode), pointer :: & phases, & phase, & diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 084cc75bc..7f30c8165 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -48,7 +48,7 @@ module subroutine elastic_init(phases) prm%C_11 = polynomial(elastic%asDict(),'C_11','T') prm%C_12 = polynomial(elastic%asDict(),'C_12','T') prm%C_44 = polynomial(elastic%asDict(),'C_44','T') - + if (any(phase_lattice(ph) == ['hP','tI'])) then prm%C_13 = polynomial(elastic%asDict(),'C_13','T') prm%C_33 = polynomial(elastic%asDict(),'C_33','T') @@ -162,7 +162,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & C66 = phase_damage_C66(phase_homogenizedC66(ph,en),ph,en) - C = math_Voigt66to3333(C66) + C = math_Voigt66to3333_stiffness(C66) E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration S = math_Voigt6to33_stress(matmul(C66,math_33toVoigt6_strain(matmul(matmul(transpose(Fi),E),Fi))))!< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 0a24b0cbf..70181d724 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -110,29 +110,35 @@ submodule(phase:mechanical) plastic end subroutine nonlocal_LpAndItsTangent - module subroutine isotropic_dotState(Mp,ph,en) + module function isotropic_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en - end subroutine isotropic_dotState + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState + end function isotropic_dotState - module subroutine phenopowerlaw_dotState(Mp,ph,en) + module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en - end subroutine phenopowerlaw_dotState + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState + end function phenopowerlaw_dotState - module subroutine plastic_kinehardening_dotState(Mp,ph,en) + module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en - end subroutine plastic_kinehardening_dotState + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState + end function plastic_kinehardening_dotState module subroutine dislotwin_dotState(Mp,T,ph,en) real(pReal), dimension(3,3), intent(in) :: & @@ -144,21 +150,20 @@ submodule(phase:mechanical) plastic en end subroutine dislotwin_dotState - module subroutine dislotungsten_dotState(Mp,T,ph,en) + module function dislotungsten_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & - T integer, intent(in) :: & ph, & en - end subroutine dislotungsten_dotState + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState + end function dislotungsten_dotState - module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,en,ip,el) + module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & - Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & @@ -285,7 +290,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & do i=1,3; do j=1,3 dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) - dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) + dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) end do; end do end if @@ -296,7 +301,7 @@ end subroutine plastic_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) +module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) integer, intent(in) :: & co, & !< component-ID of integration point @@ -308,7 +313,8 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) subdt !< timestep real(pReal), dimension(3,3) :: & Mp - logical :: broken + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then @@ -318,27 +324,28 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) plasticType: select case (phase_plasticity(ph)) case (PLASTIC_ISOTROPIC_ID) plasticType - call isotropic_dotState(Mp,ph,en) + dotState = isotropic_dotState(Mp,ph,en) case (PLASTIC_PHENOPOWERLAW_ID) plasticType - call phenopowerlaw_dotState(Mp,ph,en) + dotState = phenopowerlaw_dotState(Mp,ph,en) case (PLASTIC_KINEHARDENING_ID) plasticType - call plastic_kinehardening_dotState(Mp,ph,en) + dotState = plastic_kinehardening_dotState(Mp,ph,en) case (PLASTIC_DISLOTWIN_ID) plasticType call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en) + dotState = plasticState(ph)%dotState(:,en) case (PLASTIC_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en) + dotState = dislotungsten_dotState(Mp,ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el) + call nonlocal_dotState(Mp,subdt,ph,en,ip,el) + dotState = plasticState(ph)%dotState(:,en) + end select plasticType end if - broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,en))) - end function plastic_dotState @@ -393,6 +400,7 @@ module function plastic_deltaState(ph, en) result(broken) myOffset, & mySize + broken = .false. select case (phase_plasticity(ph)) @@ -414,10 +422,9 @@ module function plastic_deltaState(ph, en) result(broken) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) if (.not. broken) then - myOffset = plasticState(ph)%offsetDeltaState - mySize = plasticState(ph)%sizeDeltaState - plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = & - plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en) + mySize = plasticState(ph)%sizeDeltaState + plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) & + + plasticState(ph)%deltaState(1:mySize,en) end if end select diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index cd71a7fd7..5a6ac8f5c 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -43,6 +43,13 @@ submodule(phase:plastic) dislotungsten systems_sl end type tParameters !< container type for internal constitutive parameters + type :: tIndexDotState + integer, dimension(2) :: & + rho_mob, & + rho_dip, & + gamma_sl + end type tIndexDotState + type :: tDislotungstenState real(pReal), dimension(:,:), pointer :: & rho_mob, & @@ -58,10 +65,9 @@ submodule(phase:plastic) dislotungsten !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:) :: param - type(tDisloTungstenState), allocatable, dimension(:) :: & - dotState, & - state + type(tParameters), allocatable, dimension(:) :: param + type(tIndexDotState), allocatable, dimension(:) :: indexDotState + type(tDisloTungstenState), allocatable, dimension(:) :: state type(tDisloTungstenDependentState), allocatable, dimension(:) :: dependentState contains @@ -103,18 +109,17 @@ module function plastic_dislotungsten_init() result(myPlasticity) print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002' - phases => config_material%get('phase') allocate(param(phases%length)) + allocate(indexDotState(phases%length)) allocate(state(phases%length)) - allocate(dotState(phases%length)) allocate(dependentState(phases%length)) - do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), & + idx_dot => indexDotState(ph)) phase => phases%get(ph) mech => phase%get('mechanical') @@ -214,28 +219,29 @@ module function plastic_dislotungsten_init() result(myPlasticity) sizeState = sizeDotState call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl + idx_dot%rho_mob = [startIndex,endIndex] stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:) stt%rho_mob = spread(rho_mob_0,2,Nmembers) - dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl + idx_dot%rho_dip = [startIndex,endIndex] stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:) stt%rho_dip = spread(rho_dip_0,2,Nmembers) - dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl + idx_dot%gamma_sl = [startIndex,endIndex] stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:) - dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' @@ -300,15 +306,15 @@ end subroutine dislotungsten_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine dislotungsten_dotState(Mp,T,ph,en) +module function dislotungsten_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature integer, intent(in) :: & ph, & en + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos, dot_gamma_neg,& @@ -319,17 +325,22 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_rho_dip_climb, & d_hat real(pReal) :: & - mu + mu, T - associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) + + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), & + dot_rho_mob => dotState(indexDotState(ph)%rho_mob(1):indexDotState(ph)%rho_mob(2)), & + dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), & + dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2))) mu = elastic_mu(ph,en) + T = thermal_T(ph,en) call kinetics(Mp,T,ph,en,& dot_gamma_pos,dot_gamma_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg) + dot_gamma_sl = abs(dot_gamma_pos+dot_gamma_neg) where(dEq0((tau_pos+tau_neg)*0.5_pReal)) dot_rho_dip_formation = 0.0_pReal @@ -338,24 +349,24 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) d_hat = math_clip(3.0_pReal*mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), & prm%d_caron, & ! lower limit dst%Lambda_sl(:,en)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & + dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot_gamma_sl/prm%b_sl, & 0.0_pReal, & prm%dipoleformation) - v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(2.0_pReal*PI*K_B*T)) & + v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(TAU*K_B*T)) & * (1.0_pReal/(d_hat+prm%d_caron)) dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication + dot_rho_mob = dot_gamma_sl/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication - dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges - dot%rho_dip(:,en) = dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot_gamma_sl ! Spontaneous annihilation of 2 edges + dot_rho_dip = dot_rho_dip_formation & + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot_gamma_sl & ! Spontaneous annihilation of an edge with a dipole - dot_rho_dip_climb end associate -end subroutine dislotungsten_dotState +end function dislotungsten_dotState !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index bea5339c7..49efc9fae 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -37,9 +37,7 @@ submodule(phase:plastic) isotropic !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param - type(tIsotropicState), allocatable, dimension(:) :: & - dotState, & - state + type(tIsotropicState), allocatable, dimension(:) :: state contains @@ -77,16 +75,15 @@ module function plastic_isotropic_init() result(myPlasticity) phases => config_material%get('phase') allocate(param(phases%length)) allocate(state(phases%length)) - allocate(dotState(phases%length)) do ph = 1, phases%length if(.not. myPlasticity(ph)) cycle - associate(prm => param(ph), dot => dotState(ph), stt => state(ph)) + associate(prm => param(ph), stt => state(ph)) phase => phases%get(ph) - mech => phase%get('mechanical') - pl => mech%get('plastic') + mech => phase%get('mechanical') + pl => mech%get('plastic') #if defined (__GFORTRAN__) prm%output = output_as1dString(pl) @@ -125,12 +122,12 @@ module function plastic_isotropic_init() result(myPlasticity) sizeState = sizeDotState call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- ! state aliases and initialization - stt%xi => plasticState(ph)%state (1,:) - stt%xi = xi_0 - dot%xi => plasticState(ph)%dotState(1,:) + stt%xi => plasticState(ph)%state(1,:) + stt%xi = xi_0 plasticState(ph)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if (plasticState(ph)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi' @@ -178,7 +175,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) norm_Mp_dev = sqrt(squarenorm_Mp_dev) if (norm_Mp_dev > 0.0_pReal) then - dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en)))**prm%n Lp = dot_gamma * Mp_dev/norm_Mp_dev forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -242,27 +239,26 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine isotropic_dotState(Mp,ph,en) +module function isotropic_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState real(pReal) :: & dot_gamma, & !< strainrate xi_inf_star, & !< saturation xi norm_Mp !< norm of the (deviatoric) Mandel stress - associate(prm => param(ph), stt => state(ph), & - dot => dotState(ph)) + associate(prm => param(ph), stt => state(ph), dot_xi => dotState(1)) - if (prm%dilatation) then - norm_Mp = sqrt(math_tensordot(Mp,Mp)) - else - norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))) - end if + norm_Mp = merge(sqrt(math_tensordot(Mp,Mp)), & + sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))), & + prm%dilatation) dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n @@ -274,16 +270,16 @@ module subroutine isotropic_dotState(Mp,ph,en) + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) & / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n) end if - dot%xi(en) = dot_gamma & - * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) & - * sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star) + dot_xi = dot_gamma & + * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) & + * sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star) else - dot%xi(en) = 0.0_pReal + dot_xi = 0.0_pReal end if end associate -end subroutine isotropic_dotState +end function isotropic_dotState !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 03eb27f31..0f6bd53d4 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -34,6 +34,13 @@ submodule(phase:plastic) kinehardening systems_sl end type tParameters + type :: tIndexDotState + integer, dimension(2) :: & + xi, & + chi, & + gamma + end type tIndexDotState + type :: tKinehardeningState real(pReal), pointer, dimension(:,:) :: & xi, & !< resistance against plastic slip @@ -47,10 +54,8 @@ submodule(phase:plastic) kinehardening !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param - type(tKinehardeningState), allocatable, dimension(:) :: & - dotState, & - deltaState, & - state + type(tIndexDotState), allocatable, dimension(:) :: indexDotState + type(tKinehardeningState), allocatable, dimension(:) :: state, deltaState contains @@ -91,19 +96,20 @@ module function plastic_kinehardening_init() result(myPlasticity) phases => config_material%get('phase') allocate(param(phases%length)) + allocate(indexDotState(phases%length)) allocate(state(phases%length)) - allocate(dotState(phases%length)) allocate(deltaState(phases%length)) do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph)) + associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph), & + idx_dot => indexDotState(ph)) phase => phases%get(ph) - mech => phase%get('mechanical') - pl => mech%get('plastic') + mech => phase%get('mechanical') + pl => mech%get('plastic') #if defined (__GFORTRAN__) prm%output = output_as1dString(pl) @@ -173,27 +179,28 @@ module function plastic_kinehardening_init() result(myPlasticity) sizeState = sizeDotState + sizeDeltaState call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%xi => plasticState(ph)%state (startIndex:endIndex,:) + idx_dot%xi = [startIndex,endIndex] + stt%xi => plasticState(ph)%state(startIndex:endIndex,:) stt%xi = spread(xi_0, 2, Nmembers) - dot%xi => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%chi => plasticState(ph)%state (startIndex:endIndex,:) - dot%chi => plasticState(ph)%dotState(startIndex:endIndex,:) + idx_dot%chi = [startIndex,endIndex] + stt%chi => plasticState(ph)%state(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma => plasticState(ph)%dotState(startIndex:endIndex,:) + idx_dot%gamma = [startIndex,endIndex] + stt%gamma => plasticState(ph)%state(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' @@ -270,13 +277,15 @@ end subroutine kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_kinehardening_dotState(Mp,ph,en) +module function plastic_kinehardening_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState real(pReal) :: & sumGamma @@ -284,29 +293,32 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en) dot_gamma_pos,dot_gamma_neg - associate(prm => param(ph), stt => state(ph),dot => dotState(ph)) + associate(prm => param(ph), stt => state(ph), & + dot_xi => dotState(IndexDotState(ph)%xi(1):IndexDotState(ph)%xi(2)),& + dot_chi => dotState(IndexDotState(ph)%chi(1):IndexDotState(ph)%chi(2)),& + dot_gamma => dotState(IndexDotState(ph)%gamma(1):IndexDotState(ph)%gamma(2))) call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg) - dot%gamma(:,en) = abs(dot_gamma_pos+dot_gamma_neg) + dot_gamma = abs(dot_gamma_pos+dot_gamma_neg) sumGamma = sum(stt%gamma(:,en)) - dot%xi(:,en) = matmul(prm%h_sl_sl,dot%gamma(:,en)) & + dot_xi = matmul(prm%h_sl_sl,dot_gamma) & * ( prm%h_inf_f & + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) & * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) & ) - dot%chi(:,en) = stt%sgn_gamma(:,en)*dot%gamma(:,en) * & - ( prm%h_inf_b + & - (prm%h_0_b - prm%h_inf_b & + dot_chi = stt%sgn_gamma(:,en)*dot_gamma & + * ( prm%h_inf_b & + + (prm%h_0_b - prm%h_inf_b & + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi_0(:,en))*(stt%gamma(:,en)-stt%gamma_0(:,en))& ) *exp(-(stt%gamma(:,en)-stt%gamma_0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi_0(:,en))) & ) end associate -end subroutine plastic_kinehardening_dotState +end function plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 1715af2d9..d668e40fe 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -227,8 +227,8 @@ module function plastic_nonlocal_init() result(myPlasticity) st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph)) phase => phases%get(ph) - mech => phase%get('mechanical') - pl => mech%get('plastic') + mech => phase%get('mechanical') + pl => mech%get('plastic') plasticState(ph)%nonlocal = pl%get_asBool('flux',defaultVal=.True.) #if defined (__GFORTRAN__) @@ -403,7 +403,7 @@ module function plastic_nonlocal_init() result(myPlasticity) 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention allocate(geom(ph)%V_0(Nmembers)) call storeGeometry(ph) @@ -411,9 +411,6 @@ module function plastic_nonlocal_init() result(myPlasticity) if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & call IO_error(212,ext_msg='IPneighborhood does not exist') - - plasticState(ph)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention - st0%rho => plasticState(ph)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) stt%rho => plasticState(ph)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) dot%rho => plasticState(ph)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) @@ -941,13 +938,12 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -module subroutine nonlocal_dotState(Mp, Temperature,timestep, & +module subroutine nonlocal_dotState(Mp,timestep, & ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & - Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & @@ -984,7 +980,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & real(pReal) :: & D_SD, & mu, & - nu + nu, Temperature if (timestep <= 0.0_pReal) then plasticState(ph)%dotState = 0.0_pReal @@ -995,6 +991,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) + Temperature = thermal_T(ph,en) tau = 0.0_pReal dot_gamma = 0.0_pReal @@ -1195,7 +1192,6 @@ function rhoDotFlux(timestep,ph,en,ip,el) associate(prm => param(ph), & dst => dependentState(ph), & - dot => dotState(ph), & stt => state(ph)) ns = prm%sum_N_sl @@ -1206,7 +1202,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) rho0 = getRho0(ph,en) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) @@ -1217,7 +1213,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) if (plasticState(ph)%nonlocal) then !*** check CFL (Courant-Friedrichs-Lewy) condition for flux - if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... + if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... .and. prm%C_CFL * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG @@ -1395,7 +1391,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) nThresholdValues logical, dimension(param(ph)%sum_N_sl) :: & belowThreshold - type(rotation) :: mis + type(tRotation) :: mis associate(prm => param(ph)) diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 6e79968f3..0b018562c 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -47,6 +47,14 @@ submodule(phase:plastic) phenopowerlaw systems_tw end type tParameters + type :: tIndexDotState + integer, dimension(2) :: & + xi_sl, & + xi_tw, & + gamma_sl, & + gamma_tw + end type tIndexDotState + type :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & xi_sl, & @@ -56,11 +64,10 @@ submodule(phase:plastic) phenopowerlaw end type tPhenopowerlawState !-------------------------------------------------------------------------------------------------- -! containers for parameters and state +! containers for parameters, dot state index, and state type(tParameters), allocatable, dimension(:) :: param - type(tPhenopowerlawState), allocatable, dimension(:) :: & - dotState, & - state + type(tIndexDotState), allocatable, dimension(:) :: indexDotState + type(tPhenopowerlawState), allocatable, dimension(:) :: state contains @@ -101,17 +108,18 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) phases => config_material%get('phase') allocate(param(phases%length)) + allocate(indexDotState(phases%length)) allocate(state(phases%length)) - allocate(dotState(phases%length)) do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), dot => dotState(ph), stt => state(ph)) + associate(prm => param(ph), stt => state(ph), & + idx_dot => indexDotState(ph)) phase => phases%get(ph) - mech => phase%get('mechanical') - pl => mech%get('plastic') + mech => phase%get('mechanical') + pl => mech%get('plastic') !-------------------------------------------------------------------------------------------------- ! slip related parameters @@ -224,37 +232,37 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw sizeState = sizeDotState - call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) + deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%xi_sl => plasticState(ph)%state (startIndex:endIndex,:) + idx_dot%xi_sl = [startIndex,endIndex] + stt%xi_sl => plasticState(ph)%state(startIndex:endIndex,:) stt%xi_sl = spread(xi_0_sl, 2, Nmembers) - dot%xi_sl => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%xi_tw => plasticState(ph)%state (startIndex:endIndex,:) + idx_dot%xi_tw = [startIndex,endIndex] + stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:) stt%xi_tw = spread(xi_0_tw, 2, Nmembers) - dot%xi_tw => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma_sl => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) + idx_dot%gamma_sl = [startIndex,endIndex] + stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%gamma_tw => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:) + idx_dot%gamma_tw = [startIndex,endIndex] + stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) end associate @@ -324,13 +332,15 @@ end subroutine phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine phenopowerlaw_dotState(Mp,ph,en) +module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & en + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & + dotState real(pReal) :: & xi_sl_sat_offset,& @@ -340,28 +350,32 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en) right_SlipSlip - associate(prm => param(ph), stt => state(ph), dot => dotState(ph)) + associate(prm => param(ph), stt => state(ph), & + dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), & + dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), & + dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), & + dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2))) call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg) - dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg) - call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en)) + dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg) + call kinetics_tw(Mp,ph,en,dot_gamma_tw) sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char) xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF) right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, & 1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset)) - dot%xi_sl(:,en) = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) & - * matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) & - + matmul(prm%h_sl_tw,dot%gamma_tw(:,en)) + dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) & + * matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) & + + matmul(prm%h_sl_tw,dot_gamma_tw) - dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & - * matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) & - + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en)) + dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & + * matmul(prm%h_tw_sl,dot_gamma_sl) & + + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw) end associate -end subroutine phenopowerlaw_dotState +end function phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- diff --git a/src/prec.f90 b/src/prec.f90 index 8de82fee8..61fa141ba 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -45,6 +45,8 @@ module prec state, & !< state dotState, & !< rate of state change deltaState !< increment of state change + real(pReal), pointer, dimension(:,:) :: & + deltaState2 end type type, extends(tState) :: tPlasticState diff --git a/src/rotations.f90 b/src/rotations.f90 index b73fbe8da..53f1fe976 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -55,7 +55,7 @@ module rotations real(pReal), parameter :: P = -1.0_pReal !< parameter for orientation conversion. - type, public :: rotation + type, public :: tRotation real(pReal), dimension(4) :: q contains procedure, public :: asQuaternion @@ -78,10 +78,9 @@ module rotations procedure, public :: rotStiffness procedure, public :: misorientation procedure, public :: standardize - end type rotation + end type tRotation real(pReal), parameter :: & - SPI = sqrt(PI), & PREF = sqrt(6.0_pReal/PI), & A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), & AP = PI**(2.0_pReal/3.0_pReal), & @@ -118,8 +117,8 @@ end subroutine rotations_init !--------------------------------------------------------------------------------------------------- pure function asQuaternion(self) - class(rotation), intent(in) :: self - real(pReal), dimension(4) :: asQuaternion + class(tRotation), intent(in) :: self + real(pReal), dimension(4) :: asQuaternion asQuaternion = self%q @@ -127,8 +126,8 @@ end function asQuaternion !--------------------------------------------------------------------------------------------------- pure function asEulers(self) - class(rotation), intent(in) :: self - real(pReal), dimension(3) :: asEulers + class(tRotation), intent(in) :: self + real(pReal), dimension(3) :: asEulers asEulers = qu2eu(self%q) @@ -136,8 +135,8 @@ end function asEulers !--------------------------------------------------------------------------------------------------- pure function asAxisAngle(self) - class(rotation), intent(in) :: self - real(pReal), dimension(4) :: asAxisAngle + class(tRotation), intent(in) :: self + real(pReal), dimension(4) :: asAxisAngle asAxisAngle = qu2ax(self%q) @@ -145,8 +144,8 @@ end function asAxisAngle !--------------------------------------------------------------------------------------------------- pure function asMatrix(self) - class(rotation), intent(in) :: self - real(pReal), dimension(3,3) :: asMatrix + class(tRotation), intent(in) :: self + real(pReal), dimension(3,3) :: asMatrix asMatrix = qu2om(self%q) @@ -154,8 +153,8 @@ end function asMatrix !--------------------------------------------------------------------------------------------------- pure function asRodrigues(self) - class(rotation), intent(in) :: self - real(pReal), dimension(4) :: asRodrigues + class(tRotation), intent(in) :: self + real(pReal), dimension(4) :: asRodrigues asRodrigues = qu2ro(self%q) @@ -163,8 +162,8 @@ end function asRodrigues !--------------------------------------------------------------------------------------------------- pure function asHomochoric(self) - class(rotation), intent(in) :: self - real(pReal), dimension(3) :: asHomochoric + class(tRotation), intent(in) :: self + real(pReal), dimension(3) :: asHomochoric asHomochoric = qu2ho(self%q) @@ -175,7 +174,7 @@ end function asHomochoric !--------------------------------------------------------------------------------------------------- subroutine fromQuaternion(self,qu) - class(rotation), intent(out) :: self + class(tRotation), intent(out) :: self real(pReal), dimension(4), intent(in) :: qu if (dNeq(norm2(qu),1.0_pReal,1.0e-8_pReal)) call IO_error(402,ext_msg='fromQuaternion') @@ -186,7 +185,7 @@ end subroutine fromQuaternion !--------------------------------------------------------------------------------------------------- subroutine fromEulers(self,eu,degrees) - class(rotation), intent(out) :: self + class(tRotation), intent(out) :: self real(pReal), dimension(3), intent(in) :: eu logical, intent(in), optional :: degrees @@ -198,7 +197,7 @@ subroutine fromEulers(self,eu,degrees) Eulers = merge(eu*INRAD,eu,degrees) endif - if (any(Eulers<0.0_pReal) .or. any(Eulers>2.0_pReal*PI) .or. Eulers(2) > PI) & + if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) & call IO_error(402,ext_msg='fromEulers') self%q = eu2qu(Eulers) @@ -207,7 +206,7 @@ end subroutine fromEulers !--------------------------------------------------------------------------------------------------- subroutine fromAxisAngle(self,ax,degrees,P) - class(rotation), intent(out) :: self + class(tRotation), intent(out) :: self real(pReal), dimension(4), intent(in) :: ax logical, intent(in), optional :: degrees integer, intent(in), optional :: P @@ -237,7 +236,7 @@ end subroutine fromAxisAngle !--------------------------------------------------------------------------------------------------- subroutine fromMatrix(self,om) - class(rotation), intent(out) :: self + class(tRotation), intent(out) :: self real(pReal), dimension(3,3), intent(in) :: om if (dNeq(math_det33(om),1.0_pReal,tol=1.0e-5_pReal)) & @@ -254,10 +253,10 @@ end subroutine fromMatrix !--------------------------------------------------------------------------------------------------- pure elemental function rotRot__(self,R) result(rRot) - type(rotation) :: rRot - class(rotation), intent(in) :: self,R + type(tRotation) :: rRot + class(tRotation), intent(in) :: self,R - rRot = rotation(multiply_quaternion(self%q,R%q)) + rRot = tRotation(multiply_quaternion(self%q,R%q)) call rRot%standardize() end function rotRot__ @@ -268,7 +267,7 @@ end function rotRot__ !--------------------------------------------------------------------------------------------------- pure elemental subroutine standardize(self) - class(rotation), intent(inout) :: self + class(tRotation), intent(inout) :: self if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q @@ -282,7 +281,7 @@ end subroutine standardize pure function rotVector(self,v,active) result(vRot) real(pReal), dimension(3) :: vRot - class(rotation), intent(in) :: self + class(tRotation), intent(in) :: self real(pReal), intent(in), dimension(3) :: v logical, intent(in), optional :: active @@ -318,7 +317,7 @@ end function rotVector pure function rotTensor2(self,T,active) result(tRot) real(pReal), dimension(3,3) :: tRot - class(rotation), intent(in) :: self + class(tRotation), intent(in) :: self real(pReal), intent(in), dimension(3,3) :: T logical, intent(in), optional :: active @@ -347,7 +346,7 @@ end function rotTensor2 pure function rotTensor4(self,T,active) result(tRot) real(pReal), dimension(3,3,3,3) :: tRot - class(rotation), intent(in) :: self + class(tRotation), intent(in) :: self real(pReal), intent(in), dimension(3,3,3,3) :: T logical, intent(in), optional :: active @@ -379,7 +378,7 @@ end function rotTensor4 pure function rotStiffness(self,C,active) result(cRot) real(pReal), dimension(6,6) :: cRot - class(rotation), intent(in) :: self + class(tRotation), intent(in) :: self real(pReal), intent(in), dimension(6,6) :: C logical, intent(in), optional :: active @@ -416,8 +415,8 @@ end function rotStiffness !--------------------------------------------------------------------------------------------------- pure elemental function misorientation(self,other) - type(rotation) :: misorientation - class(rotation), intent(in) :: self, other + type(tRotation) :: misorientation + class(tRotation), intent(in) :: self, other misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q)) @@ -480,7 +479,7 @@ pure function qu2eu(qu) result(eu) atan2( 2.0_pReal*chi, q03-q12 ), & atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )] endif degenerated - where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) + where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU]) end function qu2eu @@ -628,7 +627,7 @@ pure function om2eu(om) result(eu) eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ] end if where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal - where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) + where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU]) end function om2eu @@ -1209,7 +1208,7 @@ pure function ho2cu(ho) result(cu) else special q2 = qxy + maxval(abs(xyz2))**2 sq2 = sqrt(q2) - q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2)) + q = (BETA/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2)) tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], & [ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], & @@ -1217,7 +1216,7 @@ pure function ho2cu(ho) result(cu) endif special ! inverse M_1 - xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc + xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / PREF ]/SC ! reverse the coordinates back to order according to the original pyramid number cu = xyz1(p(:,2)) @@ -1323,32 +1322,32 @@ pure function cu2ho(cu) result(ho) else center ! get pyramide and scale by grid parameter ratio p = GetPyramidOrder(cu) - XYZ = cu(p(:,1)) * sc + XYZ = cu(p(:,1)) * SC ! intercept all the points along the z-axis special: if (all(dEq0(XYZ(1:2)))) then - LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ] + LamXYZ = [ 0.0_pReal, 0.0_pReal, PREF * XYZ(3) ] else special order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger c = cos(q) s = sin(q) - q = prek * XYZ(order(2))/ sqrt(R2-c) + q = PREK * XYZ(order(2))/ sqrt(R2-c) T = [ (R2*c - 1.0), R2 * s] * q ! transform to sphere grid (inverse Lambert) ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero] c = sum(T**2) - s = Pi * c/(24.0*XYZ(3)**2) - c = sPi * c / sqrt(24.0_pReal) / XYZ(3) + s = PI * c/(24.0*XYZ(3)**2) + c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3) q = sqrt( 1.0 - s ) - LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ] - endif special + LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ] + end if special ! reverse the coordinates back to order according to the original pyramid number ho = LamXYZ(p(:,2)) - endif center + end if center end function cu2ho @@ -1416,7 +1415,7 @@ end function conjugate_quaternion !-------------------------------------------------------------------------------------------------- subroutine selfTest() - type(rotation) :: R + type(tRotation) :: R real(pReal), dimension(4) :: qu, ax, ro real(pReal), dimension(3) :: x, eu, ho, v3 real(pReal), dimension(3,3) :: om, t33 @@ -1437,7 +1436,7 @@ subroutine selfTest() elseif(i==2) then qu = eu2qu([0.0_pReal,0.0_pReal,0.0_pReal]) elseif(i==3) then - qu = eu2qu([2.0_pReal*PI,PI,2.0_pReal*PI]) + qu = eu2qu([TAU,PI,TAU]) elseif(i==4) then qu = [0.0_pReal,0.0_pReal,1.0_pReal,0.0_pReal] elseif(i==5) then @@ -1448,10 +1447,10 @@ subroutine selfTest() call random_number(x) A = sqrt(x(3)) B = sqrt(1-0_pReal -x(3)) - qu = [cos(2.0_pReal*PI*x(1))*A,& - sin(2.0_pReal*PI*x(2))*B,& - cos(2.0_pReal*PI*x(2))*B,& - sin(2.0_pReal*PI*x(1))*A] + qu = [cos(TAU*x(1))*A,& + sin(TAU*x(2))*B,& + cos(TAU*x(2))*B,& + sin(TAU*x(1))*A] if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal) endif @@ -1504,7 +1503,8 @@ subroutine selfTest() call random_number(C) C = C+transpose(C) - if (any(dNeq(R%rotStiffness(C),math_3333toVoigt66(R%rotate(math_Voigt66to3333(C))),1.0e-12_pReal))) & + if (any(dNeq(R%rotStiffness(C), & + math_3333toVoigt66_stiffness(R%rotate(math_Voigt66to3333_stiffness(C))),1.0e-12_pReal))) & error stop 'rotStiffness' call R%fromQuaternion(qu * (1.0_pReal + merge(+5.e-9_pReal,-5.e-9_pReal, mod(i,2) == 0))) ! allow reasonable tolerance for ASCII/YAML diff --git a/src/system_routines.f90 b/src/system_routines.f90 index 7a4e41a57..2eb0b7958 100644 --- a/src/system_routines.f90 +++ b/src/system_routines.f90 @@ -24,7 +24,7 @@ module system_routines function setCWD_C(cwd) bind(C) use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR - + integer(C_INT) :: setCWD_C character(kind=C_CHAR), dimension(*), intent(in) :: cwd end function setCWD_C @@ -150,14 +150,14 @@ function getUserName() getUserName = c_f_string(getUserName_Cstring) else getUserName = 'n/a (Error!)' - endif + end if end function getUserName !-------------------------------------------------------------------------------------------------- -!> @brief convert C string to Fortran string -!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string +!> @brief Convert C string to Fortran string. +!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string. !-------------------------------------------------------------------------------------------------- pure function c_f_string(c_string) result(f_string) @@ -174,28 +174,23 @@ pure function c_f_string(c_string) result(f_string) else f_string = f_string(:i-1) exit - endif - enddo arrayToString + end if + end do arrayToString end function c_f_string !-------------------------------------------------------------------------------------------------- -!> @brief convert Fortran string to C string -!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string +!> @brief Convert Fortran string to C string. +!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string. !-------------------------------------------------------------------------------------------------- pure function f_c_string(f_string) result(c_string) character(len=*), intent(in) :: f_string character(kind=C_CHAR), dimension(len_trim(f_string)+1) :: c_string - integer :: i - - do i=1,len_trim(f_string) - c_string(i)=f_string(i:i) - enddo - c_string(len_trim(f_string)+1) = C_NULL_CHAR + c_string = transfer(trim(f_string)//C_NULL_CHAR,c_string,size=size(c_string)) end function f_c_string