From 610c233fb68b1a87ef48d98d103ee66e4872799d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Dec 2020 08:45:08 +0100 Subject: [PATCH 01/17] ensure correct shape --- python/tests/test_Table.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_Table.py b/python/tests/test_Table.py index ac5859ecb..8f617aff5 100644 --- a/python/tests/test_Table.py +++ b/python/tests/test_Table.py @@ -22,7 +22,7 @@ class TestTable: @pytest.mark.parametrize('N',[10,40]) def test_len(self,N): - len(Table(np.random.rand(N,3),{'X':3})) == N + assert len(Table(np.random.rand(N,3),{'X':3})) == N def test_get_scalar(self,default): d = default.get('s') From 0e8082860cac913bc5739ff7ce9436295a257546 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Dec 2020 08:45:44 +0100 Subject: [PATCH 02/17] Fortran standard is 2018 will not work for older compilers --- cmake/Compiler-Intel.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 719ed885b..5b551069e 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -20,7 +20,7 @@ endif () # -assume std_mod_proc_name (included in -standard-semantics) causes problems if other modules # (PETSc, HDF5) are not compiled with this option (https://software.intel.com/en-us/forums/intel-fortran-compiler-for-linux-and-mac-os-x/topic/62172) -set (STANDARD_CHECK "-stand f15 -standard-semantics -assume nostd_mod_proc_name") +set (STANDARD_CHECK "-stand f18 -standard-semantics -assume nostd_mod_proc_name") set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel") # Link against shared Intel libraries instead of static ones From b6d00e2fb836a8324b85da3159a93234415d01ce Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Dec 2020 08:48:04 +0100 Subject: [PATCH 03/17] limit access to public variables to one function not sure if the 'volatile' attribute is needed --- src/DAMASK_interface.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index 41f421eb8..f664eb458 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -392,7 +392,7 @@ end function makeRelativePath subroutine catchSIGTERM(signal) bind(C) integer(C_INT), value :: signal - interface_SIGTERM = .true. + call interface_setSIGTERM(.true.) print'(a,i0,a)', ' received signal ',signal, ', set SIGTERM=TRUE' @@ -417,7 +417,7 @@ end subroutine interface_setSIGTERM subroutine catchSIGUSR1(signal) bind(C) integer(C_INT), value :: signal - interface_SIGUSR1 = .true. + call interface_setSIGUSR1(.true.) print'(a,i0,a)', ' received signal ',signal, ', set SIGUSR1=TRUE' @@ -442,7 +442,7 @@ end subroutine interface_setSIGUSR1 subroutine catchSIGUSR2(signal) bind(C) integer(C_INT), value :: signal - interface_SIGUSR2 = .true. + call interface_setSIGUSR2(.true.) print'(a,i0,a)', ' received signal ',signal, ', set SIGUSR2=TRUE' From 189597dbffd4a33bc64b55d4bce9ac49f638aa90 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Dec 2020 08:55:07 +0100 Subject: [PATCH 04/17] drop support for old PETSc versions --- src/DAMASK_interface.f90 | 2 +- src/mesh/mesh_mech_FEM.f90 | 29 +---------------------------- 2 files changed, 2 insertions(+), 29 deletions(-) diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index f664eb458..d38020225 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -10,7 +10,7 @@ !> and working directory. !-------------------------------------------------------------------------------------------------- #define PETSC_MAJOR 3 -#define PETSC_MINOR_MIN 10 +#define PETSC_MINOR_MIN 12 #define PETSC_MINOR_MAX 14 module DAMASK_interface diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 8aa084ac8..a4fa29204 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -147,14 +147,9 @@ subroutine FEM_mech_init(fieldBC) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) nBasis = nBasis/nc -#if (PETSC_VERSION_MINOR > 10) call DMAddField(mech_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr) call DMCreateDS(mech_mesh,ierr); CHKERRQ(ierr) -#endif call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) -#if (PETSC_VERSION_MINOR < 11) - call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr) -#endif call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr) call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr) call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr) @@ -163,11 +158,7 @@ subroutine FEM_mech_init(fieldBC) ! Setup FEM mech boundary conditions call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr) -#if (PETSC_VERSION_MINOR < 12) - call DMGetSection(mech_mesh,section,ierr); CHKERRQ(ierr) -#else call DMGetLocalSection(mech_mesh,section,ierr); CHKERRQ(ierr) -#endif allocate(pnumComp(1), source=dimPlex) allocate(pnumDof(0:dimPlex), source = 0) do topologDim = 0, dimPlex @@ -205,14 +196,8 @@ subroutine FEM_mech_init(fieldBC) endif endif enddo; enddo -#if (PETSC_VERSION_MINOR < 11) - call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, & - numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr) -#else call DMPlexCreateSection(mech_mesh,nolabel,pNumComp,pNumDof, & numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr) - -#endif CHKERRQ(ierr) call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr) do faceSet = 1, numBC @@ -267,11 +252,7 @@ subroutine FEM_mech_init(fieldBC) x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal) enddo px_scal => x_scal -#if (PETSC_VERSION_MINOR < 11) - call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr) -#else - call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,5,ierr) ! PETSc: cbee0a90b60958e5c50c89b1e41f4451dfa6008c -#endif + call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,5,ierr) CHKERRQ(ierr) enddo @@ -355,11 +336,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) allocate(pinvcellJ(dimPlex**2)) allocate(x_scal(cellDof)) -#if (PETSC_VERSION_MINOR < 12) - call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) -#else call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr) -#endif call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr) call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr) CHKERRQ(ierr) @@ -502,11 +479,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr) call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr) -#if (PETSC_VERSION_MINOR < 12) - call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) -#else call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr) -#endif call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) From 104fa167bdcbdce00b32b31170b5e82a5f1f4e55 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Dec 2020 19:30:34 +0100 Subject: [PATCH 05/17] missing rename: constituent -> phase meaningfull order --- src/material.f90 | 2 +- src/results.f90 | 169 ++++++++++++++++++++++++----------------------- 2 files changed, 87 insertions(+), 84 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index 223ea6ed8..b05979298 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -176,7 +176,7 @@ subroutine material_init(restart) if (.not. restart) then call results_openJobFile - call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,material_name_phase) + call results_mapping_phase(material_phaseAt,material_phaseMemberAt,material_name_phase) call results_mapping_homogenization(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization) call results_closeJobFile endif diff --git a/src/results.f90 b/src/results.f90 index f15ad4e4a..ea9fd62d4 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -49,7 +49,7 @@ module results results_setLink, & results_addAttribute, & results_removeLink, & - results_mapping_constituent, & + results_mapping_phase, & results_mapping_homogenization contains @@ -461,7 +461,7 @@ end subroutine results_writeTensorDataset_int !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) +subroutine results_mapping_phase(phaseAt,memberAtLocal,label) integer, dimension(:,:), intent(in) :: phaseAt !< phase section at (constituent,element) integer, dimension(:,:,:), intent(in) :: memberAtLocal !< phase member at (constituent,IP,element) @@ -491,6 +491,47 @@ subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) integer(SIZE_T) :: type_size_string, type_size_int integer :: hdferr, ierr, i +!-------------------------------------------------------------------------------------------------- +! prepare MPI communication (transparent for non-MPI runs) + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + memberOffset = 0 + do i=1, size(label) + memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAtLocal,2) ! number of points/instance of this process + enddo + writeSize = 0 + writeSize(worldrank) = size(memberAtLocal(1,:,:)) ! total number of points by this process + +!-------------------------------------------------------------------------------------------------- +! MPI settings and communication +#ifdef PETSc + 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,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process + if(ierr /= 0) error stop 'MPI error' + + call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process + if(ierr /= 0) error stop 'MPI error' +#endif + + myShape = int([size(phaseAt,1),writeSize(worldrank)], HSIZE_T) + myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T) + totalShape = int([size(phaseAt,1),sum(writeSize)], HSIZE_T) + + +!--------------------------------------------------------------------------------------------------- +! expand phaseAt to consider IPs (is not stored per IP) + do i = 1, size(phaseAtMaterialpoint,2) + phaseAtMaterialpoint(:,i,:) = phaseAt + enddo + +!--------------------------------------------------------------------------------------------------- +! renumber member from my process to all processes + do i = 1, size(label) + where(phaseAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1 ! convert to 0-based + enddo + !--------------------------------------------------------------------------------------------------- ! compound type: name of phase section + position/index within results array call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) @@ -525,34 +566,6 @@ subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) call h5tclose_f(dt_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' -!-------------------------------------------------------------------------------------------------- -! prepare MPI communication (transparent for non-MPI runs) - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - memberOffset = 0 - do i=1, size(label) - memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAtLocal,2) ! number of points/instance of this process - enddo - writeSize = 0 - writeSize(worldrank) = size(memberAtLocal(1,:,:)) ! total number of points by this process - -!-------------------------------------------------------------------------------------------------- -! MPI settings and communication -#ifdef PETSc - 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,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process - if(ierr /= 0) error stop 'MPI error' - - call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process - if(ierr /= 0) error stop 'MPI error' -#endif - - myShape = int([size(phaseAt,1),writeSize(worldrank)], HSIZE_T) - myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T) - totalShape = int([size(phaseAt,1),sum(writeSize)], HSIZE_T) - !-------------------------------------------------------------------------------------------------- ! create dataspace in memory (local shape = hyperslab) and in file (global shape) call h5screate_simple_f(2,myShape,memspace_id,hdferr,myShape) @@ -564,18 +577,6 @@ subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr) if(hdferr < 0) error stop 'HDF5 error' -!--------------------------------------------------------------------------------------------------- -! expand phaseAt to consider IPs (is not stored per IP) - do i = 1, size(phaseAtMaterialpoint,2) - phaseAtMaterialpoint(:,i,:) = phaseAt - enddo - -!--------------------------------------------------------------------------------------------------- -! renumber member from my process to all processes - do i = 1, size(label) - where(phaseAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1 ! convert to 0-based - enddo - !-------------------------------------------------------------------------------------------------- ! write the components of the compound type individually call h5pset_preserve_f(plist_id, .TRUE., hdferr) @@ -609,7 +610,7 @@ subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) if(hdferr < 0) error stop 'HDF5 error' call h5tclose_f(position_id, hdferr) -end subroutine results_mapping_constituent +end subroutine results_mapping_phase !-------------------------------------------------------------------------------------------------- @@ -645,6 +646,48 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) integer(SIZE_T) :: type_size_string, type_size_int integer :: hdferr, ierr, i + +!-------------------------------------------------------------------------------------------------- +! prepare MPI communication (transparent for non-MPI runs) + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + memberOffset = 0 + do i=1, size(label) + memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAtLocal,1) ! number of points/instance of this process + enddo + writeSize = 0 + writeSize(worldrank) = size(memberAtLocal) ! total number of points by this process + +!-------------------------------------------------------------------------------------------------- +! MPI settings and communication +#ifdef PETSc + 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,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process + if(ierr /= 0) error stop 'MPI error' + + call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process + if(ierr /= 0) error stop 'MPI error' +#endif + + myShape = int([writeSize(worldrank)], HSIZE_T) + myOffset = int([sum(writeSize(0:worldrank-1))], HSIZE_T) + totalShape = int([sum(writeSize)], HSIZE_T) + + +!--------------------------------------------------------------------------------------------------- +! expand phaseAt to consider IPs (is not stored per IP) + do i = 1, size(homogenizationAtMaterialpoint,1) + homogenizationAtMaterialpoint(i,:) = homogenizationAt + enddo + +!--------------------------------------------------------------------------------------------------- +! renumber member from my process to all processes + do i = 1, size(label) + where(homogenizationAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1 ! convert to 0-based + enddo + !--------------------------------------------------------------------------------------------------- ! compound type: name of phase section + position/index within results array call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) @@ -679,34 +722,6 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) call h5tclose_f(dt_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' -!-------------------------------------------------------------------------------------------------- -! prepare MPI communication (transparent for non-MPI runs) - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - memberOffset = 0 - do i=1, size(label) - memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAtLocal,1) ! number of points/instance of this process - enddo - writeSize = 0 - writeSize(worldrank) = size(memberAtLocal) ! total number of points by this process - -!-------------------------------------------------------------------------------------------------- -! MPI settings and communication -#ifdef PETSc - 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,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process - if(ierr /= 0) error stop 'MPI error' - - call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process - if(ierr /= 0) error stop 'MPI error' -#endif - - myShape = int([writeSize(worldrank)], HSIZE_T) - myOffset = int([sum(writeSize(0:worldrank-1))], HSIZE_T) - totalShape = int([sum(writeSize)], HSIZE_T) - !-------------------------------------------------------------------------------------------------- ! create dataspace in memory (local shape = hyperslab) and in file (global shape) call h5screate_simple_f(1,myShape,memspace_id,hdferr,myShape) @@ -718,18 +733,6 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr) if(hdferr < 0) error stop 'HDF5 error' -!--------------------------------------------------------------------------------------------------- -! expand phaseAt to consider IPs (is not stored per IP) - do i = 1, size(homogenizationAtMaterialpoint,1) - homogenizationAtMaterialpoint(i,:) = homogenizationAt - enddo - -!--------------------------------------------------------------------------------------------------- -! renumber member from my process to all processes - do i = 1, size(label) - where(homogenizationAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1 ! convert to 0-based - enddo - !-------------------------------------------------------------------------------------------------- ! write the components of the compound type individually call h5pset_preserve_f(plist_id, .TRUE., hdferr) From 8100f3ebfaa78da64f473d721016dc0e84abd8de Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Dec 2020 22:29:36 +0100 Subject: [PATCH 06/17] not required for cmake > 2.4 --- CMakeLists.txt | 2 +- PRIVATE | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index dd2348fd1..d8eb2cbf0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ ######################################################################################## # Compiler options for building DAMASK -cmake_minimum_required (VERSION 3.10.0 FATAL_ERROR) +cmake_minimum_required (VERSION 3.10.0) #--------------------------------------------------------------------------------------- # Find PETSc from system environment diff --git a/PRIVATE b/PRIVATE index 08f8aea46..6b8ba6d84 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 08f8aea465a1b5e476b584bcae7927d113919b1d +Subproject commit 6b8ba6d844b70695233b02eae85f4300315f339d From b95eb6f604dfef0bb08b4d7936f857f635342212 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Dec 2020 13:01:52 +0100 Subject: [PATCH 07/17] simplified (using pkg-config module of PETSc) might be possible to use pkg-config for FFTW and HDF5 in future --- CMakeLists.txt | 53 ++++++++++++++++---------------------------------- 1 file changed, 17 insertions(+), 36 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d8eb2cbf0..553c42400 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,18 @@ -######################################################################################## -# Compiler options for building DAMASK cmake_minimum_required (VERSION 3.10.0) +include (FindPkgConfig REQUIRED) + +# Dummy project to determine compiler names and version +project (Prerequisites LANGUAGES) +set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig") +pkg_search_module (PETSC REQUIRED PETSc>3.12.0) +pkg_get_variable (CMAKE_Fortran_COMPILER PETSc fcompiler) +pkg_get_variable (CMAKE_C_COMPILER PETSc ccompiler) + +find_program (CAT_EXECUTABLE NAMES cat) +execute_process (COMMAND ${CAT_EXECUTABLE} ${PROJECT_SOURCE_DIR}/VERSION + RESULT_VARIABLE DAMASK_VERSION_RETURN + OUTPUT_VARIABLE DAMASK_VERSION + OUTPUT_STRIP_TRAILING_WHITESPACE) #--------------------------------------------------------------------------------------- # Find PETSc from system environment @@ -28,17 +40,11 @@ include ${petsc_conf_rules} include ${petsc_conf_variables} INCLUDE_DIRS := \${PETSC_FC_INCLUDES} LIBRARIES := \${PETSC_WITH_EXTERNAL_LIB} -COMPILERF := \${FC} -COMPILERC := \${CC} LINKERNAME := \${FLINKER} includes: \t@echo \${INCLUDE_DIRS} extlibs: \t@echo \${LIBRARIES} -compilerf: -\t@echo \${COMPILERF} -compilerc: -\t@echo \${COMPILERC} linker: \t@echo \${LINKERNAME} ") @@ -57,16 +63,6 @@ execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_conf RESULT_VARIABLE PETSC_EXTERNAL_LIB_RETURN OUTPUT_VARIABLE petsc_external_lib OUTPUT_STRIP_TRAILING_WHITESPACE) -# PETSc specified fortran compiler -execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "compilerf" - RESULT_VARIABLE PETSC_MPIFC_RETURN - OUTPUT_VARIABLE PETSC_MPIFC - OUTPUT_STRIP_TRAILING_WHITESPACE) -# PETSc specified C compiler -execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "compilerc" - RESULT_VARIABLE PETSC_MPICC_RETURN - OUTPUT_VARIABLE PETSC_MPICC - OUTPUT_STRIP_TRAILING_WHITESPACE) # PETSc specified linker (Fortran compiler + PETSc linking flags) execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "linker" RESULT_VARIABLE PETSC_LINKER_RETURN @@ -91,13 +87,6 @@ message ("Found PETSC_DIR:\n${PETSC_DIR}\n" ) message ("Found PETSC_INCLUDES:\n${PETSC_INCLUDES}\n" ) message ("Found PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n") message ("Found PETSC_LINKER:\n${PETSC_LINKER}\n" ) -message ("Found MPI Fortran Compiler:\n${PETSC_MPIFC}\n" ) -message ("Found MPI C Compiler:\n${PETSC_MPICC}\n" ) - -# set compiler commands to match PETSc (needs to be done before defining the project) -# https://cmake.org/Wiki/CMake_FAQ#How_do_I_use_a_different_compiler.3F -set (CMAKE_Fortran_COMPILER "${PETSC_MPIFC}") -set (CMAKE_C_COMPILER "${PETSC_MPICC}") #--------------------------------------------------------------------------------------- # Now start to care about DAMASK @@ -115,7 +104,8 @@ elseif (DAMASK_SOLVER STREQUAL "mesh") else () message (FATAL_ERROR "Build target (DAMASK_SOLVER) is not defined") endif () -list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) +add_definitions (-DDAMASKVERSION="${DAMASK_VERSION}") +add_definitions (-DPETSc) if (CMAKE_BUILD_TYPE STREQUAL "") set (CMAKE_BUILD_TYPE "RELEASE") @@ -153,17 +143,8 @@ if (CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") set (BUILDCMD_POST "${BUILDCMD_POST} -fsyntax-only") endif () -# Parse DAMASK version from VERSION file -find_program (CAT_EXECUTABLE NAMES cat) -execute_process (COMMAND ${CAT_EXECUTABLE} ${PROJECT_SOURCE_DIR}/VERSION - RESULT_VARIABLE DAMASK_VERSION_RETURN - OUTPUT_VARIABLE DAMASK_V - OUTPUT_STRIP_TRAILING_WHITESPACE) -add_definitions (-DDAMASKVERSION="${DAMASK_V}") - -# definition of other macros -add_definitions (-DPETSc) +list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel") include (Compiler-Intel) elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") From 168d6e85e1ddd948d53b7825470b9ba1f7f7bd59 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Dec 2020 17:44:31 +0100 Subject: [PATCH 08/17] simplified --- CMakeLists.txt | 21 +++++---------------- 1 file changed, 5 insertions(+), 16 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 553c42400..8db6dd0c0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -40,13 +40,10 @@ include ${petsc_conf_rules} include ${petsc_conf_variables} INCLUDE_DIRS := \${PETSC_FC_INCLUDES} LIBRARIES := \${PETSC_WITH_EXTERNAL_LIB} -LINKERNAME := \${FLINKER} includes: \t@echo \${INCLUDE_DIRS} extlibs: \t@echo \${LIBRARIES} -linker: -\t@echo \${LINKERNAME} ") # CMake will execute each target in the ${petsc_config_makefile} @@ -58,16 +55,10 @@ execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_conf OUTPUT_VARIABLE petsc_includes OUTPUT_STRIP_TRAILING_WHITESPACE) # Find the PETSc external linking directory settings -# required for final linking, must be appended after the executable execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "extlibs" RESULT_VARIABLE PETSC_EXTERNAL_LIB_RETURN OUTPUT_VARIABLE petsc_external_lib OUTPUT_STRIP_TRAILING_WHITESPACE) -# PETSc specified linker (Fortran compiler + PETSc linking flags) -execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "linker" - RESULT_VARIABLE PETSC_LINKER_RETURN - OUTPUT_VARIABLE PETSC_LINKER - OUTPUT_STRIP_TRAILING_WHITESPACE) # Remove temporary makefile, no need to keep it anymore. file (REMOVE_RECURSE ${TEMPDIR}) @@ -86,7 +77,6 @@ endforeach (exlib) message ("Found PETSC_DIR:\n${PETSC_DIR}\n" ) message ("Found PETSC_INCLUDES:\n${PETSC_INCLUDES}\n" ) message ("Found PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n") -message ("Found PETSC_LINKER:\n${PETSC_LINKER}\n" ) #--------------------------------------------------------------------------------------- # Now start to care about DAMASK @@ -94,19 +84,19 @@ message ("Found PETSC_LINKER:\n${PETSC_LINKER}\n" ) # DAMASK solver defines project to build string(TOLOWER ${DAMASK_SOLVER} DAMASK_SOLVER) if (DAMASK_SOLVER STREQUAL "grid") - project (damask-grid Fortran C) + project (damask-grid HOMEPAGE_URL https://damask.mpie.de LANGUAGES Fortran C) add_definitions (-DGrid) - message ("Building Grid Solver\n") elseif (DAMASK_SOLVER STREQUAL "mesh") - project (damask-mesh Fortran C) + project (damask-mesh HOMEPAGE_URL https://damask.mpie.de LANGUAGES Fortran C) add_definitions (-DMesh) - message ("Building Mesh Solver\n") else () message (FATAL_ERROR "Build target (DAMASK_SOLVER) is not defined") endif () add_definitions (-DDAMASKVERSION="${DAMASK_VERSION}") add_definitions (-DPETSc) +message ("\nBuilding ${CMAKE_PROJECT_NAME}\n") + if (CMAKE_BUILD_TYPE STREQUAL "") set (CMAKE_BUILD_TYPE "RELEASE") endif () @@ -155,9 +145,8 @@ else () message (FATAL_ERROR "Compiler type (CMAKE_Fortran_COMPILER_ID) not recognized") endif () - set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}") -set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${PETSC_LINKER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}") +set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${CMAKE_Fortran_COMPILER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}") if (CMAKE_BUILD_TYPE STREQUAL "DEBUG") set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}") From 8fbadef52476d4e523acf2689e5ebdba396876df Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Dec 2020 22:37:14 +0100 Subject: [PATCH 09/17] print instead of write --- src/grid/grid_damage_spectral.f90 | 3 +-- src/grid/grid_mech_FEM.f90 | 9 ++++---- src/grid/grid_mech_spectral_basic.f90 | 9 ++++---- src/grid/grid_mech_spectral_polarisation.f90 | 9 ++++---- src/grid/grid_thermal_spectral.f90 | 3 +-- src/grid/spectral_utilities.f90 | 23 +++++++++----------- src/mesh/mesh_mech_FEM.f90 | 4 ++-- 7 files changed, 26 insertions(+), 34 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 4c014f3c0..79437945b 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -203,8 +203,7 @@ function grid_damage_spectral_solution(timeinc) result(solution) call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr) if (solution%converged) & print'(/,a)', ' ... nonlocal damage converged .....................................' - write(IO_STDOUT,'(/,a,f8.6,2x,f8.6,2x,e11.4,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& - phi_min, phi_max, stagNorm + print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm print'(/,a)', ' ===========================================================================' flush(IO_STDOUT) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 4394b6f81..8874b0cf3 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -510,11 +510,10 @@ subroutine formResidual(da_local,x_local, & newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax - if (debugRotation) & - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) + if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) + print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim =', transpose(F_aim) flush(IO_STDOUT) endif newIteration diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 563b25162..05986f32e 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -472,11 +472,10 @@ subroutine formResidual(in, F, & newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax - if (debugRotation) & - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) + if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) + print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim =', transpose(F_aim) flush(IO_STDOUT) endif newIteration diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 03780f2e0..832441a4c 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -555,11 +555,10 @@ subroutine formResidual(in, FandF_tau, & newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax - if(debugRotation) & - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) - write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) + if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) + print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + ' deformation gradient aim =', transpose(F_aim) flush(IO_STDOUT) endif newIteration diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 68a1c5ed1..f5d1a33bc 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -197,8 +197,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution) call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr) if (solution%converged) & print'(/,a)', ' ... thermal conduction converged ..................................' - write(IO_STDOUT,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',& - T_min, T_max, stagNorm + print'(/,a,f8.4,2x,f8.4,2x,f8.4)', ' Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm print'(/,a)', ' ===========================================================================' flush(IO_STDOUT) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 1e1608d7c..27fe7fd6e 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -688,8 +688,8 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) if(debugGeneral) then print'(/,a)', ' ... updating masked compliance ............................................' - write(IO_STDOUT,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& - transpose(temp99_Real)*1.0e-9_pReal + print'(/,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', & + ' Stiffness C (load) / GPa =', transpose(temp99_Real)*1.0e-9_pReal flush(IO_STDOUT) endif @@ -709,9 +709,8 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) if (debugGeneral .or. errmatinv) then write(formatString, '(i2)') size_reduced formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' - write(IO_STDOUT,trim(formatString),advance='no') ' C * S (load) ', & - transpose(matmul(c_reduced,s_reduced)) - write(IO_STDOUT,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced) + print trim(formatString), ' C * S (load) ', transpose(matmul(c_reduced,s_reduced)) + print trim(formatString), ' S (load) ', transpose(s_reduced) if(errmatinv) error stop 'matrix inversion error' endif temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9]) @@ -722,7 +721,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) utilities_maskedCompliance = math_99to3333(temp99_Real) if(debugGeneral) then - write(IO_STDOUT,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') & + print'(/,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', & ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal flush(IO_STDOUT) endif @@ -818,13 +817,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if (debugRotation) & - write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& - transpose(P_av)*1.e-6_pReal - if(present(rotation_BC)) & - P_av = rotation_BC%rotate(P_av) - write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& - transpose(P_av)*1.e-6_pReal + if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & + ' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal + if(present(rotation_BC)) P_av = rotation_BC%rotate(P_av) + print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & + ' Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pReal flush(IO_STDOUT) dPdF_max = 0.0_pReal diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index a4fa29204..14ce1f38f 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -659,8 +659,8 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), & ' @ Iteration ',PETScIter,' mechanical residual norm = ', & int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) - write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& - transpose(P_av)*1.e-6_pReal + print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & + ' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal flush(IO_STDOUT) end subroutine FEM_mech_converged From 872ceac855fa1e5e50d388f0a0fb8de992b41baf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Dec 2020 11:26:31 +0100 Subject: [PATCH 10/17] not needed --- src/crystallite.f90 | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 66c1df607..3e8b90c38 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -135,7 +135,6 @@ contains !-------------------------------------------------------------------------------------------------- subroutine crystallite_init - logical, dimension(discretization_nIPs,discretization_Nelems) :: devNull integer :: & c, & !< counter in integration point component loop i, & !< counter in integration point loop @@ -288,8 +287,6 @@ subroutine crystallite_init enddo !$OMP END PARALLEL DO - devNull = crystallite_stress() - #ifdef DEBUG if (debugCrystallite%basic) then print'(a42,1x,i10)', ' # of elements: ', eMax @@ -321,14 +318,11 @@ function crystallite_stress() subLp0,& !< plastic velocity grad at start of crystallite inc subLi0 !< intermediate velocity grad at start of crystallite inc - todo = .false. subLp0 = crystallite_partitionedLp0 subLi0 = crystallite_partitionedLi0 - - !-------------------------------------------------------------------------------------------------- ! initialize to starting condition crystallite_subStep = 0.0_pReal @@ -435,8 +429,6 @@ function crystallite_stress() ! integrate --- requires fully defined state array (basic + dependent state) where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation - - enddo cutbackLooping ! return whether converged or not @@ -471,10 +463,10 @@ subroutine crystallite_initializeRestorationPoints(i,e) crystallite_partitionedS0(1:3,1:3,c,i,e) = crystallite_S0(1:3,1:3,c,i,e) plasticState(material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = & - plasticState(material_phaseAt(c,e))%state0( :,material_phasememberAt(c,i,e)) + plasticState(material_phaseAt(c,e))%state0( :,material_phasememberAt(c,i,e)) do s = 1, phase_Nsources(material_phaseAt(c,e)) sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = & - sourceState(material_phaseAt(c,e))%p(s)%state0( :,material_phasememberAt(c,i,e)) + sourceState(material_phaseAt(c,e))%p(s)%state0( :,material_phasememberAt(c,i,e)) enddo enddo From d7f035235c5052125ddd384d952a301ed30998a7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Dec 2020 14:01:39 +0100 Subject: [PATCH 11/17] do initialization later --- src/CPFEM.f90 | 1 + src/CPFEM2.f90 | 1 + src/crystallite.f90 | 44 +++++++++++++++++++++++++++----------------- 3 files changed, 29 insertions(+), 17 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index c49ecbcb6..58c17554a 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -93,6 +93,7 @@ subroutine CPFEM_initAll call homogenization_init call CPFEM_init call config_deallocate + call crystallite_setInitialValues ! ToDo: MD More general approach needed end subroutine CPFEM_initAll diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 994859758..5d9d24149 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -67,6 +67,7 @@ subroutine CPFEM_initAll call homogenization_init call CPFEM_init call config_deallocate + call crystallite_setInitialValues ! ToDo: MD More general approach needed end subroutine CPFEM_initAll diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 3e8b90c38..0fcade113 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -115,6 +115,7 @@ module crystallite public :: & crystallite_init, & + crystallite_setInitialValues, & crystallite_stress, & crystallite_stressTangent, & crystallite_orientations, & @@ -136,9 +137,7 @@ contains subroutine crystallite_init integer :: & - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e, & !< counter in element loop + p, & cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax !< maximum number of elements @@ -237,19 +236,38 @@ subroutine crystallite_init phases => config_material%get('phase') allocate(output_constituent(phases%length)) - do c = 1, phases%length - phase => phases%get(c) + do p = 1, phases%length + phase => phases%get(p) mech => phase%get('mechanics',defaultVal = emptyDict) #if defined(__GFORTRAN__) - output_constituent(c)%label = output_asStrings(mech) + output_constituent(p)%label = output_asStrings(mech) #else - output_constituent(c)%label = mech%get_asStrings('output',defaultVal=emptyStringArray) + output_constituent(p)%label = mech%get_asStrings('output',defaultVal=emptyStringArray) #endif enddo +#ifdef DEBUG + if (debugCrystallite%basic) then + 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) + endif +#endif + +end subroutine crystallite_init + !-------------------------------------------------------------------------------------------------- -! initialize +!> @brief Set initial values +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_setInitialValues + + integer :: & + c, & !< counter in integration point component loop + i, & !< counter in integration point loop + e !< counter in element loop + !$OMP PARALLEL DO PRIVATE(i,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Nconstituents(material_homogenizationAt(e)) @@ -287,16 +305,8 @@ subroutine crystallite_init enddo !$OMP END PARALLEL DO -#ifdef DEBUG - if (debugCrystallite%basic) then - 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) - endif -#endif -end subroutine crystallite_init +end subroutine crystallite_setInitialValues !-------------------------------------------------------------------------------------------------- From f8756ad95ae71128c7c9fca0239f9cf9887c72bb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Dec 2020 17:45:11 +0100 Subject: [PATCH 12/17] simplifying no extral handling for homogeneous temperature (the memory that was saved was consumed by the extra mapping) --- src/CPFEM.f90 | 2 +- src/constitutive.f90 | 24 ++++++++++------------- src/constitutive_mech.f90 | 10 ++++++++-- src/damage_local.f90 | 8 +++----- src/damage_none.f90 | 5 +++-- src/damage_nonlocal.f90 | 8 +++----- src/grid/grid_thermal_spectral.f90 | 3 +-- src/kinematics_thermal_expansion.f90 | 4 ++-- src/material.f90 | 29 +++------------------------- src/thermal_adiabatic.f90 | 11 ++++------- src/thermal_conduction.f90 | 5 +---- src/thermal_isothermal.f90 | 7 +++---- 12 files changed, 42 insertions(+), 74 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 58c17554a..eb3243576 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -182,7 +182,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) case (THERMAL_conduction_ID) chosenThermal1 - temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & + temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & temperature_inp end select chosenThermal1 homogenization_F0(1:3,1:3,ip,elCP) = ffn diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 470ac4dc7..a60352f81 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -19,25 +19,21 @@ module constitutive implicit none private - integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable :: & !ToDo: old intel compiler complains about protected - phase_elasticity !< elasticity of each phase - - integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & !ToDo: old intel compiler complains about protected + integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & phase_plasticity !< plasticity of each phase - integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & ! ToDo: old intel compiler complains about protected + integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & phase_source, & !< active sources mechanisms of each phase - phase_kinematics, & !< active kinematic mechanisms of each phase - phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase + phase_kinematics !< active kinematic mechanisms of each phase - integer, dimension(:), allocatable, public :: & ! ToDo: old intel compiler complains about protected + integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) phase_Nsources, & !< number of source mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase phase_plasticityInstance, & !< instance of particular plasticity of each phase phase_elasticityInstance !< instance of particular elasticity of each phase - logical, dimension(:), allocatable, public :: & ! ToDo: old intel compiler complains about protected + logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) phase_localPlasticity !< flags phases with local constitutive law type(tPlasticState), allocatable, dimension(:), public :: & @@ -634,10 +630,10 @@ pure function constitutive_initialFi(ipc, ip, el) KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption kinematicsType: select case (phase_kinematics(k,phase)) case (KINEMATICS_thermal_expansion_ID) kinematicsType - homog = material_homogenizationAt(el) - offset = thermalMapping(homog)%p(ip,el) - constitutive_initialFi = & - constitutive_initialFi + kinematics_thermal_expansion_initialStrain(homog,phase,offset) + homog = material_homogenizationAt(el) + offset = material_homogenizationMemberAt(ip,el) + constitutive_initialFi = constitutive_initialFi & + + kinematics_thermal_expansion_initialStrain(homog,phase,offset) end select kinematicsType enddo KinematicsLoop @@ -674,7 +670,7 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el logical :: broken ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) + tme = material_homogenizationMemberAt(ip,el) instance = phase_plasticityInstance(phase) Mp = matmul(matmul(transpose(Fi),Fi),S) diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index dc3a935e3..d64ea327c 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -3,6 +3,12 @@ !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_mech + integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable :: & + phase_elasticity !< elasticity of each phase + integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & + phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase + + interface module function plastic_none_init() result(myPlasticity) @@ -360,7 +366,7 @@ module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) instance, of ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) + tme = material_homogenizationMemberAt(ip,el) of = material_phasememberAt(ipc,ip,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) @@ -407,7 +413,7 @@ module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & i, j, instance, of ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) + tme = material_homogenizationMemberAt(ip,el) Mp = matmul(matmul(transpose(Fi),Fi),S) of = material_phasememberAt(ipc,ip,el) diff --git a/src/damage_local.f90 b/src/damage_local.f90 index e63db90b0..0003bb2a8 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -75,13 +75,11 @@ subroutine damage_local_init Nmaterialpoints = count(material_homogenizationAt == h) damageState(h)%sizeState = 1 - allocate(damageState(h)%state0 (1,Nmaterialpoints), source=damage_initialPhi(h)) - allocate(damageState(h)%subState0(1,Nmaterialpoints), source=damage_initialPhi(h)) - allocate(damageState(h)%state (1,Nmaterialpoints), source=damage_initialPhi(h)) + allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal) - nullify(damageMapping(h)%p) damageMapping(h)%p => material_homogenizationMemberAt - deallocate(damage(h)%p) damage(h)%p => damageState(h)%state(1,:) end associate diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 2279bc06b..1cf3de6de 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -3,6 +3,7 @@ !> @brief material subroutine for constant damage field !-------------------------------------------------------------------------------------------------- module damage_none + use prec use config use material @@ -29,8 +30,8 @@ subroutine damage_none_init allocate(damageState(h)%subState0(0,Nmaterialpoints)) allocate(damageState(h)%state (0,Nmaterialpoints)) - deallocate(damage(h)%p) - allocate (damage(h)%p(1), source=damage_initialPhi(h)) + damageMapping(h)%p => material_homogenizationMemberAt + allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal) enddo diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 24a51cf54..ce97ec24a 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -78,13 +78,11 @@ subroutine damage_nonlocal_init Nmaterialpoints = count(material_homogenizationAt == h) damageState(h)%sizeState = 1 - allocate(damageState(h)%state0 (1,Nmaterialpoints), source=damage_initialPhi(h)) - allocate(damageState(h)%subState0(1,Nmaterialpoints), source=damage_initialPhi(h)) - allocate(damageState(h)%state (1,Nmaterialpoints), source=damage_initialPhi(h)) + allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal) - nullify(damageMapping(h)%p) damageMapping(h)%p => material_homogenizationMemberAt - deallocate(damage(h)%p) damage(h)%p => damageState(h)%state(1,:) end associate diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index f5d1a33bc..259b45f33 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -131,8 +131,7 @@ subroutine grid_thermal_spectral_init cell = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 - T_current(i,j,k) = temperature(material_homogenizationAt(cell))% & - p(thermalMapping(material_homogenizationAt(cell))%p(1,cell)) + T_current(i,j,k) = temperature(material_homogenizationAt(cell))%p(material_homogenizationMemberAt(1,cell)) T_lastInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k) enddo; enddo; enddo diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 2b8b04d85..36a882a48 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -126,8 +126,8 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, i phase = material_phaseAt(ipc,el) homog = material_homogenizationAt(el) - T = temperature(homog)%p(thermalMapping(homog)%p(ip,el)) - TDot = temperatureRate(homog)%p(thermalMapping(homog)%p(ip,el)) + T = temperature(homog)%p(material_homogenizationMemberAt(ip,el)) + TDot = temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el)) associate(prm => param(kinematics_thermal_expansion_instance(phase))) Li = TDot * ( & diff --git a/src/material.f90 b/src/material.f90 index b05979298..33b64e6df 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -64,17 +64,16 @@ module material homogenization_type !< type of each homogenization integer, public, protected :: & - homogenization_maxNconstituents !< max number of grains in any USED homogenization + homogenization_maxNconstituents !< max number of grains in any USED homogenization integer, dimension(:), allocatable, public, protected :: & - homogenization_Nconstituents, & !< number of grains in each homogenization + homogenization_Nconstituents, & !< number of grains in each homogenization homogenization_typeInstance, & !< instance of particular type of each homogenization thermal_typeInstance, & !< instance of particular type of each thermal transport damage_typeInstance !< instance of particular type of each nonlocal damage real(pReal), dimension(:), allocatable, public, protected :: & - thermal_initialT, & !< initial temperature per each homogenization - damage_initialPhi !< initial damage per each homogenization + thermal_initialT !< initial temperature per each homogenization integer, dimension(:), allocatable, public, protected :: & ! (elem) material_homogenizationAt !< homogenization ID of each element @@ -93,12 +92,7 @@ module material type(Rotation), dimension(:,:,:), allocatable, public, protected :: & material_orientation0 !< initial orientation of each grain,IP,element -! BEGIN DEPRECATED - integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field -! END DEPRECATED - type(tHomogMapping), allocatable, dimension(:), public :: & - thermalMapping, & !< mapping for thermal state/fields damageMapping !< mapping for damage state/fields type(group_float), allocatable, dimension(:), public :: & @@ -165,7 +159,6 @@ subroutine material_init(restart) allocate(thermalState (size(material_name_homogenization))) allocate(damageState (size(material_name_homogenization))) - allocate(thermalMapping (size(material_name_homogenization))) allocate(damageMapping (size(material_name_homogenization))) allocate(temperature (size(material_name_homogenization))) @@ -181,20 +174,6 @@ subroutine material_init(restart) call results_closeJobFile endif -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! BEGIN DEPRECATED - allocate(mappingHomogenizationConst( discretization_nIPs,discretization_Nelems),source=1) - -! hack needed to initialize field values used during constitutive initialization - do myHomog = 1, size(material_name_homogenization) - thermalMapping (myHomog)%p => mappingHomogenizationConst - damageMapping (myHomog)%p => mappingHomogenizationConst - allocate(temperature (myHomog)%p(1), source=thermal_initialT(myHomog)) - allocate(damage (myHomog)%p(1), source=damage_initialPhi(myHomog)) - allocate(temperatureRate (myHomog)%p(1), source=0.0_pReal) - enddo -! END DEPRECATED - end subroutine material_init @@ -222,7 +201,6 @@ subroutine material_parseHomogenization allocate(thermal_typeInstance(size(material_name_homogenization)), source=0) allocate(damage_typeInstance(size(material_name_homogenization)), source=0) allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal) - allocate(damage_initialPhi(size(material_name_homogenization)), source=1.0_pReal) do h=1, size(material_name_homogenization) homog => material_homogenization%get(h) @@ -258,7 +236,6 @@ subroutine material_parseHomogenization if(homog%contains('damage')) then homogDamage => homog%get('damage') - damage_initialPhi(h) = homogDamage%get_asFloat('phi_0',defaultVal=1.0_pReal) select case (homogDamage%get_asString('type')) case('none') damage_type(h) = DAMAGE_none_ID diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index aa807924c..c67d004bf 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -72,12 +72,9 @@ subroutine thermal_adiabatic_init allocate(thermalState(h)%state0 (1,Nmaterialpoints), source=thermal_initialT(h)) allocate(thermalState(h)%subState0(1,Nmaterialpoints), source=thermal_initialT(h)) allocate(thermalState(h)%state (1,Nmaterialpoints), source=thermal_initialT(h)) - - thermalMapping(h)%p => material_homogenizationMemberAt - deallocate(temperature(h)%p) + temperature(h)%p => thermalState(h)%state(1,:) - deallocate(temperatureRate(h)%p) - allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal) + allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal) end associate enddo @@ -117,8 +114,8 @@ function thermal_adiabatic_updateState(subdt, ip, el) <= 1.0e-6_pReal*abs(thermalState(homog)%state(1,offset)), & .true.] - temperature (homog)%p(thermalMapping(homog)%p(ip,el)) = T - temperatureRate(homog)%p(thermalMapping(homog)%p(ip,el)) = & + temperature (homog)%p(material_homogenizationMemberAt(ip,el)) = T + temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el)) = & (thermalState(homog)%state(1,offset) - thermalState(homog)%subState0(1,offset))/(subdt+tiny(0.0_pReal)) end function thermal_adiabatic_updateState diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index daa7391a9..602bdab35 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -71,10 +71,7 @@ subroutine thermal_conduction_init allocate(thermalState(h)%subState0(0,Nmaterialpoints)) allocate(thermalState(h)%state (0,Nmaterialpoints)) - thermalMapping(h)%p => material_homogenizationMemberAt - deallocate(temperature (h)%p) allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h)) - deallocate(temperatureRate(h)%p) allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal) end associate @@ -205,7 +202,7 @@ subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el) offset homog = material_homogenizationAt(el) - offset = thermalMapping(homog)%p(ip,el) + offset = material_homogenizationMemberAt(ip,el) temperature (homog)%p(offset) = T temperatureRate(homog)%p(offset) = Tdot diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index 39c8efe91..adf2257de 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -3,6 +3,7 @@ !> @brief material subroutine for isothermal temperature field !-------------------------------------------------------------------------------------------------- module thermal_isothermal + use prec use config use material @@ -29,10 +30,8 @@ subroutine thermal_isothermal_init allocate(thermalState(h)%subState0(0,Nmaterialpoints)) allocate(thermalState(h)%state (0,Nmaterialpoints)) - deallocate(temperature (h)%p) - allocate (temperature (h)%p(1), source=thermal_initialT(h)) - deallocate(temperatureRate(h)%p) - allocate (temperatureRate(h)%p(1)) + allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h)) + allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal) enddo From f8e3cfe91d0334073d7a2a1c193883e5c913a13e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Dec 2020 19:41:47 +0100 Subject: [PATCH 13/17] not needed (was stored as restart data) --- src/CPFEM2.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 5d9d24149..a76f018f7 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -67,7 +67,8 @@ subroutine CPFEM_initAll call homogenization_init call CPFEM_init call config_deallocate - call crystallite_setInitialValues ! ToDo: MD More general approach needed + if (interface_restartInc==0) & + call crystallite_setInitialValues ! ToDo: MD More general approach needed end subroutine CPFEM_initAll From 710c217d8a472e9496e6ffd073965a5088c341df Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Dec 2020 19:55:55 +0100 Subject: [PATCH 14/17] no extra mapping for damage --- src/constitutive_mech.f90 | 2 +- src/damage_local.f90 | 1 - src/damage_none.f90 | 1 - src/damage_nonlocal.f90 | 3 +-- src/kinematics_cleavage_opening.f90 | 2 +- src/kinematics_slipplane_opening.f90 | 2 +- src/material.f90 | 8 +------- src/prec.f90 | 6 +----- src/source_damage_anisoBrittle.f90 | 2 +- src/source_damage_anisoDuctile.f90 | 2 +- src/source_damage_isoDuctile.f90 | 2 +- 11 files changed, 9 insertions(+), 22 deletions(-) diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index d64ea327c..6b3c6fce6 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -332,7 +332,7 @@ module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el)) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2 + C = C * damage(ho)%p(material_homogenizationMemberAt(ip,el))**2 end select degradationType enddo DegradationLoop diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 0003bb2a8..97eaf9a8c 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -79,7 +79,6 @@ subroutine damage_local_init allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal) - damageMapping(h)%p => material_homogenizationMemberAt damage(h)%p => damageState(h)%state(1,:) end associate diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 1cf3de6de..3f1144833 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -30,7 +30,6 @@ subroutine damage_none_init allocate(damageState(h)%subState0(0,Nmaterialpoints)) allocate(damageState(h)%state (0,Nmaterialpoints)) - damageMapping(h)%p => material_homogenizationMemberAt allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal) enddo diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index ce97ec24a..c4426f185 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -82,7 +82,6 @@ subroutine damage_nonlocal_init allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal) - damageMapping(h)%p => material_homogenizationMemberAt damage(h)%p => damageState(h)%state(1,:) end associate @@ -179,7 +178,7 @@ subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) offset homog = material_homogenizationAt(el) - offset = damageMapping(homog)%p(ip,el) + offset = material_homogenizationMemberAt(ip,el) damage(homog)%p(offset) = phi end subroutine damage_nonlocal_putNonLocalDamage diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index 66bce6a92..6fe8ed7f6 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -120,7 +120,7 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt homog = material_homogenizationAt(el) - damageOffset = damageMapping(homog)%p(ip,el) + damageOffset = material_homogenizationMemberAt(ip,el) Ld = 0.0_pReal dLd_dTstar = 0.0_pReal diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index aa0bdfbde..b7adb6807 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -141,7 +141,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S phase = material_phaseAt(ipc,el) instance = kinematics_slipplane_opening_instance(phase) homog = material_homogenizationAt(el) - damageOffset = damageMapping(homog)%p(ip,el) + damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(instance)) Ld = 0.0_pReal diff --git a/src/material.f90 b/src/material.f90 index 33b64e6df..574da0d51 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -77,7 +77,7 @@ module material integer, dimension(:), allocatable, public, protected :: & ! (elem) material_homogenizationAt !< homogenization ID of each element - integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack + integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem) material_homogenizationMemberAt !< position of the element within its homogenization instance integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) material_phaseAt !< phase ID of each element @@ -92,9 +92,6 @@ module material type(Rotation), dimension(:,:,:), allocatable, public, protected :: & material_orientation0 !< initial orientation of each grain,IP,element - type(tHomogMapping), allocatable, dimension(:), public :: & - damageMapping !< mapping for damage state/fields - type(group_float), allocatable, dimension(:), public :: & temperature, & !< temperature field damage, & !< damage field @@ -159,11 +156,8 @@ subroutine material_init(restart) allocate(thermalState (size(material_name_homogenization))) allocate(damageState (size(material_name_homogenization))) - allocate(damageMapping (size(material_name_homogenization))) - allocate(temperature (size(material_name_homogenization))) allocate(damage (size(material_name_homogenization))) - allocate(temperatureRate (size(material_name_homogenization))) diff --git a/src/prec.f90 b/src/prec.f90 index 738775e3b..95b1116cd 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -54,7 +54,7 @@ module prec type, extends(tState) :: tPlasticState logical :: & nonlocal = .false. - real(pReal), pointer, dimension(:,:) :: & + real(pReal), pointer, dimension(:,:) :: & slipRate !< slip rate end type @@ -62,10 +62,6 @@ module prec type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase end type - type :: tHomogMapping - integer, pointer, dimension(:,:) :: p - end type - real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal) !< minimum positive number such that 1.0 + EPSILON /= 1.0. real(pReal), private, parameter :: PREAL_MIN = tiny(0.0_pReal) !< smallest normalized floating point number diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index ca8d6ec2b..55d5546fc 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -143,7 +143,7 @@ module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) constituent = material_phasememberAt(ipc,ip,el) sourceOffset = source_damage_anisoBrittle_offset(phase) homog = material_homogenizationAt(el) - damageOffset = damageMapping(homog)%p(ip,el) + damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(source_damage_anisoBrittle_instance(phase))) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 601ec2531..912fe1387 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -125,7 +125,7 @@ module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) constituent = material_phasememberAt(ipc,ip,el) sourceOffset = source_damage_anisoDuctile_offset(phase) homog = material_homogenizationAt(el) - damageOffset = damageMapping(homog)%p(ip,el) + damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(source_damage_anisoDuctile_instance(phase))) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 1bff20570..b66e220d9 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -116,7 +116,7 @@ module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) constituent = material_phasememberAt(ipc,ip,el) sourceOffset = source_damage_isoDuctile_offset(phase) homog = material_homogenizationAt(el) - damageOffset = damageMapping(homog)%p(ip,el) + damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(source_damage_isoDuctile_instance(phase))) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & From d7889aff12f6aefcd7270fb31756124fff0bbed2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Dec 2020 09:13:13 +0100 Subject: [PATCH 15/17] extra function not (yet) needed --- src/CPFEM.f90 | 3 +-- src/CPFEM2.f90 | 4 +--- src/crystallite.f90 | 23 ++++++----------------- 3 files changed, 8 insertions(+), 22 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index eb3243576..76522fb16 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -89,11 +89,10 @@ subroutine CPFEM_initAll call lattice_init call material_init(.false.) call constitutive_init - call crystallite_init call homogenization_init + call crystallite_init call CPFEM_init call config_deallocate - call crystallite_setInitialValues ! ToDo: MD More general approach needed end subroutine CPFEM_initAll diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index a76f018f7..54e381d34 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -63,12 +63,10 @@ subroutine CPFEM_initAll #endif call material_init(restart=interface_restartInc>0) call constitutive_init - call crystallite_init call homogenization_init + call crystallite_init call CPFEM_init call config_deallocate - if (interface_restartInc==0) & - call crystallite_setInitialValues ! ToDo: MD More general approach needed end subroutine CPFEM_initAll diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0fcade113..31a6bde2d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -69,9 +69,9 @@ module crystallite real(pReal), dimension(:,:,:,:,:), allocatable, public :: & crystallite_partitionedF !< def grad to be reached at end of homog inc - logical, dimension(:,:,:), allocatable, public :: & + logical, dimension(:,:,:), allocatable, public :: & crystallite_requested !< used by upper level (homogenization) to request crystallite calculation - logical, dimension(:,:,:), allocatable :: & + logical, dimension(:,:,:), allocatable :: & crystallite_converged !< convergence flag type :: tOutput !< new requested output (per phase) @@ -115,7 +115,6 @@ module crystallite public :: & crystallite_init, & - crystallite_setInitialValues, & crystallite_stress, & crystallite_stressTangent, & crystallite_orientations, & @@ -138,6 +137,9 @@ subroutine crystallite_init integer :: & p, & + c, & !< counter in integration point component loop + i, & !< counter in integration point loop + e, & !< counter in element loop cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax !< maximum number of elements @@ -255,19 +257,6 @@ subroutine crystallite_init endif #endif -end subroutine crystallite_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief Set initial values -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_setInitialValues - - integer :: & - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e !< counter in element loop - !$OMP PARALLEL DO PRIVATE(i,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Nconstituents(material_homogenizationAt(e)) @@ -306,7 +295,7 @@ subroutine crystallite_setInitialValues !$OMP END PARALLEL DO -end subroutine crystallite_setInitialValues +end subroutine crystallite_init !-------------------------------------------------------------------------------------------------- From 5d9c931008f56a5c01a4fb8b6d97a0ccee86372e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Dec 2020 11:21:24 +0100 Subject: [PATCH 16/17] code follows structure --- src/commercialFEM_fileList.f90 | 1 + src/homogenization.f90 | 214 ++++++-------------------- src/homogenization_mech.f90 | 199 ++++++++++++++++++++++++ src/homogenization_mech_RGC.f90 | 2 +- src/homogenization_mech_isostrain.f90 | 2 +- src/homogenization_mech_none.f90 | 4 +- 6 files changed, 252 insertions(+), 170 deletions(-) create mode 100644 src/homogenization_mech.f90 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index d161b36eb..beabfcae1 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -52,6 +52,7 @@ #include "damage_local.f90" #include "damage_nonlocal.f90" #include "homogenization.f90" +#include "homogenization_mech.f90" #include "homogenization_mech_none.f90" #include "homogenization_mech_isostrain.f90" #include "homogenization_mech_RGC.f90" diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 347634212..5958f35fa 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -31,13 +31,15 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point real(pReal), dimension(:,:,:,:), allocatable, public :: & - homogenization_F0, & !< def grad of IP at start of FE increment - homogenization_F !< def grad of IP to be reached at end of FE increment - real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & - homogenization_P !< first P--K stress of IP - real(pReal), dimension(:,:,:,:,:,:), allocatable, public, protected :: & - homogenization_dPdF !< tangent of first P--K stress at IP + homogenization_F0, & !< def grad of IP at start of FE increment + homogenization_F !< def grad of IP to be reached at end of FE increment + real(pReal), dimension(:,:,:,:), allocatable, public :: & !, protected :: & ! Issue with ifort + homogenization_P !< first P--K stress of IP + real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: & !, protected :: & + homogenization_dPdF !< tangent of first P--K stress at IP + +!-------------------------------------------------------------------------------------------------- type :: tNumerics integer :: & nMPstate !< materialpoint state loop limit @@ -62,52 +64,37 @@ module homogenization type(tDebugOptions) :: debugHomog + +!-------------------------------------------------------------------------------------------------- interface - module subroutine mech_none_init - end subroutine mech_none_init - - module subroutine mech_isostrain_init - end subroutine mech_isostrain_init - - module subroutine mech_RGC_init(num_homogMech) + module subroutine mech_init(num_homog) class(tNode), pointer, intent(in) :: & - num_homogMech !< pointer to mechanical homogenization numerics data - end subroutine mech_RGC_init + num_homog !< pointer to mechanical homogenization numerics data + end subroutine mech_init + module subroutine mech_partition(subF,ip,el) + real(pReal), intent(in), dimension(3,3) :: & + subF + integer, intent(in) :: & + ip, & !< integration point + el !< element number + end subroutine mech_partition - module subroutine mech_isostrain_partitionDeformation(F,avgF) - real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient - real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point - end subroutine mech_isostrain_partitionDeformation + module subroutine mech_homogenize(ip,el) + integer, intent(in) :: & + ip, & !< integration point + el !< element number + end subroutine mech_homogenize - module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) - real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient - real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point - integer, intent(in) :: & - instance, & - of - end subroutine mech_RGC_partitionDeformation + module subroutine mech_results(group_base,h) + character(len=*), intent(in) :: group_base + integer, intent(in) :: h - module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) - real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point - real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point - - real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses - real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses - integer, intent(in) :: instance - end subroutine mech_isostrain_averageStressAndItsTangent - - module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) - real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point - real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point - - real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses - real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses - integer, intent(in) :: instance - end subroutine mech_RGC_averageStressAndItsTangent + end subroutine mech_results +! -------- ToDo --------------------------------------------------------- module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) logical, dimension(2) :: mech_RGC_updateState real(pReal), dimension(:,:,:), intent(in) :: & @@ -122,13 +109,8 @@ module homogenization el !< element number end function mech_RGC_updateState - - module subroutine mech_RGC_results(instance,group) - integer, intent(in) :: instance !< homogenization instance - character(len=*), intent(in) :: group !< group name in HDF5 file - end subroutine mech_RGC_results - end interface +! ----------------------------------------------------------------------- public :: & homogenization_init, & @@ -145,10 +127,11 @@ subroutine homogenization_init class (tNode) , pointer :: & num_homog, & - num_homogMech, & num_homogGeneric, & debug_homogenization + print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) + debug_homogenization => config_debug%get('homogenization', defaultVal=emptyList) debugHomog%basic = debug_homogenization%contains('basic') debugHomog%extensive = debug_homogenization%contains('extensive') @@ -163,31 +146,8 @@ subroutine homogenization_init num_homog => config_numerics%get('homogenization',defaultVal=emptyDict) - num_homogMech => num_homog%get('mech',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) - if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init - if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init - if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init(num_homogMech) - - if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init - if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init - if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init - - if (any(damage_type == DAMAGE_none_ID)) call damage_none_init - if (any(damage_type == DAMAGE_local_ID)) call damage_local_init - if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init - - -!-------------------------------------------------------------------------------------------------- -! allocate and initialize global variables - allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) - homogenization_F0 = spread(spread(math_I3,3,discretization_nIPs),4,discretization_Nelems) ! initialize to identity - homogenization_F = homogenization_F0 ! initialize to identity - allocate(homogenization_P(3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) - - print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) - num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10) num%subStepMinHomog = num_homogGeneric%get_asFloat('subStepMin', defaultVal=1.0e-3_pReal) num%subStepSizeHomog = num_homogGeneric%get_asFloat('subStepSize', defaultVal=0.25_pReal) @@ -198,6 +158,18 @@ subroutine homogenization_init if (num%subStepSizeHomog <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeHomog') if (num%stepIncreaseHomog <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseHomog') + + call mech_init(num_homog) + + if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init + if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init + if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init + + if (any(damage_type == DAMAGE_none_ID)) call damage_none_init + if (any(damage_type == DAMAGE_local_ID)) call damage_local_init + if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init + + end subroutine homogenization_init @@ -330,7 +302,7 @@ subroutine materialpoint_stressAndItsTangent(dt) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e)) IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2) if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done - call partitionDeformation(homogenization_F0(1:3,1:3,i,e) & + call mech_partition(homogenization_F0(1:3,1:3,i,e) & + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e))& *(subStep(i,e)+subFrac(i,e)), & i,e) @@ -379,7 +351,7 @@ subroutine materialpoint_stressAndItsTangent(dt) !$OMP PARALLEL DO elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) IpLooping4: do i = FEsolving_execIP(1),FEsolving_execIP(2) - call averageStressAndItsTangent(i,e) + call mech_homogenize(i,e) enddo IpLooping4 enddo elementLooping4 !$OMP END PARALLEL DO @@ -390,38 +362,6 @@ subroutine materialpoint_stressAndItsTangent(dt) end subroutine materialpoint_stressAndItsTangent -!-------------------------------------------------------------------------------------------------- -!> @brief partition material point def grad onto constituents -!-------------------------------------------------------------------------------------------------- -subroutine partitionDeformation(subF,ip,el) - - real(pReal), intent(in), dimension(3,3) :: & - subF - integer, intent(in) :: & - ip, & !< integration point - el !< element number - - chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) - - case (HOMOGENIZATION_NONE_ID) chosenHomogenization - crystallite_partitionedF(1:3,1:3,1,ip,el) = subF - - case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - call mech_isostrain_partitionDeformation(& - crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - subF) - - case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call mech_RGC_partitionDeformation(& - crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - subF,& - ip, & - el) - end select chosenHomogenization - -end subroutine partitionDeformation - - !-------------------------------------------------------------------------------------------------- !> @brief update the internal state of the homogenization scheme and tell whether "done" and !> "happy" with result @@ -478,49 +418,6 @@ function updateState(subdt,subF,ip,el) end function updateState -!-------------------------------------------------------------------------------------------------- -!> @brief derive average stress and stiffness from constituent quantities -!-------------------------------------------------------------------------------------------------- -subroutine averageStressAndItsTangent(ip,el) - - integer, intent(in) :: & - ip, & !< integration point - el !< element number - integer :: c - real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) - - - chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) - case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) - homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el) - - case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) - enddo - call mech_isostrain_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,ip,el), & - homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& - crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - dPdFs, & - homogenization_typeInstance(material_homogenizationAt(el))) - - case (HOMOGENIZATION_RGC_ID) chosenHomogenization - do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) - enddo - call mech_RGC_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,ip,el), & - homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& - crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - dPdFs, & - homogenization_typeInstance(material_homogenizationAt(el))) - end select chosenHomogenization - -end subroutine averageStressAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief writes homogenization results to HDF5 output file !-------------------------------------------------------------------------------------------------- @@ -531,27 +428,12 @@ subroutine homogenization_results integer :: p character(len=:), allocatable :: group_base,group - !real(pReal), dimension(:,:,:), allocatable :: temp do p=1,size(material_name_homogenization) group_base = 'current/homogenization/'//trim(material_name_homogenization(p)) call results_closeGroup(results_addGroup(group_base)) - group = trim(group_base)//'/generic' - call results_closeGroup(results_addGroup(group)) - !temp = reshape(homogenization_F,[3,3,discretization_nIPs*discretization_Nelems]) - !call results_writeDataset(group,temp,'F',& - ! 'deformation gradient','1') - !temp = reshape(homogenization_P,[3,3,discretization_nIPs*discretization_Nelems]) - !call results_writeDataset(group,temp,'P',& - ! '1st Piola-Kirchhoff stress','Pa') - - group = trim(group_base)//'/mech' - call results_closeGroup(results_addGroup(group)) - select case(material_homogenization_type(p)) - case(HOMOGENIZATION_rgc_ID) - call mech_RGC_results(homogenization_typeInstance(p),group) - end select + call mech_results(group_base,p) group = trim(group_base)//'/damage' call results_closeGroup(results_addGroup(group)) diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 new file mode 100644 index 000000000..40d26df51 --- /dev/null +++ b/src/homogenization_mech.f90 @@ -0,0 +1,199 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief Partition F and homogenize P/dPdF +!-------------------------------------------------------------------------------------------------- +submodule(homogenization) homogenization_mech + + interface + + module subroutine mech_none_init + end subroutine mech_none_init + + module subroutine mech_isostrain_init + end subroutine mech_isostrain_init + + module subroutine mech_RGC_init(num_homogMech) + class(tNode), pointer, intent(in) :: & + num_homogMech !< pointer to mechanical homogenization numerics data + end subroutine mech_RGC_init + + + module subroutine mech_isostrain_partitionDeformation(F,avgF) + real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient + real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point + end subroutine mech_isostrain_partitionDeformation + + module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) + real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient + real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point + integer, intent(in) :: & + instance, & + of + end subroutine mech_RGC_partitionDeformation + + + module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) + real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point + real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point + + real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses + real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + integer, intent(in) :: instance + end subroutine mech_isostrain_averageStressAndItsTangent + + module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) + real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point + real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point + + real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses + real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + integer, intent(in) :: instance + end subroutine mech_RGC_averageStressAndItsTangent + + + module subroutine mech_RGC_results(instance,group) + integer, intent(in) :: instance !< homogenization instance + character(len=*), intent(in) :: group !< group name in HDF5 file + end subroutine mech_RGC_results + + end interface + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief Allocate variables and set parameters. +!-------------------------------------------------------------------------------------------------- +module subroutine mech_init(num_homog) + + class(tNode), pointer, intent(in) :: & + num_homog + + class(tNode), pointer :: & + num_homogMech + + print'(/,a)', ' <<<+- homogenization_mech init -+>>>' + + allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) + homogenization_F0 = spread(spread(math_I3,3,discretization_nIPs),4,discretization_Nelems) ! initialize to identity + homogenization_F = homogenization_F0 ! initialize to identity + allocate(homogenization_P(3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) + + num_homogMech => num_homog%get('mech',defaultVal=emptyDict) + if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init + if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init + if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init(num_homogMech) + +end subroutine mech_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief Partition F onto the individual constituents. +!-------------------------------------------------------------------------------------------------- +module subroutine mech_partition(subF,ip,el) + + real(pReal), intent(in), dimension(3,3) :: & + subF + integer, intent(in) :: & + ip, & !< integration point + el !< element number + + chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) + + case (HOMOGENIZATION_NONE_ID) chosenHomogenization + crystallite_partitionedF(1:3,1:3,1,ip,el) = subF + + case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization + call mech_isostrain_partitionDeformation(& + crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & + subF) + + case (HOMOGENIZATION_RGC_ID) chosenHomogenization + call mech_RGC_partitionDeformation(& + crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & + subF,& + ip, & + el) + + end select chosenHomogenization + +end subroutine mech_partition + + +!-------------------------------------------------------------------------------------------------- +!> @brief Average P and dPdF from the individual constituents. +!-------------------------------------------------------------------------------------------------- +module subroutine mech_homogenize(ip,el) + + integer, intent(in) :: & + ip, & !< integration point + el !< element number + integer :: c + real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + + + chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) + + case (HOMOGENIZATION_NONE_ID) chosenHomogenization + homogenization_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el) + + case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization + do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) + enddo + call mech_isostrain_averageStressAndItsTangent(& + homogenization_P(1:3,1:3,ip,el), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& + crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & + dPdFs, & + homogenization_typeInstance(material_homogenizationAt(el))) + + case (HOMOGENIZATION_RGC_ID) chosenHomogenization + do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) + enddo + call mech_RGC_averageStressAndItsTangent(& + homogenization_P(1:3,1:3,ip,el), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& + crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & + dPdFs, & + homogenization_typeInstance(material_homogenizationAt(el))) + + end select chosenHomogenization + +end subroutine mech_homogenize + + +!-------------------------------------------------------------------------------------------------- +!> @brief Write results to file. +!-------------------------------------------------------------------------------------------------- +module subroutine mech_results(group_base,h) + use material, only: & + material_homogenization_type => homogenization_type + + character(len=*), intent(in) :: group_base + integer, intent(in) :: h + + character(len=:), allocatable :: group + + group = trim(group_base)//'/mech' + call results_closeGroup(results_addGroup(group)) + + select case(material_homogenization_type(h)) + + case(HOMOGENIZATION_rgc_ID) + call mech_RGC_results(homogenization_typeInstance(h),group) + + end select + + !temp = reshape(homogenization_F,[3,3,discretization_nIPs*discretization_Nelems]) + !call results_writeDataset(group,temp,'F',& + ! 'deformation gradient','1') + !temp = reshape(homogenization_P,[3,3,discretization_nIPs*discretization_Nelems]) + !call results_writeDataset(group,temp,'P',& + ! '1st Piola-Kirchhoff stress','Pa') + +end subroutine mech_results + + +end submodule homogenization_mech diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 585752469..0a9d0ac92 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -6,7 +6,7 @@ !> @brief Relaxed grain cluster (RGC) homogenization scheme !> N_constituents is defined as p x q x r (cluster) !-------------------------------------------------------------------------------------------------- -submodule(homogenization) homogenization_mech_RGC +submodule(homogenization:homogenization_mech) homogenization_mech_RGC use rotations type :: tParameters diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index 751518e09..a56104647 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Isostrain (full constraint Taylor assuption) homogenization scheme !-------------------------------------------------------------------------------------------------- -submodule(homogenization) homogenization_mech_isostrain +submodule(homogenization:homogenization_mech) homogenization_mech_isostrain enum, bind(c); enumerator :: & parallel_ID, & diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90 index 5b12247cd..d434d1ca0 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mech_none.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief dummy homogenization homogenization scheme for 1 constituent per material point !-------------------------------------------------------------------------------------------------- -submodule(homogenization) homogenization_mech_none +submodule(homogenization:homogenization_mech) homogenization_mech_none contains @@ -28,7 +28,7 @@ module subroutine mech_none_init if(homogenization_Nconstituents(h) /= 1) & call IO_error(211,ext_msg='N_constituents (mech_none)') - + Nmaterialpoints = count(material_homogenizationAt == h) homogState(h)%sizeState = 0 allocate(homogState(h)%state0 (0,Nmaterialpoints)) From 3884549e194c5f5e25e71c7d510d659132c7035e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Dec 2020 12:48:45 +0100 Subject: [PATCH 17/17] store field variables as 1D array first step of simplifying layout: 1) Solver translates from ip,el tuple (FEM) or cells(1),cells(2),cells(3) triple to list. 2) DAMASK iterates over all points 3) homogenization knows mapping (point,constituent) -> (instance,member) --- src/CPFEM.f90 | 18 ++++---- src/grid/grid_mech_FEM.f90 | 30 +++++++------- src/grid/grid_mech_spectral_basic.f90 | 4 +- src/grid/grid_mech_spectral_polarisation.f90 | 6 +-- src/grid/spectral_utilities.f90 | 16 ++++---- src/homogenization.f90 | 23 ++++++----- src/homogenization_mech.f90 | 23 ++++++----- src/material.f90 | 1 - src/mesh/FEM_utilities.f90 | 2 +- src/mesh/mesh_mech_FEM.f90 | 43 ++++++++++---------- 10 files changed, 86 insertions(+), 80 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 76522fb16..a19a70432 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -153,7 +153,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS H integer(pInt) elCP, & ! crystal plasticity element number - i, j, k, l, m, n, ph, homog, mySource + i, j, k, l, m, n, ph, homog, mySource,ma real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll @@ -161,6 +161,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS elCP = mesh_FEM2DAMASK_elem(elFE) + ma = (elCP-1) * discretization_nIPs + ip + if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then print'(/,a)', '#############################################' print'(a1,a22,1x,i8,a13)', '#','element', elCP, '#' @@ -184,8 +186,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & temperature_inp end select chosenThermal1 - homogenization_F0(1:3,1:3,ip,elCP) = ffn - homogenization_F(1:3,1:3,ip,elCP) = ffn1 + homogenization_F0(1:3,1:3,ma) = ffn + homogenization_F(1:3,1:3,ma) = ffn1 if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then @@ -212,17 +214,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS else terminalIllness ! translate from P to sigma - Kirchhoff = matmul(homogenization_P(1:3,1:3,ip,elCP), transpose(homogenization_F(1:3,1:3,ip,elCP))) - J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ip,elCP)) + Kirchhoff = matmul(homogenization_P(1:3,1:3,ma), transpose(homogenization_F(1:3,1:3,ma))) + J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ma)) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) ! translate from dP/dF to dCS/dE H = 0.0_pReal do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 H(i,j,k,l) = H(i,j,k,l) & - + homogenization_F(j,m,ip,elCP) * homogenization_F(l,n,ip,elCP) & - * homogenization_dPdF(i,m,k,n,ip,elCP) & - - math_delta(j,l) * homogenization_F(i,m,ip,elCP) * homogenization_P(k,m,ip,elCP) & + + homogenization_F(j,m,ma) * homogenization_F(l,n,ma) & + * homogenization_dPdF(i,m,k,n,ma) & + - math_delta(j,l) * homogenization_F(i,m,ma) * homogenization_P(k,m,ma) & + 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) & + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k)) enddo; enddo; enddo; enddo; enddo; enddo diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 741ce404a..cdf806b35 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -238,7 +238,7 @@ subroutine grid_mech_FEM_init F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent 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 @@ -359,7 +359,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,& F_lastInc = F - homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F, [3,3,product(grid(1:2))*grid3]) endif !-------------------------------------------------------------------------------------------------- @@ -557,9 +557,9 @@ subroutine formResidual(da_local,x_local, & ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 ele = ele + 1 f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & - matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,1,ele) + & - homogenization_dPdF(2,2,2,2,1,ele) + & - homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal + 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 ctr = 0 do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 ctr = ctr + 1 @@ -636,18 +636,18 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr) row = col ele = ele + 1 K_ele = 0.0 - K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + & - homogenization_dPdF(2,2,2,2,1,ele) + & - homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + & - homogenization_dPdF(2,2,2,2,1,ele) + & - homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + & - homogenization_dPdF(2,2,2,2,1,ele) + & - homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,ele) + & + homogenization_dPdF(3,3,3,3,ele))/3.0_pReal + K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,ele) + & + homogenization_dPdF(3,3,3,3,ele))/3.0_pReal + K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,ele) + & + homogenization_dPdF(3,3,3,3,ele))/3.0_pReal K_ele = K_ele + & matmul(transpose(BMatFull), & - matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,1,ele), & + matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,ele), & shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) CHKERRQ(ierr) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index d87a22fb7..ebaaf3b55 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -199,7 +199,7 @@ subroutine grid_mech_spectral_basic_init F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent 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 @@ -319,7 +319,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_ rotation_BC%rotate(F_aimDot,active=.true.)) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) - homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) endif !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 9bbb40a53..9f2a17c97 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -225,7 +225,7 @@ subroutine grid_mech_spectral_polarisation_init F_tau_lastInc = 2.0_pReal*F_lastInc endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent 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 @@ -359,7 +359,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) - homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) endif !-------------------------------------------------------------------------------------------------- @@ -604,7 +604,7 @@ subroutine formResidual(in, FandF_tau, & do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) e = e + 1 residual_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & + math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + residual_F_tau(1:3,1:3,i,j,k) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 989448dc3..c0c84233d 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -810,7 +810,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& print'(/,a)', ' ... evaluating constitutive response ......................................' flush(IO_STDOUT) - homogenization_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field call materialpoint_stressAndItsTangent(timeinc) ! calculate P field @@ -829,13 +829,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& dPdF_min = huge(1.0_pReal) dPdF_norm_min = huge(1.0_pReal) do i = 1, product(grid(1:2))*grid3 - if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then - dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i) - dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) + if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) 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.0_pReal) endif - if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then - dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i) - dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) + if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) then + dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i) + dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal) endif end do @@ -853,7 +853,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min) - C_volAvg = sum(sum(homogenization_dPdF,dim=6),dim=5) + C_volAvg = sum(homogenization_dPdF,dim=5) call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) if (ierr /= 0) error stop 'MPI error' C_volAvg = C_volAvg * wgt diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 5958f35fa..57478e039 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -30,12 +30,12 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point - real(pReal), dimension(:,:,:,:), allocatable, public :: & + real(pReal), dimension(:,:,:), allocatable, public :: & homogenization_F0, & !< def grad of IP at start of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment - real(pReal), dimension(:,:,:,:), allocatable, public :: & !, protected :: & ! Issue with ifort + real(pReal), dimension(:,:,:), allocatable, public :: & !, protected :: & Issue with ifort homogenization_P !< first P--K stress of IP - real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: & !, protected :: & + real(pReal), dimension(:,:,:,:,:), allocatable, public :: & !, protected :: & homogenization_dPdF !< tangent of first P--K stress at IP @@ -193,6 +193,7 @@ subroutine materialpoint_stressAndItsTangent(dt) converged logical, dimension(2,discretization_nIPs,discretization_Nelems) :: & doneAndHappy + integer :: m !-------------------------------------------------------------------------------------------------- @@ -227,7 +228,7 @@ subroutine materialpoint_stressAndItsTangent(dt) any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),& FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) - !$OMP PARALLEL DO + !$OMP PARALLEL DO PRIVATE(m) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e)) IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -297,13 +298,14 @@ subroutine materialpoint_stressAndItsTangent(dt) !-------------------------------------------------------------------------------------------------- ! deformation partitioning - !$OMP PARALLEL DO PRIVATE(myNgrains) + !$OMP PARALLEL DO PRIVATE(myNgrains,m) elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e)) IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2) if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done - call mech_partition(homogenization_F0(1:3,1:3,i,e) & - + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e))& + m = (e-1)*discretization_nIPs + i + call mech_partition(homogenization_F0(1:3,1:3,m) & + + (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m))& *(subStep(i,e)+subFrac(i,e)), & i,e) crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains @@ -321,16 +323,17 @@ subroutine materialpoint_stressAndItsTangent(dt) !-------------------------------------------------------------------------------------------------- ! state update - !$OMP PARALLEL DO + !$OMP PARALLEL DO PRIVATE(m) elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2) if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then if (.not. converged(i,e)) then doneAndHappy(1:2,i,e) = [.true.,.false.] else + m = (e-1)*discretization_nIPs + i doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), & - homogenization_F0(1:3,1:3,i,e) & - + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e)) & + homogenization_F0(1:3,1:3,m) & + + (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m)) & *(subStep(i,e)+subFrac(i,e)), & i,e) converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 40d26df51..b0641be07 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -73,10 +73,10 @@ module subroutine mech_init(num_homog) print'(/,a)', ' <<<+- homogenization_mech init -+>>>' - allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) - homogenization_F0 = spread(spread(math_I3,3,discretization_nIPs),4,discretization_Nelems) ! initialize to identity - homogenization_F = homogenization_F0 ! initialize to identity - allocate(homogenization_P(3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) + allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) + homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity + homogenization_F = homogenization_F0 ! initialize to identity + allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) num_homogMech => num_homog%get('mech',defaultVal=emptyDict) if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init @@ -127,23 +127,24 @@ module subroutine mech_homogenize(ip,el) integer, intent(in) :: & ip, & !< integration point el !< element number - integer :: c + integer :: c,m real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + m = (el-1)* discretization_nIPs + ip chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) - homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el) + homogenization_P(1:3,1:3,m) = crystallite_P(1:3,1:3,1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,m) = crystallite_stressTangent(1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,ip,el), & - homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& + homogenization_P(1:3,1:3,m), & + homogenization_dPdF(1:3,1:3,1:3,1:3,m),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) @@ -153,8 +154,8 @@ module subroutine mech_homogenize(ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,ip,el), & - homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& + homogenization_P(1:3,1:3,m), & + homogenization_dPdF(1:3,1:3,1:3,1:3,m),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) diff --git a/src/material.f90 b/src/material.f90 index 574da0d51..bb5f484f6 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -140,7 +140,6 @@ contains subroutine material_init(restart) logical, intent(in) :: restart - integer :: myHomog print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 118735e89..cb81f1f0c 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -164,7 +164,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) cutBack = .false. ! reset cutBack status - P_av = sum(sum(homogenization_P,dim=4),dim=3) * wgt ! average of P + P_av = sum(homogenization_P,dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) end subroutine utilities_constitutiveResponse diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index eb5f862c2..e19c35998 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -316,16 +316,16 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) Vec :: x_local, f_local, xx_local PetscSection :: section PetscScalar, dimension(:), pointer :: x_scal, pf_scal - PetscScalar, target :: f_scal(cellDof) - PetscReal :: detJ, IcellJMat(dimPlex,dimPlex) - PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer + PetscScalar, dimension(cellDof), target :: f_scal + PetscReal :: IcellJMat(dimPlex,dimPlex) + PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer PetscInt :: cellStart, cellEnd, cell, field, face, & qPt, basis, comp, cidx, & - numFields - PetscReal :: detFAvg - PetscReal :: BMat(dimPlex*dimPlex,cellDof) + numFields, & + bcSize,m + PetscReal :: detFAvg, detJ + PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat - PetscInt :: bcSize IS :: bcPoints @@ -366,6 +366,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) CHKERRQ(ierr) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) do qPt = 0, nQuadrature-1 + m = cell*nQuadrature + qPt+1 BMat = 0.0 do basis = 0, nBasis-1 do comp = 0, dimPlex-1 @@ -375,15 +376,14 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) enddo enddo - homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & - reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) + homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) enddo if (num%BBarStabilisation) then - detFAvg = math_det33(sum(homogenization_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) - do qPt = 1, nQuadrature - homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1) = & - homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1)* & - (detFAvg/math_det33(homogenization_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex)) + detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature)) + do qPt = 0, nQuadrature-1 + m = cell*nQuadrature + qPt+1 + homogenization_F(1:dimPlex,1:dimPlex,m) = homogenization_F(1:dimPlex,1:dimPlex,m) & + * (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0/real(dimPlex)) enddo endif @@ -407,6 +407,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) f_scal = 0.0 do qPt = 0, nQuadrature-1 + m = cell*nQuadrature + qPt+1 BMat = 0.0 do basis = 0, nBasis-1 do comp = 0, dimPlex-1 @@ -418,7 +419,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) enddo f_scal = f_scal + & matmul(transpose(BMat), & - reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), & + reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), & shape=[dimPlex*dimPlex]))*qWeights(qPt+1) enddo f_scal = f_scal*abs(detJ) @@ -463,7 +464,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) K_eB PetscInt :: cellStart, cellEnd, cell, field, face, & - qPt, basis, comp, cidx,bcSize + qPt, basis, comp, cidx,bcSize, m IS :: bcPoints @@ -506,6 +507,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) FAvg = 0.0 BMatAvg = 0.0 do qPt = 0, nQuadrature-1 + m = cell*nQuadrature + qPt + 1 BMat = 0.0 do basis = 0, nBasis-1 do comp = 0, dimPlex-1 @@ -516,7 +518,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) enddo enddo - MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & + MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) if (num%BBarStabilisation) then @@ -524,12 +526,11 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) FInv = math_inv33(F) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) K_eB = K_eB - & - matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), & - shape=[dimPlex*dimPlex,1]), & + matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex*dimPlex,1]), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), & shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA) - MatB = MatB + & - matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA) + MatB = MatB & + + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1,dimPlex*dimPlex]),MatA) FAvg = FAvg + F BMatAvg = BMatAvg + BMat else