From da9fdf53d20e9781e980e32f4750f8a1b805f2da Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 15 Nov 2021 12:35:44 -0500 Subject: [PATCH 1/5] consistent indentation and line-spacings in reporting --- src/CPFEM.f90 | 2 +- src/CPFEM2.f90 | 2 +- src/DAMASK_Marc.f90 | 2 +- src/DAMASK_interface.f90 | 2 +- src/HDF5_utilities.f90 | 2 +- src/IO.f90 | 2 +- src/YAML_types.f90 | 2 +- src/base64.f90 | 2 +- src/config.f90 | 25 ++-- src/discretization.f90 | 2 +- src/element.f90 | 2 +- src/grid/DAMASK_grid.f90 | 52 +++---- src/grid/discretization_grid.f90 | 80 +++++----- src/grid/grid_damage_spectral.f90 | 28 ++-- src/grid/grid_mech_FEM.f90 | 24 +-- src/grid/grid_mech_spectral_basic.f90 | 60 ++++---- src/grid/grid_mech_spectral_polarisation.f90 | 74 ++++----- src/grid/grid_thermal_spectral.f90 | 28 ++-- src/grid/spectral_utilities.f90 | 108 +++++++------- src/homogenization.f90 | 2 +- src/homogenization_damage.f90 | 2 +- src/homogenization_damage_pass.f90 | 2 +- src/homogenization_mechanical.f90 | 12 +- src/homogenization_mechanical_RGC.f90 | 90 +++++------ src/homogenization_mechanical_isostrain.f90 | 6 +- src/homogenization_mechanical_pass.f90 | 10 +- src/homogenization_thermal.f90 | 16 +- src/homogenization_thermal_pass.f90 | 2 +- src/lattice.f90 | 2 +- src/material.f90 | 30 ++-- src/math.f90 | 8 +- src/mesh/DAMASK_mesh.f90 | 88 +++++------ src/mesh/FEM_quadrature.f90 | 6 +- src/mesh/FEM_utilities.f90 | 10 +- src/mesh/discretization_mesh.f90 | 17 ++- src/mesh/mesh_mech_FEM.f90 | 26 ++-- src/parallelization.f90 | 10 +- src/phase.f90 | 52 +++---- src/phase_damage.f90 | 38 ++--- src/phase_damage_anisobrittle.f90 | 26 ++-- src/phase_damage_isobrittle.f90 | 16 +- src/phase_mechanical.f90 | 2 +- src/phase_mechanical_eigen.f90 | 20 +-- ...phase_mechanical_eigen_cleavageopening.f90 | 4 +- ...hase_mechanical_eigen_thermalexpansion.f90 | 24 +-- src/phase_mechanical_elastic.f90 | 14 +- src/phase_mechanical_plastic.f90 | 12 +- ...phase_mechanical_plastic_dislotungsten.f90 | 14 +- src/phase_mechanical_plastic_dislotwin.f90 | 72 ++++----- src/phase_mechanical_plastic_isotropic.f90 | 20 +-- ...phase_mechanical_plastic_kinehardening.f90 | 22 +-- src/phase_mechanical_plastic_none.f90 | 10 +- src/phase_mechanical_plastic_nonlocal.f90 | 140 +++++++++--------- ...phase_mechanical_plastic_phenopowerlaw.f90 | 30 ++-- src/phase_thermal.f90 | 2 +- src/phase_thermal_dissipation.f90 | 14 +- src/phase_thermal_externalheat.f90 | 16 +- src/prec.f90 | 16 +- src/results.f90 | 42 +++--- src/rotations.f90 | 6 +- 60 files changed, 725 insertions(+), 725 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 4123af37e..8f1217a88 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -103,7 +103,7 @@ subroutine CPFEM_init class(tNode), pointer :: & debug_CPFEM - print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- CPFEM init -+>>>'; flush(IO_STDOUT) allocate(CPFEM_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal) allocate(CPFEM_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal) diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 2bb8420b9..647befd30 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -81,7 +81,7 @@ subroutine CPFEM_init integer(HID_T) :: fileHandle - print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- CPFEM init -+>>>'; flush(IO_STDOUT) if (interface_restartInc > 0) then diff --git a/src/DAMASK_Marc.f90 b/src/DAMASK_Marc.f90 index e4859b4ab..0525db323 100644 --- a/src/DAMASK_Marc.f90 +++ b/src/DAMASK_Marc.f90 @@ -46,7 +46,7 @@ subroutine DAMASK_interface_init integer :: ierr character(len=pPathLen) :: wd - print'(/,a)', ' <<<+- DAMASK_marc init -+>>>' + print'(/,1x,a)', '<<<+- DAMASK_marc init -+>>>' print*, 'Roters et al., Computational Materials Science 158:420–478, 2019' print*, 'https://doi.org/10.1016/j.commatsci.2018.04.030' diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index c0d149457..f5233f2f0 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -70,7 +70,7 @@ subroutine DAMASK_interface_init external :: & quit - print'(/,a)', ' <<<+- DAMASK_interface init -+>>>' + print'(/,1x,a)', '<<<+- DAMASK_interface init -+>>>' if(worldrank == 0) open(OUTPUT_UNIT, encoding='UTF-8') ! for special characters in output diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 540ce6781..5e0e40646 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -109,7 +109,7 @@ subroutine HDF5_utilities_init integer(SIZE_T) :: typeSize - print'(/,a)', ' <<<+- HDF5_Utilities init -+>>>' + print'(/,1x,a)', '<<<+- HDF5_Utilities init -+>>>' call h5open_f(hdferr) diff --git a/src/IO.f90 b/src/IO.f90 index 717493006..ec7cf9704 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -56,7 +56,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine IO_init - print'(/,a)', ' <<<+- IO init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- IO init -+>>>'; flush(IO_STDOUT) call selfTest diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index 8f81923a8..127f380cd 100644 --- a/src/YAML_types.f90 +++ b/src/YAML_types.f90 @@ -187,7 +187,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine YAML_types_init - print'(/,a)', ' <<<+- YAML_types init -+>>>' + print'(/,1x,a)', '<<<+- YAML_types init -+>>>' call selfTest diff --git a/src/base64.f90 b/src/base64.f90 index 2f91334b7..7bc00a1f9 100644 --- a/src/base64.f90 +++ b/src/base64.f90 @@ -27,7 +27,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine base64_init - print'(/,a)', ' <<<+- base64 init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- base64 init -+>>>'; flush(IO_STDOUT) call selfTest diff --git a/src/config.f90 b/src/config.f90 index 0dfb8c878..aaf6fee1a 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -30,8 +30,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine config_init - print'(/,a)', ' <<<+- config init -+>>>'; flush(IO_STDOUT) - + print'(/,1x,a)', '<<<+- config init -+>>>'; flush(IO_STDOUT) call parse_material call parse_numerics @@ -50,15 +49,15 @@ subroutine parse_material() inquire(file='material.yaml',exist=fileExists) - if(.not. fileExists) call IO_error(100,ext_msg='material.yaml') + if (.not. fileExists) call IO_error(100,ext_msg='material.yaml') if (worldrank == 0) then - print*, 'reading material.yaml'; flush(IO_STDOUT) + print'(/,1x,a)', 'reading material.yaml'; flush(IO_STDOUT) fileContent = IO_read('material.yaml') call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup','material.yaml','main configuration') call results_closeJobFile - endif + end if call parallelization_bcast_str(fileContent) config_material => YAML_parse_str(fileContent) @@ -81,19 +80,19 @@ subroutine parse_numerics() if (fileExists) then if (worldrank == 0) then - print*, 'reading numerics.yaml'; flush(IO_STDOUT) + print'(1x,a)', 'reading numerics.yaml'; flush(IO_STDOUT) fileContent = IO_read('numerics.yaml') if (len(fileContent) > 0) then call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration') call results_closeJobFile - endif - endif + end if + end if call parallelization_bcast_str(fileContent) config_numerics => YAML_parse_str(fileContent) - endif + end if end subroutine parse_numerics @@ -113,19 +112,19 @@ subroutine parse_debug() if (fileExists) then if (worldrank == 0) then - print*, 'reading debug.yaml'; flush(IO_STDOUT) + print'(1x,a)', 'reading debug.yaml'; flush(IO_STDOUT) fileContent = IO_read('debug.yaml') if (len(fileContent) > 0) then call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration') call results_closeJobFile - endif - endif + end if + end if call parallelization_bcast_str(fileContent) config_debug => YAML_parse_str(fileContent) - endif + end if end subroutine parse_debug diff --git a/src/discretization.f90 b/src/discretization.f90 index 17077d5aa..c79c72b5e 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -49,7 +49,7 @@ subroutine discretization_init(materialAt,& integer, optional, intent(in) :: & sharedNodesBegin !< index of first node shared among different processes (MPI) - print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6) + print'(/,1x,a)', '<<<+- discretization init -+>>>'; flush(6) discretization_Nelems = size(materialAt,1) discretization_nIPs = size(IPcoords0,2)/discretization_Nelems diff --git a/src/element.f90 b/src/element.f90 index b061f942d..c8e3b0416 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -923,7 +923,7 @@ subroutine tElement_init(self,elemType) self%nIPneighbors = size(self%IPneighbor,1) - print'(/,a)', ' <<<+- element_init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- element_init -+>>>'; flush(IO_STDOUT) print*, 'element type: ',self%elemType print*, ' geom type: ',self%geomType diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 7ffa6b80a..4e79bb6a8 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -113,10 +113,10 @@ program DAMASK_grid ! init DAMASK (all modules) call CPFEM_initAll - print'(/,a)', ' <<<+- DAMASK_grid init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- DAMASK_grid init -+>>>'; flush(IO_STDOUT) - print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' - print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' + print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' + print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80' !------------------------------------------------------------------------------------------------- @@ -227,12 +227,12 @@ program DAMASK_grid loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1) reportAndCheck: if (worldrank == 0) then - print'(/,a,i0)', ' load case: ', l - print*, ' estimate_rate:', loadCases(l)%estimate_rate + print'(/,1x,a,1x,i0)', 'load case:', l + print'(2x,a,1x,l1)', 'estimate_rate:', loadCases(l)%estimate_rate if (loadCases(l)%deformation%myType == 'F') then - print*, ' F:' + print'(2x,a)', 'F:' else - print*, ' '//loadCases(l)%deformation%myType//' / 1/s:' + print'(2x,a)', loadCases(l)%deformation%myType//' / 1/s:' endif do i = 1, 3; do j = 1, 3 if (loadCases(l)%deformation%mask(i,j)) then @@ -246,8 +246,8 @@ program DAMASK_grid if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) & errorID = 838 ! no rotation is allowed by stress BC - if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:' - if (loadCases(l)%stress%myType == 'dot_P') print*, ' dot_P / MPa/s:' + if (loadCases(l)%stress%myType == 'P') print'(2x,a)', 'P / MPa:' + if (loadCases(l)%stress%myType == 'dot_P') print'(2x,a)', 'dot_P / MPa/s:' if (loadCases(l)%stress%myType /= '') then do i = 1, 3; do j = 1, 3 @@ -270,15 +270,15 @@ program DAMASK_grid if (loadCases(l)%f_restart < 1) errorID = 839 if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then - print'(a)', ' r: 1 (constant step width)' + print'(2x,a)', 'r: 1 (constant step width)' else - print'(a,f0.3)', ' r: ', loadCases(l)%r + print'(2x,a,1x,f0.3)', 'r:', loadCases(l)%r endif - print'(a,f0.3)', ' t: ', loadCases(l)%t - print'(a,i0)', ' N: ', loadCases(l)%N - print'(a,i0)', ' f_out: ', loadCases(l)%f_out + print'(2x,a,1x,f0.3)', 't:', loadCases(l)%t + print'(2x,a,1x,i0)', 'N:', loadCases(l)%N + print'(2x,a,1x,i0)', 'f_out:', loadCases(l)%f_out if (loadCases(l)%f_restart < huge(0)) & - print'(a,i0)', ' f_restart: ', loadCases(l)%f_restart + print'(2x,a,1x,i0)', 'f_restart:', loadCases(l)%f_restart if (errorID > 0) call IO_error(error_ID = errorID, el = l) @@ -317,7 +317,7 @@ program DAMASK_grid endif writeUndeformed: if (interface_restartInc < 1) then - print'(/,a)', ' ... writing initial configuration to file ........................' + print'(/,1x,a)', '... writing initial configuration to file ........................' flush(IO_STDOUT) call CPFEM_results(0,0.0_pReal) endif writeUndeformed @@ -353,8 +353,8 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! report begin of new step - print'(/,a)', ' ###########################################################################' - print'(1x,a,es12.5,6(a,i0))', & + print'(/,1x,a)', '###########################################################################' + print'(1x,a,1x,es12.5,6(a,i0))', & 'Time', t, & 's: Increment ', inc,'/',loadCases(l)%N,& '-', stepFraction,'/',subStepFactor**cutBackLevel,& @@ -379,7 +379,7 @@ program DAMASK_grid case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) end select enddo - if(.not. cutBack) call CPFEM_forward + if (.not. cutBack) call CPFEM_forward !-------------------------------------------------------------------------------------------------- ! solve fields @@ -425,7 +425,7 @@ program DAMASK_grid cutBackLevel = cutBackLevel + 1 t = t - Delta_t Delta_t = Delta_t/real(subStepFactor,pReal) ! cut timestep - print'(/,a)', ' cutting back ' + print'(/,1x,a)', 'cutting back ' else ! no more options to continue if (worldrank == 0) close(statUnit) call IO_error(950) @@ -436,26 +436,26 @@ program DAMASK_grid cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc if (all(solres(:)%converged)) then - print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged' + print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' converged' else - print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged' + print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' NOT converged' endif; flush(IO_STDOUT) call MPI_Allreduce(interface_SIGUSR1,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) if (ierr /= 0) error stop 'MPI error' if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then - print'(1/,a)', ' ... writing results to file ......................................' + print'(/,1x,a)', '... writing results to file ......................................' flush(IO_STDOUT) call CPFEM_results(totalIncsCounter,t) endif - if(signal) call interface_setSIGUSR1(.false.) + if (signal) call interface_setSIGUSR1(.false.) call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) if (ierr /= 0) error stop 'MPI error' if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then call mechanical_restartWrite call CPFEM_restartWrite endif - if(signal) call interface_setSIGUSR2(.false.) + if (signal) call interface_setSIGUSR2(.false.) call MPI_Allreduce(interface_SIGTERM,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) if (ierr /= 0) error stop 'MPI error' if (signal) exit loadCaseLooping @@ -468,7 +468,7 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation - print'(/,a)', ' ###########################################################################' + print'(/,1x,a)', '###########################################################################' if (worldrank == 0) close(statUnit) call quit(0) ! no complains ;) diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index e16001feb..871be2c74 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -72,10 +72,10 @@ subroutine discretization_grid_init(restart) fileContent, fname - print'(/,a)', ' <<<+- discretization_grid init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- discretization_grid init -+>>>'; flush(IO_STDOUT) - if(worldrank == 0) then + if (worldrank == 0) then fileContent = IO_read(interface_geomFile) call readVTI(grid,geomSize,origin,materialAt_global,fileContent) fname = interface_geomFile @@ -85,7 +85,7 @@ subroutine discretization_grid_init(restart) call results_closeJobFile else allocate(materialAt_global(0)) ! needed for IntelMPI - endif + end if call MPI_Bcast(grid,3,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) @@ -96,11 +96,11 @@ subroutine discretization_grid_init(restart) call MPI_Bcast(origin,3,MPI_DOUBLE,0,MPI_COMM_WORLD, ierr) if (ierr /= 0) error stop 'MPI error' - print'(/,a,3(i12 ))', ' cells a b c: ', grid - print'(a,3(es12.5))', ' size x y z: ', geomSize - print'(a,3(es12.5))', ' origin x y z: ', origin + print'(/,1x,a,3(i12,1x))', 'cells a b c: ', grid + print '(1x,a,3(es12.5,1x))', 'size x y z: ', geomSize + print '(1x,a,3(es12.5,1x))', 'origin x y z: ', origin - if(worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)') + if (worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)') call fftw_mpi_init devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & @@ -109,7 +109,7 @@ subroutine discretization_grid_init(restart) PETSC_COMM_WORLD, & z, & ! domain grid size along z z_offset) ! domain grid offset along z - if(z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes') + if (z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes') grid3 = int(z) grid3Offset = int(z_offset) @@ -136,14 +136,14 @@ subroutine discretization_grid_init(restart) !-------------------------------------------------------------------------------------------------- ! store geometry information for post processing - if(.not. restart) then + if (.not. restart) then call results_openJobFile call results_closeGroup(results_addGroup('geometry')) call results_addAttribute('cells', grid, '/geometry') call results_addAttribute('size', geomSize,'/geometry') call results_addAttribute('origin',origin, '/geometry') call results_closeJobFile - endif + end if !-------------------------------------------------------------------------------------------------- ! geometry information required by the nonlocal CP model @@ -202,18 +202,18 @@ subroutine readVTI(grid,geomSize,origin,material, & if (endPos < startPos) endPos = len(fileContent,kind=pI64) ! end of file without new line if (.not. inFile) then - if(index(fileContent(startPos:endPos),'',kind=pI64) /= 0_pI64) then gotCellData = .true. @@ -230,25 +230,25 @@ subroutine readVTI(grid,geomSize,origin,material, & s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64 ! start (no leading whitespace) material = as_Int(fileContent(s:endPos),headerType,compressed,dataType) exit - endif + end if startPos = endPos + 2_pI64 endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64 - enddo - endif - endif - endif + end do + end if + end if + end if if (gotCellData) exit startPos = endPos + 2_pI64 - enddo + end do - if(.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found') - if(size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)') - if(any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size') - if(any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid') + if (.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found') + if (size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)') + if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size') + if (any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid') material = material + 1 - if(any(material<1)) call IO_error(error_ID = 844, ext_msg='material ID < 0') + if (any(material<1)) call IO_error(error_ID = 844, ext_msg='material ID < 0') contains @@ -352,11 +352,11 @@ subroutine readVTI(grid,geomSize,origin,material, & integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes - if(compressed) then + if (compressed) then bytes = asBytes_compressed(base64_str,headerType) else bytes = asBytes_uncompressed(base64_str,headerType) - endif + end if end function asBytes @@ -384,12 +384,12 @@ subroutine readVTI(grid,geomSize,origin,material, & nBlock = int(temp(1),pI64) headerLen = 4_pI64 * (3_pI64 + nBlock) temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64) - elseif(headerType == 'UInt64') then + elseif (headerType == 'UInt64') then temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(8_pI64)))),pI64) nBlock = int(temp(1),pI64) headerLen = 8_pI64 * (3_pI64 + nBlock) temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64) - endif + end if allocate(size_inflated(nBlock),source=temp(2)) size_inflated(nBlock) = merge(temp(3),temp(2),temp(3)/=0_pI64) @@ -402,7 +402,7 @@ subroutine readVTI(grid,geomSize,origin,material, & s = e + 1_pI64 e = s + size_deflated(b) - 1_pI64 bytes(sum(size_inflated(:b-1))+1_pI64:sum(size_inflated(:b))) = zlib_inflate(bytes_inflated(s:e),size_inflated(b)) - enddo + end do end function asBytes_compressed @@ -429,14 +429,14 @@ subroutine readVTI(grid,geomSize,origin,material, & nByte = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64)))),pI64) bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64+nByte(1))),5_pI64)] s = s + base64_nChar(4_pI64+nByte(1)) - enddo - elseif(headerType == 'UInt64') then + end do + elseif (headerType == 'UInt64') then do while(s+base64_nChar(8_pI64)<(len(base64_str,pI64))) nByte = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64)))),pI64) bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64+nByte(1))),9_pI64)] s = s + base64_nChar(8_pI64+nByte(1)) - enddo - endif + end do + end if end function asBytes_uncompressed @@ -455,11 +455,11 @@ subroutine readVTI(grid,geomSize,origin,material, & #endif s = index(line," "//key,back=.true.) - if(s==0) then + if (s==0) then getXMLValue = '' else e = s + 1 + scan(line(s+1:),"'"//'"') - if(scan(line(s:e-2),'=') == 0) then + if (scan(line(s:e-2),'=') == 0) then getXMLValue = '' else s = e @@ -471,8 +471,8 @@ subroutine readVTI(grid,geomSize,origin,material, & e = s + index(line(s:),merge("'",'"',line(s-1:s-1)=="'")) - 1 #endif getXMLValue = line(s:e-1) - endif - endif + end if + end if end function @@ -514,7 +514,7 @@ function IPcoordinates0(grid,geomSize,grid3Offset) do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1) i = i + 1 IPcoordinates0(1:3,i) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal) - enddo; enddo; enddo + end do; end do; end do end function IPcoordinates0 @@ -538,7 +538,7 @@ pure function nodes0(grid,geomSize,grid3Offset) do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1) n = n + 1 nodes0(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal) - enddo; enddo; enddo + end do; end do; end do end function nodes0 @@ -631,7 +631,7 @@ pure function IPneighborhood(grid) IPneighborhood(3,5,1,e) = 6 IPneighborhood(3,6,1,e) = 5 - enddo; enddo; enddo + end do; end do; end do end function IPneighborhood diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index f43a6cb18..187a6f552 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -76,10 +76,10 @@ subroutine grid_damage_spectral_init() character(len=pStringLen) :: & snes_type - print'(/,a)', ' <<<+- grid_spectral_damage init -+>>>' + print'(/,1x,a)', '<<<+- grid_spectral_damage init -+>>>' - print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' - print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' + print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' + print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80' !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks @@ -137,7 +137,7 @@ subroutine grid_damage_spectral_init() call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) ! variable bounds for variational inequalities like contact mechanics, damage etc. call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) - endif + end if !-------------------------------------------------------------------------------------------------- ! init fields @@ -187,7 +187,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution) else solution%converged = .true. solution%iterationsNeeded = totalIter - endif + end if stagNorm = maxval(abs(phi_current - phi_stagInc)) solnNorm = maxval(abs(phi_current)) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr) @@ -201,14 +201,14 @@ function grid_damage_spectral_solution(Delta_t) result(solution) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 call homogenization_set_phi(phi_current(i,j,k),ce) - enddo; enddo; enddo + end do; end do; end do call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr) if (solution%converged) & - print'(/,a)', ' ... nonlocal damage converged .....................................' - print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm - print'(/,a)', ' ===========================================================================' + print'(/,1x,a)', '... nonlocal damage converged .....................................' + print'(/,1x,a,f8.6,2x,f8.6,2x,e11.4)', 'Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm + print'(/,1x,a)', '===========================================================================' flush(IO_STDOUT) end function grid_damage_spectral_solution @@ -238,11 +238,11 @@ subroutine grid_damage_spectral_forward(cutBack) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 call homogenization_set_phi(phi_current(i,j,k),ce) - enddo; enddo; enddo + end do; end do; end do else phi_lastInc = phi_current call updateReference - endif + end if end subroutine grid_damage_spectral_forward @@ -277,7 +277,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k)) - enddo; enddo; enddo + end do; end do; end do call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field call utilities_FFTscalarBackward @@ -287,7 +287,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + mu_ref*phi_current(i,j,k) - enddo; enddo; enddo + end do; end do; end do !-------------------------------------------------------------------------------------------------- ! convolution of damage field with green operator @@ -320,7 +320,7 @@ subroutine updateReference() do ce = 1, product(grid(1:2))*grid3 K_ref = K_ref + homogenization_K_phi(ce) mu_ref = mu_ref + homogenization_mu_phi(ce) - enddo + end do K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index ca616dffb..9befd7ade 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -116,7 +116,7 @@ subroutine grid_mechanical_FEM_init num_grid, & debug_grid - print'(/,a)', ' <<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT) !------------------------------------------------------------------------------------------------- ! debugging options @@ -234,7 +234,7 @@ subroutine grid_mechanical_FEM_init !-------------------------------------------------------------------------------------------------- ! init fields restartRead: if (interface_restartInc > 0) then - print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file' + print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file' fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r') groupHandle = HDF5_openGroup(fileHandle,'solver') @@ -272,7 +272,7 @@ subroutine grid_mechanical_FEM_init CHKERRQ(ierr) restartRead2: if (interface_restartInc > 0) then - print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file' + print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file' call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.) call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) if(ierr /=0) error stop 'MPI error' @@ -442,7 +442,7 @@ subroutine grid_mechanical_FEM_restartWrite call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr) CHKERRQ(ierr) - print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) + print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT) fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w') groupHandle = HDF5_addGroup(fileHandle,'solver') @@ -506,12 +506,12 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i reason = 0 endif - print'(1/,a)', ' ... reporting .............................................................' - print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & + print'(/,1x,a)', '... reporting .............................................................' + print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', & err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' - print'(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', & + print'(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error stress BC = ', & err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - print'(/,a)', ' ===========================================================================' + print'(/,1x,a)', '===========================================================================' flush(IO_STDOUT) end subroutine converged @@ -547,10 +547,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) 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) + if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + 'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) + print'(/,1x,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 b033a7b29..66b703ed6 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -111,13 +111,13 @@ subroutine grid_mechanical_spectral_basic_init num_grid, & debug_grid - print'(/,a)', ' <<<+- grid_mechanical_spectral_basic init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- grid_mechanical_spectral_basic init -+>>>'; flush(IO_STDOUT) - print*, 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' - print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL + print'(/,1x,a)', 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' + print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL - print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' - print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' + print'( 1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' + print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006' !------------------------------------------------------------------------------------------------- ! debugging options @@ -186,30 +186,30 @@ subroutine grid_mechanical_spectral_basic_init call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! places pointer on PETSc data restartRead: if (interface_restartInc > 0) then - print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file' + print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file' fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r') groupHandle = HDF5_openGroup(fileHandle,'solver') call HDF5_read(P_aim,groupHandle,'P_aim',.false.) call MPI_Bcast(P_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(F_aim,groupHandle,'F_aim',.false.) call MPI_Bcast(F_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.) call MPI_Bcast(F_aim_lastInc,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.) call MPI_Bcast(F_aimDot,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(F,groupHandle,'F') call HDF5_read(F_lastInc,groupHandle,'F_lastInc') elseif (interface_restartInc == 0) then restartRead F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) - endif restartRead + end if restartRead homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) @@ -219,13 +219,13 @@ subroutine grid_mechanical_spectral_basic_init call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer restartRead2: if (interface_restartInc > 0) then - print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file' + print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file' call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.) call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.) call MPI_Bcast(C_volAvgLastInc,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) @@ -234,7 +234,7 @@ subroutine grid_mechanical_spectral_basic_init MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,ierr) call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr) call MPI_File_close(fileUnit,ierr) - endif restartRead2 + end if restartRead2 call utilities_updateGamma(C_minMaxAvg) call utilities_saveReferenceStiffness @@ -263,7 +263,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) - if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) + if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg) !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -328,7 +328,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed F_aimDot = F_aimDot & + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask) - endif + end if Fdot = utilities_calculateRate(guess, & F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, & @@ -336,7 +336,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) - endif + end if !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients @@ -385,7 +385,7 @@ subroutine grid_mechanical_spectral_basic_restartWrite call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) - print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) + print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT) fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w') groupHandle = HDF5_addGroup(fileHandle,'solver') @@ -406,7 +406,7 @@ subroutine grid_mechanical_spectral_basic_restartWrite call HDF5_write(C_minMaxAvg,groupHandle,'C_minMaxAvg',.false.) call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) - endif + end if if (num%update_gamma) call utilities_saveReferenceStiffness @@ -443,14 +443,14 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm reason = -1 else reason = 0 - endif + end if - print'(1/,a)', ' ... reporting .............................................................' - print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & + print'(/,1x,a)', '... reporting .............................................................' + print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', & err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' - print'(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', & + print'(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error stress BC = ', & err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - print'(/,a)', ' ===========================================================================' + print'(/,1x,a)', '===========================================================================' flush(IO_STDOUT) end subroutine converged @@ -485,12 +485,12 @@ 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) 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) + if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + 'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) + print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + 'deformation gradient aim =', transpose(F_aim) flush(IO_STDOUT) - endif newIteration + end if newIteration !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 8258ad43d..3d524029c 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -124,10 +124,10 @@ subroutine grid_mechanical_spectral_polarisation_init num_grid, & debug_grid - print'(/,a)', ' <<<+- grid_mechanical_spectral_polarization init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- grid_mechanical_spectral_polarization init -+>>>'; flush(IO_STDOUT) - print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' - print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' + print'(/,1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' + print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006' !------------------------------------------------------------------------------------------------- ! debugging options @@ -208,23 +208,23 @@ subroutine grid_mechanical_spectral_polarisation_init F_tau => FandF_tau(9:17,:,:,:) restartRead: if (interface_restartInc > 0) then - print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file' + print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file' fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r') groupHandle = HDF5_openGroup(fileHandle,'solver') call HDF5_read(P_aim,groupHandle,'P_aim',.false.) call MPI_Bcast(P_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(F_aim,groupHandle,'F_aim',.false.) call MPI_Bcast(F_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.) call MPI_Bcast(F_aim_lastInc,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.) call MPI_Bcast(F_aimDot,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(F,groupHandle,'F') call HDF5_read(F_lastInc,groupHandle,'F_lastInc') call HDF5_read(F_tau,groupHandle,'F_tau') @@ -235,7 +235,7 @@ subroutine grid_mechanical_spectral_polarisation_init F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) F_tau = 2.0_pReal*F F_tau_lastInc = 2.0_pReal*F_lastInc - endif restartRead + end if restartRead homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) @@ -245,13 +245,13 @@ subroutine grid_mechanical_spectral_polarisation_init call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer restartRead2: if (interface_restartInc > 0) then - print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file' + print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file' call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.) call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.) call MPI_Bcast(C_volAvgLastInc,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) @@ -260,7 +260,7 @@ subroutine grid_mechanical_spectral_polarisation_init MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,ierr) call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr) call MPI_File_close(fileUnit,ierr) - endif restartRead2 + end if restartRead2 call utilities_updateGamma(C_minMaxAvg) call utilities_saveReferenceStiffness @@ -291,11 +291,11 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) - if(num%update_gamma) then + if (num%update_gamma) then call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) - endif + end if !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -364,7 +364,7 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed F_aimDot = F_aimDot & + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask) - endif + end if Fdot = utilities_calculateRate(guess, & F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, & @@ -376,14 +376,14 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) - endif + end if !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * Delta_t - if(stress_BC%myType=='P') P_aim = P_aim & + if (stress_BC%myType=='P') P_aim = P_aim & + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t - if(stress_BC%myType=='dot_P') P_aim = P_aim & + if (stress_BC%myType=='dot_P') P_aim = P_aim & + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average @@ -399,8 +399,8 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D + math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, & math_mul3333xx33(C_scale,matmul(transpose(F_lambda33),F_lambda33)-math_I3))) F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) - enddo; enddo; enddo - endif + end do; end do; end do + end if call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) @@ -441,7 +441,7 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite F => FandF_tau(0: 8,:,:,:) F_tau => FandF_tau(9:17,:,:,:) - print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) + print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT) fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w') groupHandle = HDF5_addGroup(fileHandle,'solver') @@ -463,9 +463,9 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite call HDF5_write(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.) call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) - endif + end if - if(num%update_gamma) call utilities_saveReferenceStiffness + if (num%update_gamma) call utilities_saveReferenceStiffness call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) @@ -502,16 +502,16 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm reason = -1 else reason = 0 - endif + end if - print'(1/,a)', ' ... reporting .............................................................' - print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & + print'(/,1x,a)', '... reporting .............................................................' + print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', & err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' - print '(a,f12.2,a,es8.2,a,es9.2,a)', ' error curl = ', & + print '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error curl = ', & err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' - print '(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', & + print '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error stress BC = ', & err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - print'(/,a)', ' ===========================================================================' + print'(/,1x,a)', '===========================================================================' flush(IO_STDOUT) end subroutine converged @@ -565,12 +565,12 @@ 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) 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) + if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + 'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) + print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & + 'deformation gradient aim =', transpose(F_aim) flush(IO_STDOUT) - endif newIteration + end if newIteration !-------------------------------------------------------------------------------------------------- ! @@ -580,7 +580,7 @@ subroutine formResidual(in, FandF_tau, & num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& num%alpha*matmul(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) - enddo; enddo; enddo + end do; end do; end do !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space @@ -621,7 +621,7 @@ subroutine formResidual(in, FandF_tau, & 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) - enddo; enddo; enddo + end do; end do; end do !-------------------------------------------------------------------------------------------------- ! calculating curl diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 41fca2133..868d3101e 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -75,10 +75,10 @@ subroutine grid_thermal_spectral_init(T_0) class(tNode), pointer :: & num_grid - print'(/,a)', ' <<<+- grid_thermal_spectral init -+>>>' + print'(/,1x,a)', '<<<+- grid_thermal_spectral init -+>>>' - print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' - print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' + print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' + print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80' !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks @@ -141,7 +141,7 @@ subroutine grid_thermal_spectral_init(T_0) T_lastInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k) call homogenization_thermal_setField(T_0,0.0_pReal,ce) - enddo; enddo; enddo + end do; end do; end do call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,ierr); CHKERRQ(ierr) T_PETSc(xstart:xend,ystart:yend,zstart:zend) = T_current @@ -182,7 +182,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) else solution%converged = .true. solution%iterationsNeeded = totalIter - endif + end if stagNorm = maxval(abs(T_current - T_stagInc)) solnNorm = maxval(abs(T_current)) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr) @@ -196,14 +196,14 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) - enddo; enddo; enddo + end do; end do; end do call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr) if (solution%converged) & - print'(/,a)', ' ... thermal conduction converged ..................................' - print'(/,a,f8.4,2x,f8.4,2x,f8.4)', ' Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm - print'(/,a)', ' ===========================================================================' + print'(/,1x,a)', '... thermal conduction converged ..................................' + print'(/,1x,a,f8.4,2x,f8.4,2x,f8.4)', 'Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm + print'(/,1x,a)', '===========================================================================' flush(IO_STDOUT) end function grid_thermal_spectral_solution @@ -234,11 +234,11 @@ subroutine grid_thermal_spectral_forward(cutBack) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) - enddo; enddo; enddo + end do; end do; end do else T_lastInc = T_current call updateReference - endif + end if end subroutine grid_thermal_spectral_forward @@ -272,7 +272,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k)) - enddo; enddo; enddo + end do; end do; end do call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field call utilities_FFTscalarBackward @@ -282,7 +282,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + mu_ref*T_current(i,j,k) - enddo; enddo; enddo + end do; end do; end do !-------------------------------------------------------------------------------------------------- ! convolution of temperature field with green operator @@ -310,7 +310,7 @@ subroutine updateReference() do ce = 1, product(grid(1:2))*grid3 K_ref = K_ref + homogenization_K_T(ce) mu_ref = mu_ref + homogenization_mu_T(ce) - enddo + end do K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 579397104..da8c3f05e 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -177,19 +177,19 @@ subroutine spectral_utilities_init num_grid, & debug_grid ! pointer to grid debug options - print'(/,a)', ' <<<+- spectral_utilities init -+>>>' + print'(/,1x,a)', '<<<+- spectral_utilities init -+>>>' - print*, 'M. Diehl, Diploma Thesis TU München, 2010' - print*, 'https://doi.org/10.13140/2.1.3234.3840'//IO_EOL + print'(/,1x,a)', 'M. Diehl, Diploma Thesis TU München, 2010' + print'( 1x,a)', 'https://doi.org/10.13140/2.1.3234.3840'//IO_EOL - print*, 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' - print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL + print'( 1x,a)', 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' + print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL - print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' - print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'//IO_EOL + print'( 1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' + print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'//IO_EOL - print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' - print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' + print'( 1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' + print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80' !-------------------------------------------------------------------------------------------------- ! set debugging parameters @@ -200,15 +200,15 @@ subroutine spectral_utilities_init debugRotation = debug_grid%contains('rotation') debugPETSc = debug_grid%contains('PETSc') - if(debugPETSc) print'(3(/,a),/)', & - ' Initializing PETSc with debug options: ', & + if (debugPETSc) print'(3(/,1x,a),/)', & + 'Initializing PETSc with debug options: ', & trim(PETScDebug), & - ' add more using the "PETSc_options" keyword in numerics.yaml' + 'add more using the "PETSc_options" keyword in numerics.yaml' flush(IO_STDOUT) call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr) CHKERRQ(ierr) - if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) + if (debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,& num_grid%get_asString('PETSc_options',defaultVal=''),ierr) @@ -271,7 +271,7 @@ subroutine spectral_utilities_init if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match' call fftw_set_timelimit(num_grid%get_asFloat('fftw_timelimit',defaultVal=-1.0_pReal)) - print*, 'FFTW initialized'; flush(IO_STDOUT) + print'(/,1x,a)', 'FFTW initialized'; flush(IO_STDOUT) !-------------------------------------------------------------------------------------------------- ! MPI allocation @@ -342,10 +342,10 @@ subroutine spectral_utilities_init ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) do k = grid3Offset+1, grid3Offset+grid3 k_s(3) = k - 1 - if(k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + if (k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 do j = 1, grid(2) k_s(2) = j - 1 - if(j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + if (j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 do i = 1, grid1Red k_s(1) = i - 1 ! symmetry, junst running from 0,1,...,N/2,N/2+1 xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s) @@ -357,7 +357,7 @@ subroutine spectral_utilities_init endwhere enddo; enddo; enddo - if(num%memory_efficient) then ! allocate just single fourth order tensor + if (num%memory_efficient) then ! allocate just single fourth order tensor allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal)) else ! precalculation of gamma_hat field allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal)) @@ -384,7 +384,7 @@ subroutine utilities_updateGamma(C) C_ref = C - if(.not. num%memory_efficient) then + if (.not. num%memory_efficient) then gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 @@ -497,12 +497,12 @@ subroutine utilities_fourierGammaConvolution(fieldAim) logical :: err - print'(/,a)', ' ... doing gamma convolution ...............................................' + print'(/,1x,a)', '... doing gamma convolution ...............................................' flush(IO_STDOUT) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation (mechanical equilibrium) - memoryEfficient: if(num%memory_efficient) then + memoryEfficient: if (num%memory_efficient) then do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 forall(l = 1:3, m = 1:3) & @@ -567,7 +567,7 @@ real(pReal) function utilities_divergenceRMS() integer :: i, j, k, ierr complex(pReal), dimension(3) :: rescaledGeom - print'(/,a)', ' ... calculating divergence ................................................' + print'(/,1x,a)', '... calculating divergence ................................................' flush(IO_STDOUT) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) @@ -593,9 +593,9 @@ real(pReal) function utilities_divergenceRMS() + sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), & conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal) enddo; enddo - if(grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + if (grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space end function utilities_divergenceRMS @@ -610,7 +610,7 @@ real(pReal) function utilities_curlRMS() complex(pReal), dimension(3,3) :: curl_fourier complex(pReal), dimension(3) :: rescaledGeom - print'(/,a)', ' ... calculating curl ......................................................' + print'(/,1x,a)', '... calculating curl ......................................................' flush(IO_STDOUT) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) @@ -655,9 +655,9 @@ real(pReal) function utilities_curlRMS() enddo; enddo call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' utilities_curlRMS = sqrt(utilities_curlRMS) * wgt - if(grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + if (grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 end function utilities_curlRMS @@ -686,13 +686,13 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) mask_stressVector = .not. reshape(transpose(mask_stress), [9]) size_reduced = count(mask_stressVector) - if(size_reduced > 0) then + if (size_reduced > 0) then temp99_real = math_3333to99(rot_BC%rotate(C)) - if(debugGeneral) then - print'(/,a)', ' ... updating masked compliance ............................................' - print'(/,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', & - ' Stiffness C (load) / GPa =', transpose(temp99_Real)*1.0e-9_pReal + if (debugGeneral) then + print'(/,1x,a)', '... updating masked compliance ............................................' + print'(/,1x,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 @@ -711,10 +711,10 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pReal)) if (debugGeneral .or. errmatinv) then write(formatString, '(i2)') size_reduced - formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' - 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' + formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' + 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]) else @@ -723,9 +723,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) utilities_maskedCompliance = math_99to3333(temp99_Real) - if(debugGeneral) then - print'(/,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', & - ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal + if (debugGeneral) then + print'(/,1x,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 @@ -810,7 +810,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& real(pReal) :: dPdF_norm_max, dPdF_norm_min real(pReal), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF - print'(/,a)', ' ... evaluating constitutive response ......................................' + print'(/,1x,a)', '... evaluating constitutive response ......................................' flush(IO_STDOUT) homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field @@ -824,11 +824,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 call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) - 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 + if (debugRotation) print'(/,1x,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'(/,1x,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 @@ -1017,7 +1017,7 @@ subroutine utilities_updateCoords(F) call utilities_FFTtensorForward() do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red - if(any([i,j,k+grid3Offset] /= 1)) then + if (any([i,j,k+grid3Offset] /= 1)) then vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) & / sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal) else @@ -1031,7 +1031,7 @@ subroutine utilities_updateCoords(F) ! average F if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt call MPI_Bcast(Favg,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' !-------------------------------------------------------------------------------------------------- ! pad cell center fluctuations along z-direction (needed when running MPI simulation) @@ -1042,22 +1042,22 @@ subroutine utilities_updateCoords(F) ! send bottom layer to process below call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0,MPI_COMM_WORLD,request(1),ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0,MPI_COMM_WORLD,request(2),ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' ! send top layer to process above call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,1,MPI_COMM_WORLD,request(3),ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,1,MPI_COMM_WORLD,request(4),ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' call MPI_Waitall(4,request,status,ierr) - if(ierr /=0) error stop 'MPI error' + if (ierr /=0) error stop 'MPI error' #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) ! ToDo #else - if(any(status(MPI_ERROR,:) /= 0)) error stop 'MPI error' + if (any(status(MPI_ERROR,:) /= 0)) error stop 'MPI error' #endif !-------------------------------------------------------------------------------------------------- @@ -1094,10 +1094,10 @@ subroutine utilities_saveReferenceStiffness fileUnit,ierr if (worldrank == 0) then - print'(a)', ' writing reference stiffness data required for restart to file'; flush(IO_STDOUT) + print'(1x,a)', 'writing reference stiffness data required for restart to file'; flush(IO_STDOUT) open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',& status='replace',access='stream',action='write',iostat=ierr) - if(ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref') + if (ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref') write(fileUnit) C_ref close(fileUnit) endif diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 924e45989..cc29ed619 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -199,7 +199,7 @@ subroutine homogenization_init() num_homog, & num_homogGeneric - print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- homogenization init -+>>>'; flush(IO_STDOUT) allocate(homogState (size(material_name_homogenization))) diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 251cc72dc..6ae4ac07a 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -41,7 +41,7 @@ module subroutine damage_init() integer :: ho,Nmembers - print'(/,a)', ' <<<+- homogenization:damage init -+>>>' + print'(/,1x,a)', '<<<+- homogenization:damage init -+>>>' configHomogenizations => config_material%get('homogenization') diff --git a/src/homogenization_damage_pass.f90 b/src/homogenization_damage_pass.f90 index bdb44b822..41573ff49 100644 --- a/src/homogenization_damage_pass.f90 +++ b/src/homogenization_damage_pass.f90 @@ -8,7 +8,7 @@ contains module subroutine pass_init() - print'(/,a)', ' <<<+- homogenization:damage:pass init -+>>>' + print'(/,1x,a)', '<<<+- homogenization:damage:pass init -+>>>' end subroutine pass_init diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index 0138d2468..d2ec08394 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -68,7 +68,7 @@ module subroutine mechanical_init(num_homog) class(tNode), pointer :: & num_homogMech - print'(/,a)', ' <<<+- homogenization:mechanical init -+>>>' + print'(/,1x,a)', '<<<+- homogenization:mechanical init -+>>>' call material_parseHomogenization2() @@ -114,7 +114,7 @@ module subroutine mechanical_partition(subF,ce) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) call phase_set_F(Fs(1:3,1:3,co),co,ce) - enddo + end do end subroutine mechanical_partition @@ -138,7 +138,7 @@ module subroutine mechanical_homogenize(Delta_t,ce) + phase_P(co,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) & + phase_mechanical_dPdF(Delta_t,co,ce) - enddo + end do homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) & / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) @@ -173,11 +173,11 @@ module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce) Fs(:,:,co) = phase_F(co,ce) Ps(:,:,co) = phase_P(co,ce) - enddo + end do doneAndHappy = RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) else doneAndHappy = .true. - endif + end if end function mechanical_updateState @@ -241,7 +241,7 @@ subroutine material_parseHomogenization2() case default call IO_error(500,ext_msg=homogMech%get_asString('type')) end select - enddo + end do end subroutine material_parseHomogenization2 diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index cd10aa5e2..b02e48498 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -87,16 +87,16 @@ module subroutine RGC_init(num_homogMech) homog, & homogMech - print'(/,a)', ' <<<+- homogenization:mechanical:RGC init -+>>>' + print'(/,1x,a)', '<<<+- homogenization:mechanical:RGC init -+>>>' - print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID) + print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID) flush(IO_STDOUT) - print*, 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' - print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL + print'(/,1x,a)', 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' + print'( 1x,a)', 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL - print*, 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010' - print*, 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL + print'(/,1x,a)', 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010' + print'( 1x,a)', 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL material_homogenization => config_material%get('homogenization') @@ -186,7 +186,7 @@ module subroutine RGC_init(num_homogMech) end associate - enddo + end do end subroutine RGC_init @@ -222,9 +222,9 @@ module subroutine RGC_partitionDeformation(F,avgF,ce) nVect = interfaceNormal(intFace,ho,en) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation - enddo + end do F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient - enddo + end do end associate @@ -257,10 +257,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax - zeroTimeStep: if(dEq0(dt)) then + zeroTimeStep: if (dEq0(dt)) then doneAndHappy = .true. ! pretend everything is fine and return return - endif zeroTimeStep + end if zeroTimeStep ho = material_homogenizationID(ce) en = material_homogenizationEntry(ce) @@ -319,10 +319,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) & ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface + (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j) resid(i+3*(iNum-1)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array - enddo - enddo + end do + end do - enddo + end do !-------------------------------------------------------------------------------------------------- ! convergence check for stress residual @@ -347,7 +347,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then ! try to restart when residual blows up exceeding maximum bound doneAndHappy = [.true.,.false.] ! with direct cut-back return - endif + end if !--------------------------------------------------------------------------------------------------- ! construct the global Jacobian matrix for updating the global relaxation vector array when @@ -373,11 +373,11 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) do i=1,3; do j=1,3; do k=1,3; do l=1,3 smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l) - enddo;enddo;enddo;enddo + end do;enddo;enddo;enddo ! projecting the material tangent dPdF into the interface ! to obtain the Jacobian matrix contribution of dPdF - endif - enddo + end if + end do !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) @@ -394,10 +394,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) do i=1,3; do j=1,3; do k=1,3; do l=1,3 smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l) - enddo;enddo;enddo;enddo - endif - enddo - enddo + end do;enddo;enddo;enddo + end if + end do + end do !-------------------------------------------------------------------------------------------------- ! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical @@ -443,10 +443,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) & + (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) & + (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j) - enddo; enddo - enddo + end do; end do + end do pmatrix(:,ipert) = p_resid/num%pPert - enddo + end do !-------------------------------------------------------------------------------------------------- @@ -455,7 +455,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) do i=1,3*nIntFaceTot rmatrix(i,i) = num%viscModus*num%viscPower/(num%refRelaxRate*dt)* & ! tangent due to numerical viscosity traction appears (abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal) ! only in the main diagonal term - enddo + end do !-------------------------------------------------------------------------------------------------- @@ -472,7 +472,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) drelax = 0.0_pReal do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable - enddo; enddo + end do; end do stt%relaxationVector(:,en) = relax + drelax ! Updateing the state variable for the next iteration if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large doneAndHappy = [.true.,.false.] @@ -481,7 +481,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax)) flush(IO_STDOUT) !$OMP END CRITICAL (write2out) - endif + end if end associate @@ -545,9 +545,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) do i = 1,3; do j = 1,3 do k = 1,3; do l = 1,3 nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient - enddo; enddo + end do; end do nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor - enddo; enddo + end do; end do nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity) nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces) @@ -560,11 +560,11 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) *cosh(prm%c_alpha*nDefNorm) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) & *tanh(nDefNorm/num%xSmoo) - enddo; enddo;enddo; enddo - enddo interfaceLoop + end do; end do;enddo; end do + end do interfaceLoop - enddo grainLoop + end do grainLoop end associate @@ -594,7 +594,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) gVol(i) = math_det33(fDef(1:3,1:3,i)) ! compute the volume of individual grains vDiscrep = vDiscrep - gVol(i)/real(nGrain,pReal) ! calculate the difference/dicrepancy between ! the volume of the cluster and the the total volume of grains - enddo + end do !---------------------------------------------------------------------------------------------- ! calculate the stress and penalty due to volume discrepancy @@ -603,7 +603,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* & sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* & gVol(i)*transpose(math_inv33(fDef(:,:,i))) - enddo + end do end subroutine volumePenalty @@ -633,9 +633,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) nVect = interfaceNormal([iBase,1,1,1],ho,en) do i = 1,3; do j = 1,3 surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal - enddo; enddo + end do; end do surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF ! get the surface correction factor (area contraction/enlargement) - enddo + end do end function surfaceCorrection @@ -690,9 +690,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) nVect = interfaceNormal(intFace,ho,en) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations - enddo + end do F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient - enddo + end do end associate @@ -727,7 +727,7 @@ module subroutine RGC_results(ho,group) call results_writeDataset(dst%relaxationrate_avg,group,trim(prm%output(o)), & 'average relaxation rate','m/s') end select - enddo outputsLoop + end do outputsLoop end associate end subroutine RGC_results @@ -756,7 +756,7 @@ pure function relaxationVector(intFace,ho,en) relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en) else relaxationVector = 0.0_pReal - endif + end if end associate @@ -855,7 +855,7 @@ integer pure function interface4to1(iFace4D, nGDim) else interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) & + nGDim(2)*nGDim(3)*(iFace4D(2)-1) - endif + end if case(2) if ((iFace4D(3) == 0) .or. (iFace4D(3) == nGDim(2))) then @@ -864,7 +864,7 @@ integer pure function interface4to1(iFace4D, nGDim) interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) & + nGDim(3)*nGDim(1)*(iFace4D(3)-1) & + (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total # of interfaces normal || e1 - endif + end if case(3) if ((iFace4D(4) == 0) .or. (iFace4D(4) == nGDim(3))) then @@ -874,7 +874,7 @@ integer pure function interface4to1(iFace4D, nGDim) + nGDim(1)*nGDim(2)*(iFace4D(4)-1) & + (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total # of interfaces normal || e1 + nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total # of interfaces normal || e2 - endif + end if case default interface4to1 = -1 @@ -918,7 +918,7 @@ pure function interface1to4(iFace1D, nGDim) interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1 interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1 interface1to4(4) = int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)/real(nGDim(2),pReal))+1 - endif + end if end function interface1to4 diff --git a/src/homogenization_mechanical_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90 index b35aea6e4..81c2cd3b9 100644 --- a/src/homogenization_mechanical_isostrain.f90 +++ b/src/homogenization_mechanical_isostrain.f90 @@ -17,9 +17,9 @@ module subroutine isostrain_init ho, & Nmembers - print'(/,a)', ' <<<+- homogenization:mechanical:isostrain init -+>>>' + print'(/,1x,a)', '<<<+- homogenization:mechanical:isostrain init -+>>>' - print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) + print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) flush(IO_STDOUT) do ho = 1, size(homogenization_type) @@ -30,7 +30,7 @@ module subroutine isostrain_init allocate(homogState(ho)%state0(0,Nmembers)) allocate(homogState(ho)%state (0,Nmembers)) - enddo + end do end subroutine isostrain_init diff --git a/src/homogenization_mechanical_pass.f90 b/src/homogenization_mechanical_pass.f90 index dda50c409..1bcdf817a 100644 --- a/src/homogenization_mechanical_pass.f90 +++ b/src/homogenization_mechanical_pass.f90 @@ -17,15 +17,15 @@ module subroutine pass_init ho, & Nmembers - print'(/,a)', ' <<<+- homogenization:mechanical:pass init -+>>>' + print'(/,1x,a)', '<<<+- homogenization:mechanical:pass init -+>>>' - print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_NONE_ID) + print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_NONE_ID) flush(IO_STDOUT) do ho = 1, size(homogenization_type) - if(homogenization_type(ho) /= HOMOGENIZATION_NONE_ID) cycle + if (homogenization_type(ho) /= HOMOGENIZATION_NONE_ID) cycle - if(homogenization_Nconstituents(ho) /= 1) & + if (homogenization_Nconstituents(ho) /= 1) & call IO_error(211,ext_msg='N_constituents (pass)') Nmembers = count(material_homogenizationID == ho) @@ -33,7 +33,7 @@ module subroutine pass_init allocate(homogState(ho)%state0(0,Nmembers)) allocate(homogState(ho)%state (0,Nmembers)) - enddo + end do end subroutine pass_init diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 2d0adc855..c5490d869 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -44,7 +44,7 @@ module subroutine thermal_init() integer :: ho - print'(/,a)', ' <<<+- homogenization:thermal init -+>>>' + print'(/,1x,a)', '<<<+- homogenization:thermal init -+>>>' configHomogenizations => config_material%get('homogenization') @@ -65,9 +65,9 @@ module subroutine thermal_init() #endif else prm%output = emptyStringArray - endif + end if end associate - enddo + end do call pass_init() @@ -89,7 +89,7 @@ module subroutine thermal_partition(ce) dot_T = current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) call phase_thermal_setField(T,dot_T,co,ce) - enddo + end do end subroutine thermal_partition @@ -108,7 +108,7 @@ module function homogenization_mu_T(ce) result(mu) mu = phase_mu_T(1,ce) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) mu = mu + phase_mu_T(co,ce) - enddo + end do mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) @@ -129,7 +129,7 @@ module function homogenization_K_T(ce) result(K) K = phase_K_T(1,ce) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) K = K + phase_K_T(co,ce) - enddo + end do K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) @@ -150,7 +150,7 @@ module function homogenization_f_T(ce) result(f) f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce)) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce)) - enddo + end do f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) @@ -189,7 +189,7 @@ module subroutine thermal_results(ho,group) case('T') call results_writeDataset(current(ho)%T,group,'T','temperature','K') end select - enddo outputsLoop + end do outputsLoop end associate end subroutine thermal_results diff --git a/src/homogenization_thermal_pass.f90 b/src/homogenization_thermal_pass.f90 index 87020d8d2..2673c2789 100644 --- a/src/homogenization_thermal_pass.f90 +++ b/src/homogenization_thermal_pass.f90 @@ -8,7 +8,7 @@ contains module subroutine pass_init() - print'(/,a)', ' <<<+- homogenization:thermal:pass init -+>>>' + print'(/,1x,a)', '<<<+- homogenization:thermal:pass init -+>>>' end subroutine pass_init diff --git a/src/lattice.f90 b/src/lattice.f90 index 65950774e..11d1dca0f 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -409,7 +409,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine lattice_init - print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- lattice init -+>>>'; flush(IO_STDOUT) call selfTest diff --git a/src/material.f90 b/src/material.f90 index 652cc30fb..c83cd468c 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -61,11 +61,11 @@ subroutine material_init(restart) logical, intent(in) :: restart - print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- material init -+>>>'; flush(IO_STDOUT) call parse - print*, 'parsed material.yaml' + print'(/,1x,a)', 'parsed material.yaml' if (.not. restart) then @@ -73,7 +73,7 @@ subroutine material_init(restart) call results_mapping_phase(material_phaseID,material_phaseEntry,material_name_phase) call results_mapping_homogenization(material_homogenizationID,material_homogenizationEntry,material_name_homogenization) call results_closeJobFile - endif + end if end subroutine material_init @@ -113,7 +113,7 @@ subroutine parse() do h=1, homogenizations%length homogenization => homogenizations%get(h) homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents') - enddo + end do homogenization_maxNconstituents = maxval(homogenization_Nconstituents) allocate(counterPhase(phases%length),source=0) @@ -137,7 +137,7 @@ subroutine parse() material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization')) counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1 material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce)) - enddo + end do frac = 0.0_pReal do co = 1, constituents%length @@ -151,12 +151,12 @@ subroutine parse() material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el)) material_phaseID(co,ce) = material_phaseAt(co,el) - enddo + end do - enddo + end do if (dNeq(frac,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent') - enddo + end do allocate(material_O_0(materials%length)) @@ -167,8 +167,8 @@ subroutine parse() do co = 1, constituents%length constituent => constituents%get(co) call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) - enddo - enddo + end do + end do end subroutine parse @@ -186,15 +186,15 @@ subroutine sanityCheck(materials,homogenizations) constituents integer :: m - if(maxval(discretization_materialAt) > materials%length) & + if (maxval(discretization_materialAt) > materials%length) & call IO_error(155,ext_msg='More materials requested than found in material.yaml') do m = 1, materials%length material => materials%get(m) constituents => material%get('constituents') homogenization => homogenizations%get(material%get_asString('homogenization')) - if(constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148) - enddo + if (constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148) + end do end subroutine sanityCheck @@ -215,12 +215,12 @@ function getKeys(dict) do i=1, dict%length temp(i) = dict%getKey(i) l = max(len_trim(temp(i)),l) - enddo + end do allocate(character(l)::getKeys(dict%length)) do i=1, dict%length getKeys(i) = trim(temp(i)) - enddo + end do end function getKeys diff --git a/src/math.f90 b/src/math.f90 index 8e31d3fcd..31bd3b952 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -91,7 +91,7 @@ subroutine math_init class(tNode), pointer :: & num_generic - print'(/,a)', ' <<<+- math init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- math init -+>>>'; flush(IO_STDOUT) num_generic => config_numerics%get('generic',defaultVal=emptyDict) randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0) @@ -109,9 +109,9 @@ subroutine math_init call random_seed(put = randInit) call random_number(randTest) - print'(a,i2)', ' size of random seed: ', randSize - print'(a,i0)', ' value of random seed: ', randInit(1) - print'(a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest + print'(/,a,i2)', ' size of random seed: ', randSize + print'( a,i0)', ' value of random seed: ', randInit(1) + print'( a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest call random_seed(put = randInit) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index a8446d9fa..372cc26ca 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -86,7 +86,7 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) call CPFEM_initAll - print'(/,a)', ' <<<+- DAMASK_mesh init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- DAMASK_mesh init -+>>>'; flush(IO_STDOUT) !--------------------------------------------------------------------- ! reading field information from numerics file and do sanity checks @@ -115,16 +115,16 @@ program DAMASK_mesh case('$Loadcase') N_def = N_def + 1 end select - enddo ! count all identifiers to allocate memory and do sanity check - enddo + end do ! count all identifiers to allocate memory and do sanity check + end do - if(N_def < 1) call IO_error(error_ID = 837) + if (N_def < 1) call IO_error(error_ID = 837) allocate(loadCases(N_def)) do i = 1, size(loadCases) allocate(loadCases(i)%fieldBC(1)) loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID - enddo + end do do i = 1, size(loadCases) loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements @@ -138,12 +138,12 @@ program DAMASK_mesh case (3) loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID end select - enddo + end do do component = 1, loadCases(i)%fieldBC(1)%nComponents allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal) allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) - enddo - enddo + end do + end do !-------------------------------------------------------------------------------------------------- ! reading the load case and assign values to the allocated data structure @@ -163,7 +163,7 @@ program DAMASK_mesh currentFaceSet = -1 do faceSet = 1, mesh_Nboundaries if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet - enddo + end do if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') case('t') loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1) @@ -192,11 +192,11 @@ program DAMASK_mesh .true. loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = & IO_floatValue(line,chunkPos,i+1) - endif - enddo + end if + end do end select - enddo - enddo + end do + end do !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case @@ -204,28 +204,28 @@ program DAMASK_mesh errorID = 0 checkLoadcases: do currentLoadCase = 1, size(loadCases) write (loadcase_string, '(i0)' ) currentLoadCase - print'(a,i0)', ' load case: ', currentLoadCase + print'(/,1x,a,i0)', 'load case: ', currentLoadCase if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & - print'(a)', ' drop guessing along trajectory' - print'(a)', ' Field '//trim(FIELD_MECH_label) + print'(2x,a)', 'drop guessing along trajectory' + print'(2x,a)', 'Field '//trim(FIELD_MECH_label) do faceSet = 1, mesh_Nboundaries do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) & - print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), & - ' Component ', component, & - ' Value ', loadCases(currentLoadCase)%fieldBC(1)% & - componentBC(component)%Value(faceSet) - enddo - enddo - print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time + print'(a,i2,a,i2,a,f12.7)', & + ' Face ', mesh_boundaries(faceSet), & + ' Component ', component, & + ' Value ', loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(faceSet) + end do + end do + print'(2x,a,f12.6)', 'time: ', loadCases(currentLoadCase)%time if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count - print'(a,i5)', ' increments: ', loadCases(currentLoadCase)%incs + print'(2x,a,i5)', 'increments: ', loadCases(currentLoadCase)%incs if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency - print'(a,i5)', ' output frequency: ', & + print'(2x,a,i5)', 'output frequency: ', & loadCases(currentLoadCase)%outputfrequency if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message - enddo checkLoadcases + end do checkLoadcases !-------------------------------------------------------------------------------------------------- ! doing initialization depending on active solvers @@ -235,9 +235,9 @@ program DAMASK_mesh if (worldrank == 0) then open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file - endif + end if - print'(/,a)', ' ... writing initial configuration to file ........................' + print'(/,1x,a)', '... writing initial configuration to file ........................' flush(IO_STDOUT) call CPFEM_results(0,0.0_pReal) @@ -262,7 +262,7 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! report begin of new step - print'(/,a)', ' ###########################################################################' + print'(/,1x,a)', '###########################################################################' print'(1x,a,es12.5,6(a,i0))',& 'Time', time, & 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& @@ -281,60 +281,60 @@ program DAMASK_mesh stagIterate = .true. do while (stagIterate) solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) - if(.not. solres(1)%converged) exit + if (.not. solres(1)%converged) exit stagIter = stagIter + 1 stagIterate = stagIter < stagItMax & .and. all(solres(:)%converged) & .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration - enddo + end do ! check solution cutBack = .False. - if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found + if (.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if (cutBackLevel < maxCutBack) then ! do cut back cutBack = .True. stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator cutBackLevel = cutBackLevel + 1 time = time - timeinc ! rewind time timeinc = timeinc/2.0_pReal - print'(/,a)', ' cutting back' + print'(/,1x,a)', 'cutting back' else ! default behavior, exit if spectral solver does not converge if (worldrank == 0) close(statUnit) call IO_error(950) - endif + end if else guess = .true. ! start guessing after first converged (sub)inc timeIncOld = timeinc - endif + end if if (.not. cutBack .and. worldrank == 0) & write(statUnit,*) totalIncsCounter, time, cutBackLevel, & solres%converged, solres%iterationsNeeded ! write statistics about accepted solution - enddo subStepLooping + end do subStepLooping cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc if (all(solres(:)%converged)) then - print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged' + print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' converged' else - print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged' - endif; flush(IO_STDOUT) + print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' NOT converged' + end if; flush(IO_STDOUT) if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency - print'(/,a)', ' ... writing results to file ......................................' + print'(/,1x,a)', '... writing results to file ......................................' call FEM_mechanical_updateCoords call CPFEM_results(totalIncsCounter,time) - endif + end if - enddo incLooping + end do incLooping - enddo loadCaseLooping + end do loadCaseLooping !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation - print'(/,a)', ' ###########################################################################' + print'(/,1x,a)', '###########################################################################' if (worldrank == 0) close(statUnit) call quit(0) ! no complains ;) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 0c470ca6b..f519dc377 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -41,10 +41,10 @@ contains !-------------------------------------------------------------------------------------------------- subroutine FEM_quadrature_init() - print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6) + print'(/,1x,a)', '<<<+- FEM_quadrature init -+>>>'; flush(6) - print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009' - print*, 'https://www.jstor.org/stable/43693493' + print'(/,1x,a)', 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009' + print'( 1x,a)', 'https://www.jstor.org/stable/43693493' !-------------------------------------------------------------------------------------------------- ! 2D linear diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index a3856ccaa..d5196a9f4 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -96,7 +96,7 @@ subroutine FEM_utilities_init logical :: debugPETSc !< use some in debug defined options for more verbose PETSc solution - print'(/,a)', ' <<<+- FEM_utilities init -+>>>' + print'(/,1x,a)', '<<<+- FEM_utilities init -+>>>' num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) @@ -111,10 +111,10 @@ subroutine FEM_utilities_init debug_mesh => config_debug%get('mesh',defaultVal=emptyList) debugPETSc = debug_mesh%contains('PETSc') - if(debugPETSc) print'(3(/,a),/)', & - ' Initializing PETSc with debug options: ', & + if(debugPETSc) print'(3(/,1x,a),/)', & + 'Initializing PETSc with debug options: ', & trim(PETScDebug), & - ' add more using the "PETSc_options" keyword in numerics.yaml' + 'add more using the "PETSc_options" keyword in numerics.yaml' flush(IO_STDOUT) call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr) CHKERRQ(ierr) @@ -151,7 +151,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) PetscErrorCode :: ierr - print'(/,a)', ' ... evaluating constitutive response ......................................' + print'(/,1x,a)', '... evaluating constitutive response ......................................' call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field if (.not. terminallyIll) & diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index f1d38760d..2446bc0de 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -51,7 +51,7 @@ module discretization_mesh real(pReal), pointer, dimension(:) :: & mesh_node0_temp - + real(pReal), dimension(:,:,:), allocatable :: & mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) @@ -88,7 +88,7 @@ subroutine discretization_mesh_init(restart) integer :: p_i !< integration order (quadrature rule) type(tvec) :: coords_node0 - print'(/,a)', ' <<<+- discretization_mesh init -+>>>' + print'(/,1x,a)', '<<<+- discretization_mesh init -+>>>' !-------------------------------------------------------------------------------- ! read numerics parameter @@ -106,6 +106,7 @@ subroutine discretization_mesh_init(restart) CHKERRQ(ierr) call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr) CHKERRQ(ierr) + print'()' call DMView(globalMesh, PETSC_VIEWER_STDOUT_WORLD,ierr) CHKERRQ(ierr) @@ -173,7 +174,7 @@ subroutine discretization_mesh_init(restart) reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) - call writeGeometry(reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]),mesh_node0) + call writeGeometry(reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]),mesh_node0) end subroutine discretization_mesh_init @@ -249,18 +250,18 @@ subroutine writeGeometry(coordinates_points,coordinates_nodes) real(pReal), dimension(:,:), intent(in) :: & coordinates_nodes, & coordinates_points - + call results_openJobFile call results_closeGroup(results_addGroup('geometry')) - + call results_writeDataset(coordinates_nodes,'geometry','x_n', & 'initial coordinates of the nodes','m') - + call results_writeDataset(coordinates_points,'geometry','x_p', & 'initial coordinates of the materialpoints (cell centers)','m') - + call results_closeJobFile - + end subroutine writeGeometry end module discretization_mesh diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 496a82fc5..c1c0d8e3d 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -48,9 +48,9 @@ module mesh_mechanical_FEM real(pReal) :: & eps_struct_atol, & !< absolute tolerance for mechanical equilibrium eps_struct_rtol !< relative tolerance for mechanical equilibrium - end type tNumerics + end type tNumerics - type(tNumerics), private :: num + type(tNumerics), private :: num !-------------------------------------------------------------------------------------------------- ! PETSc data SNES :: mechanical_snes @@ -112,8 +112,8 @@ subroutine FEM_mechanical_init(fieldBC) real(pReal), dimension(3,3) :: devNull class(tNode), pointer :: & num_mesh - - print'(/,a)', ' <<<+- FEM_mech init -+>>>'; flush(IO_STDOUT) + + print'(/,1x,a)', '<<<+- FEM_mech init -+>>>'; flush(IO_STDOUT) !----------------------------------------------------------------------------- ! read numerical parametes and do sanity checks @@ -123,7 +123,7 @@ subroutine FEM_mechanical_init(fieldBC) num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.) num%eps_struct_atol = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal) num%eps_struct_rtol = num_mesh%get_asFloat('eps_struct_rtol', defaultVal = 1.0e-4_pReal) - + if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%eps_struct_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_rtol') if (num%eps_struct_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_atol') @@ -302,7 +302,7 @@ type(tSolutionState) function FEM_mechanical_solution( & CHKERRQ(ierr) endif - print'(/,a)', ' ===========================================================================' + print'(/,1x,a)', '===========================================================================' flush(IO_STDOUT) end function FEM_mechanical_solution @@ -362,7 +362,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate field derivatives do cell = cellStart, cellEnd-1 !< loop over all elements - + call PetscSectionGetNumFields(section,numFields,ierr) CHKERRQ(ierr) call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element @@ -663,8 +663,8 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), & ' @ Iteration ',PETScIter,' mechanical residual norm = ', & int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) - print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & - ' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal + print'(/,1x,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_mechanical_converged @@ -679,7 +679,7 @@ subroutine FEM_mechanical_updateCoords() nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes) real(pReal), pointer, dimension(:,:) :: & nodeCoords !< nodal coordinates (3,Nnodes) - real(pReal), pointer, dimension(:,:,:) :: & + real(pReal), pointer, dimension(:,:,:) :: & ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems) integer :: & @@ -720,7 +720,7 @@ subroutine FEM_mechanical_updateCoords() call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call PetscDSGetTabulation(mechQuad,0,basisField,basisFieldDer,ierr); CHKERRQ(ierr) allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal) - do c=cellStart,cellEnd-1 + do c=cellStart,cellEnd-1 qOffset=0 call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) !< get nodal coordinates of each element do qPt=0,nQuadrature-1 @@ -737,8 +737,8 @@ subroutine FEM_mechanical_updateCoords() enddo enddo enddo - call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) - end do + call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) + end do call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature])) call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) diff --git a/src/parallelization.f90 b/src/parallelization.f90 index 8ea22ef0a..ebb9573d4 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -76,11 +76,11 @@ subroutine parallelization_init call MPI_Comm_rank(MPI_COMM_WORLD,worldrank,err) if (err /= 0) error stop 'Could not determine worldrank' - if (worldrank == 0) print'(/,a)', ' <<<+- parallelization init -+>>>' + if (worldrank == 0) print'(/,1x,a)', '<<<+- parallelization init -+>>>' call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err) if (err /= 0) error stop 'Could not determine worldsize' - if (worldrank == 0) print'(a,i3)', ' MPI processes: ',worldsize + if (worldrank == 0) print'(/,1x,a,i3)', 'MPI processes: ',worldsize call MPI_Type_size(MPI_INTEGER,typeSize,err) if (err /= 0) error stop 'Could not determine MPI integer size' @@ -97,16 +97,16 @@ subroutine parallelization_init !$ call get_environment_variable(name='OMP_NUM_THREADS',value=NumThreadsString,STATUS=got_env) !$ if(got_env /= 0) then -!$ print*, 'Could not get $OMP_NUM_THREADS, using default' +!$ print'(1x,a)', 'Could not get $OMP_NUM_THREADS, using default' !$ OMP_NUM_THREADS = 4_pI32 !$ else !$ read(NumThreadsString,'(i6)') OMP_NUM_THREADS !$ if (OMP_NUM_THREADS < 1_pI32) then -!$ print*, 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default' +!$ print'(1x,a)', 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default' !$ OMP_NUM_THREADS = 4_pI32 !$ endif !$ endif -!$ print'(a,i2)', ' OMP_NUM_THREADS: ',OMP_NUM_THREADS +!$ print'(1x,a,1x,i2)', 'OMP_NUM_THREADS:',OMP_NUM_THREADS !$ call omp_set_num_threads(OMP_NUM_THREADS) end subroutine parallelization_init diff --git a/src/phase.f90 b/src/phase.f90 index 7b19254b7..a36fc26b9 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -343,7 +343,7 @@ subroutine phase_init phase - print'(/,a)', ' <<<+- phase init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase init -+>>>'; flush(IO_STDOUT) debug_constitutive => config_debug%get('phase', defaultVal=emptyList) debugConstitutive%basic = debug_constitutive%contains('basic') @@ -371,20 +371,20 @@ subroutine phase_init phase_cOverA(ph) = phase%get_asFloat('c/a') phase_rho(ph) = phase%get_asFloat('rho',defaultVal=0.0_pReal) allocate(phase_O_0(ph)%data(count(material_phaseID==ph))) - enddo + end do do ce = 1, size(material_phaseID,2) ma = discretization_materialAt((ce-1)/discretization_nIPs+1) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) ph = material_phaseID(co,ce) phase_O_0(ph)%data(material_phaseEntry(co,ce)) = material_O_0(ma)%data(co) - enddo - enddo + end do + end do allocate(phase_O(phases%length)) do ph = 1,phases%length phase_O(ph)%data = phase_O_0(ph)%data - enddo + end do call mechanical_init(phases) call damage_init @@ -471,7 +471,7 @@ subroutine phase_results() call mechanical_results(group,ph) call damage_results(group,ph) - enddo + end do end subroutine phase_results @@ -495,7 +495,7 @@ subroutine crystallite_init() phases - print'(/,a)', ' <<<+- crystallite init -+>>>' + print'(/,1x,a)', '<<<+- crystallite init -+>>>' cMax = homogenization_maxNconstituents iMax = discretization_nIPs @@ -515,28 +515,28 @@ subroutine crystallite_init() num%nState = num_crystallite%get_asInt ('nState', defaultVal=20) num%nStress = num_crystallite%get_asInt ('nStress', defaultVal=40) - if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst') - if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst') - if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst') + if (num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst') + if (num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst') + if (num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst') - if(num%subStepSizeLp <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLp') - if(num%subStepSizeLi <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLi') + if (num%subStepSizeLp <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLp') + if (num%subStepSizeLi <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLi') - if(num%rtol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState') - if(num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress') - if(num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress') + if (num%rtol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState') + if (num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress') + if (num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress') - if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') + if (num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') - if(num%nState < 1) call IO_error(301,ext_msg='nState') - if(num%nStress< 1) call IO_error(301,ext_msg='nStress') + if (num%nState < 1) call IO_error(301,ext_msg='nState') + if (num%nStress< 1) call IO_error(301,ext_msg='nStress') phases => config_material%get('phase') - print'(a42,1x,i10)', ' # of elements: ', eMax - print'(a42,1x,i10)', ' # of integration points/element: ', iMax - print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax + 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) @@ -547,9 +547,9 @@ subroutine crystallite_init() do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) call crystallite_orientations(co,ip,el) call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states - enddo - enddo - enddo + end do + end do + end do !$OMP END PARALLEL DO @@ -642,7 +642,7 @@ subroutine phase_restartWrite(fileHandle) call HDF5_closeGroup(groupHandle(2)) - enddo + end do call HDF5_closeGroup(groupHandle(1)) @@ -670,7 +670,7 @@ subroutine phase_restartRead(fileHandle) call HDF5_closeGroup(groupHandle(2)) - enddo + end do call HDF5_closeGroup(groupHandle(1)) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index b3f2c7c22..49b9f9528 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -83,7 +83,7 @@ module subroutine damage_init source logical:: damage_active - print'(/,a)', ' <<<+- phase:damage init -+>>>' + print'(/,1x,a)', '<<<+- phase:damage init -+>>>' phases => config_material%get('phase') @@ -108,16 +108,16 @@ module subroutine damage_init param(ph)%D(1,1) = source%get_asFloat('D_11') if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33') param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph)) - endif + end if - enddo + end do allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID) if (damage_active) then where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID - endif + end if phase_damage_maxSizeDotState = maxval(damageState%sizeDotState) @@ -181,7 +181,7 @@ module subroutine damage_restore(ce) if (damageState(material_phaseID(co,ce))%sizeState > 0) & damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = & damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce)) - enddo + end do end subroutine damage_restore @@ -242,11 +242,11 @@ function integrateDamageState(Delta_t,ph,en) result(broken) if (damageState(ph)%sizeState == 0) then broken = .false. return - endif + end if converged_ = .true. broken = phase_damage_collectDotState(ph,en) - if(broken) return + if (broken) return size_so = damageState(ph)%sizeDotState damageState(ph)%state(1:size_so,en) = damageState(ph)%state0 (1:size_so,en) & @@ -255,11 +255,11 @@ function integrateDamageState(Delta_t,ph,en) result(broken) iteration: do NiterationState = 1, num%nState - if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) + if (nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en) broken = phase_damage_collectDotState(ph,en) - if(broken) exit iteration + if (broken) exit iteration zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) @@ -274,12 +274,12 @@ function integrateDamageState(Delta_t,ph,en) result(broken) damageState(ph)%atol(1:size_so)) - if(converged_) then + if (converged_) then broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en) exit iteration - endif + end if - enddo iteration + end do iteration broken = broken .or. .not. converged_ @@ -302,7 +302,7 @@ function integrateDamageState(Delta_t,ph,en) result(broken) damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) else damper = 1.0_pReal - endif + end if end function damper @@ -358,7 +358,7 @@ function phase_damage_collectDotState(ph,en) result(broken) broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en))) - endif + end if end function phase_damage_collectDotState @@ -385,7 +385,7 @@ module function phase_K_phi(co,ce) result(K) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: K real(pReal), parameter :: l = 1.0_pReal - + K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%D) & * l**2.0_pReal @@ -419,12 +419,12 @@ function phase_damage_deltaState(Fe, ph, en) result(broken) case (DAMAGE_ISOBRITTLE_ID) sourceType call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) - if(.not. broken) then + if (.not. broken) then myOffset = damageState(ph)%offsetDeltaState mySize = damageState(ph)%sizeDeltaState damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = & damageState(ph)%state(myOffset + 1: myOffset + mySize,en) + damageState(ph)%deltaState(1:mySize,en) - endif + end if end select sourceType @@ -454,7 +454,7 @@ function source_active(source_label) result(active_source) sources => phase%get('damage',defaultVal=emptyList) src => sources%get(1) active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label - enddo + end do end function source_active @@ -497,7 +497,7 @@ module subroutine damage_forward() do ph = 1, size(damageState) if (damageState(ph)%sizeState > 0) & damageState(ph)%state0 = damageState(ph)%state - enddo + end do end subroutine damage_forward diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index e8f4e7f98..f6facbe48 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -46,10 +46,10 @@ module function anisobrittle_init() result(mySources) mySources = source_active('anisobrittle') - if(count(mySources) == 0) return + if (count(mySources) == 0) return - print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>' - print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:damage:anisobrittle init -+>>>' + print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) phases => config_material%get('phase') @@ -57,7 +57,7 @@ module function anisobrittle_init() result(mySources) do ph = 1, phases%length - if(mySources(ph)) then + if (mySources(ph)) then phase => phases%get(ph) sources => phase%get('damage') @@ -94,14 +94,14 @@ module function anisobrittle_init() result(mySources) Nmembers = count(material_phaseID==ph) call phase_allocateState(damageState(ph),Nmembers,1,1,0) damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) - if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' + if (any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' end associate if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)') - endif + end if - enddo + end do end function anisobrittle_init @@ -136,7 +136,7 @@ module subroutine anisobrittle_dotState(S, ph,en) * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q) - enddo + end do end associate end subroutine anisobrittle_dotState @@ -159,7 +159,7 @@ module subroutine anisobrittle_results(phase,group) case ('f_phi') call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') end select - enddo outputsLoop + end do outputsLoop end associate end subroutine anisobrittle_results @@ -200,7 +200,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i) - endif + end if traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) if (abs(traction_t) > traction_crit + tol_math_check) then @@ -210,7 +210,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i) - endif + end if traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) if (abs(traction_n) > traction_crit + tol_math_check) then @@ -220,8 +220,8 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i) - endif - enddo + end if + end do end associate end subroutine damage_anisobrittle_LiAndItsTangent diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 066cab47e..25ce1d3e9 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -44,10 +44,10 @@ module function isobrittle_init() result(mySources) mySources = source_active('isobrittle') - if(count(mySources) == 0) return + if (count(mySources) == 0) return - print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>' - print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:damage:isobrittle init -+>>>' + print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) phases => config_material%get('phase') @@ -56,7 +56,7 @@ module function isobrittle_init() result(mySources) allocate(deltaState(phases%length)) do ph = 1, phases%length - if(mySources(ph)) then + if (mySources(ph)) then phase => phases%get(ph) sources => phase%get('damage') @@ -77,7 +77,7 @@ module function isobrittle_init() result(mySources) Nmembers = count(material_phaseID==ph) call phase_allocateState(damageState(ph),Nmembers,1,1,1) damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) - if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' + if (any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' stt%r_W => damageState(ph)%state(1,:) dlt%r_W => damageState(ph)%deltaState(1,:) @@ -86,9 +86,9 @@ module function isobrittle_init() result(mySources) if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)') - endif + end if - enddo + end do end function isobrittle_init @@ -141,7 +141,7 @@ module subroutine isobrittle_results(phase,group) case ('f_phi') call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') ! Wrong, this is dimensionless end select - enddo outputsLoop + end do outputsLoop end associate diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 6ad96665c..4186b11cf 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -214,7 +214,7 @@ module subroutine mechanical_init(phases) phase, & mech - print'(/,a)', ' <<<+- phase:mechanical init -+>>>' + print'(/,1x,a)', '<<<+- phase:mechanical init -+>>>' !------------------------------------------------------------------------------------------------- allocate(output_constituent(phases%length)) diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index d33adc5d1..c73665d17 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -44,7 +44,7 @@ module subroutine eigen_init(phases) kinematics, & mechanics - print'(/,a)', ' <<<+- phase:mechanical:eigen init -+>>>' + print'(/,1x,a)', '<<<+- phase:mechanical:eigen init -+>>>' !-------------------------------------------------------------------------------------------------- ! explicit eigen mechanisms @@ -55,11 +55,11 @@ module subroutine eigen_init(phases) mechanics => phase%get('mechanical') kinematics => mechanics%get('eigen',defaultVal=emptyList) Nmodels(ph) = kinematics%length - enddo + end do allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID) - if(maxval(Nmodels) /= 0) then + if (maxval(Nmodels) /= 0) then where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID endif @@ -97,8 +97,8 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki do k = 1, kinematics%length kinematics_type => kinematics%get(k) active_kinematics(k,ph) = kinematics_type%get_asString('type') == kinematics_label - enddo - enddo + end do + end do end function kinematics_active @@ -125,11 +125,11 @@ function kinematics_active2(kinematics_label) result(active_kinematics) do ph = 1, phases%length phase => phases%get(ph) kinematics => phase%get('damage',defaultVal=emptyList) - if(kinematics%length < 1) return + if (kinematics%length < 1) return kinematics_type => kinematics%get(1) if (.not. kinematics_type%contains('type')) continue active_kinematics(ph) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label - enddo + end do end function kinematics_active2 @@ -188,7 +188,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & dLi_dS = dLi_dS + my_dLi_dS active = .true. end select kinematicsType - enddo KinematicsLoop + end do KinematicsLoop select case (model_damage(ph)) case (KINEMATICS_cleavage_opening_ID) @@ -198,7 +198,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & active = .true. end select - if(.not. active) return + if (.not. active) return FiInv = math_inv33(Fi) detFi = math_det33(Fi) @@ -209,7 +209,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) - enddo; enddo + end do; end do end subroutine phase_LiAndItsTangents diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index 4243d09a4..1bf231c2c 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -21,8 +21,8 @@ module function damage_anisobrittle_init() result(myKinematics) myKinematics = kinematics_active2('anisobrittle') if(count(myKinematics) == 0) return - print'(/,a)', ' <<<+- phase:mechanical:eigen:cleavageopening init -+>>>' - print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:mechanical:eigen:cleavageopening init -+>>>' + print'(/,a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) end function damage_anisobrittle_init diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index 0976c1dfa..1c08140ae 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -36,25 +36,25 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) kinematics, & kinematic_type - print'(/,a)', ' <<<+- phase:mechanical:eigen:thermalexpansion init -+>>>' + print'(/,1x,a)', '<<<+- phase:mechanical:eigen:thermalexpansion init -+>>>' myKinematics = kinematics_active('thermalexpansion',kinematics_length) Ninstances = count(myKinematics) - print'(a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + print'(/,a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT) + if (Ninstances == 0) return phases => config_material%get('phase') allocate(param(Ninstances)) allocate(kinematics_thermal_expansion_instance(phases%length), source=0) do p = 1, phases%length - if(any(myKinematics(:,p))) kinematics_thermal_expansion_instance(p) = count(myKinematics(:,1:p)) + if (any(myKinematics(:,p))) kinematics_thermal_expansion_instance(p) = count(myKinematics(:,1:p)) phase => phases%get(p) - if(count(myKinematics(:,p)) == 0) cycle + if (count(myKinematics(:,p)) == 0) cycle mech => phase%get('mechanical') kinematics => mech%get('eigen') do k = 1, kinematics%length - if(myKinematics(k,p)) then + if (myKinematics(k,p)) then associate(prm => param(kinematics_thermal_expansion_instance(p))) kinematic_type => kinematics%get(k) @@ -67,14 +67,14 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) prm%A(3,3,1) = kinematic_type%get_asFloat('A_33') prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T',defaultVal=0.0_pReal) prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) - endif + end if do i=1, size(prm%A,3) prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p)) - enddo + end do end associate - endif - enddo - enddo + end if + end do + end do end function thermalexpansion_init @@ -92,7 +92,7 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) real(pReal) :: T, dot_T - + T = thermal_T(ph,me) dot_T = thermal_dot_T(ph,me) diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index f253a0f5d..8b05adc31 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -26,10 +26,10 @@ module subroutine elastic_init(phases) elastic - print'(/,a)', ' <<<+- phase:mechanical:elastic init -+>>>' - print'(/,a)', ' <<<+- phase:mechanical:elastic:Hooke init -+>>>' + print'(/,1x,a)', '<<<+- phase:mechanical:elastic init -+>>>' + print'(/,1x,a)', '<<<+- phase:mechanical:elastic:Hooke init -+>>>' - print'(a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT) + print'(/,a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT) allocate(param(phases%length)) @@ -48,18 +48,18 @@ module subroutine elastic_init(phases) if (any(phase_lattice(ph) == ['hP','tI'])) then prm%C66(1,3) = elastic%get_asFloat('C_13') prm%C66(3,3) = elastic%get_asFloat('C_33') - endif + end if if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66') prm%C66 = lattice_symmetrize_C66(prm%C66,phase_lattice(ph)) prm%nu = lattice_equivalent_nu(prm%C66,'voigt') prm%mu = lattice_equivalent_mu(prm%C66,'voigt') - + prm%C66 = math_sym3333to66(math_Voigt66to3333(prm%C66)) ! Literature data is in Voigt notation end associate - enddo + end do end subroutine elastic_init @@ -98,7 +98,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & do i =1, 3;do j=1,3 dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn - enddo; enddo + end do; end do end subroutine phase_hooke_SandItsTangents diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 75da4fcdd..d0cc3a261 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -222,7 +222,7 @@ contains module subroutine plastic_init - print'(/,a)', ' <<<+- phase:mechanical:plastic init -+>>>' + print'(/,1x,a)', '<<<+- phase:mechanical:plastic init -+>>>' where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID @@ -295,7 +295,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) - enddo; enddo + end do; end do end subroutine plastic_LpAndItsTangents @@ -418,7 +418,7 @@ module function plastic_deltaState(ph, en) result(broken) end select plasticType - if(.not. broken) then + if (.not. broken) then select case(phase_plasticity(ph)) case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) @@ -427,7 +427,7 @@ module function plastic_deltaState(ph, en) result(broken) plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = & plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en) end select - endif + end if end function plastic_deltaState @@ -453,8 +453,8 @@ function plastic_active(plastic_label) result(active_plastic) phase => phases%get(ph) mech => phase%get('mechanical') pl => mech%get('plastic',defaultVal = emptyDict) - if(pl%get_asString('type',defaultVal='none') == plastic_label) active_plastic(ph) = .true. - enddo + if (pl%get_asString('type',defaultVal='none') == plastic_label) active_plastic(ph) = .true. + end do end function plastic_active diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index c759cdaad..33bec98cc 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -99,13 +99,13 @@ module function plastic_dislotungsten_init() result(myPlasticity) myPlasticity = plastic_active('dislotungsten') - if(count(myPlasticity) == 0) return + if (count(myPlasticity) == 0) return - print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotungsten init -+>>>' - print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotungsten init -+>>>' + print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) - print*, 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016' - print*, 'https://doi.org/10.1016/j.ijplas.2015.09.002' + print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016' + print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002' phases => config_material%get('phase') @@ -116,7 +116,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) do ph = 1, phases%length - if(.not. myPlasticity(ph)) cycle + if (.not. myPlasticity(ph)) cycle associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) @@ -243,7 +243,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:) dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) allocate(dst%tau_pass(prm%sum_N_sl,Nmembers), source=0.0_pReal) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index de73cee04..d7ddff80a 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -150,17 +150,17 @@ module function plastic_dislotwin_init() result(myPlasticity) myPlasticity = plastic_active('dislotwin') if(count(myPlasticity) == 0) return - print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotwin init -+>>>' - print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotwin init -+>>>' + print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) - print*, 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004' - print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL + print'(/,1x,a)', 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004' + print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL - print*, 'F. Roters et al., Computational Materials Science 39:91–95, 2007' - print*, 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL + print'(/,1x,a)', 'F. Roters et al., Computational Materials Science 39:91–95, 2007' + print'( 1x,a)', 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL - print*, 'S.L. Wong et al., Acta Materialia 118:140–151, 2016' - print*, 'https://doi.org/10.1016/j.actamat.2016.07.032' + print'(/,1x,a)', 'S.L. Wong et al., Acta Materialia 118:140–151, 2016' + print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2016.07.032' phases => config_material%get('phase') @@ -256,7 +256,7 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray) allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0)) - endif slipActive + end if slipActive !-------------------------------------------------------------------------------------------------- ! twin related parameters @@ -283,7 +283,7 @@ module function plastic_dislotwin_init() result(myPlasticity) if (.not. prm%fccTwinTransNucleation) then prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw') prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,N_tw) - endif + end if ! expand: family => system prm%b_tw = math_expand(prm%b_tw,N_tw) @@ -299,11 +299,11 @@ module function plastic_dislotwin_init() result(myPlasticity) if (any(prm%r < 0.0_pReal)) extmsg = trim(extmsg)//' p_tw' if (.not. prm%fccTwinTransNucleation) then if (any(prm%dot_N_0_tw < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tw' - endif + end if else twinActive allocate(prm%gamma_char,prm%b_tw,prm%dot_N_0_tw,prm%t_tw,prm%r,source=emptyRealArray) allocate(prm%h_tw_tw(0,0)) - endif twinActive + end if twinActive !-------------------------------------------------------------------------------------------------- ! transformation related parameters @@ -335,7 +335,7 @@ module function plastic_dislotwin_init() result(myPlasticity) if (phase_lattice(ph) /= 'cF') then prm%dot_N_0_tr = pl%get_as1dFloat('dot_N_0_tr') prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr) - endif + end if prm%t_tr = pl%get_as1dFloat('t_tr') prm%t_tr = math_expand(prm%t_tr,N_tr) prm%s = pl%get_as1dFloat('p_tr',defaultVal=[0.0_pReal]) @@ -349,11 +349,11 @@ module function plastic_dislotwin_init() result(myPlasticity) if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr' if (phase_lattice(ph) /= 'cF') then if (any(prm%dot_N_0_tr < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tr' - endif + end if else transActive allocate(prm%s,prm%b_tr,prm%t_tr,prm%dot_N_0_tr,source=emptyRealArray) allocate(prm%h_tr_tr(0,0)) - endif transActive + end if transActive !-------------------------------------------------------------------------------------------------- ! shearband related parameters @@ -369,7 +369,7 @@ module function plastic_dislotwin_init() result(myPlasticity) if (prm%E_sb < 0.0_pReal) extmsg = trim(extmsg)//' Q_sb' if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb' if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb' - endif + end if !-------------------------------------------------------------------------------------------------- ! parameters required for several mechanisms and their interactions @@ -383,19 +383,19 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%T_ref = pl%get_asFloat('T_ref') prm%Gamma_sf(1) = pl%get_asFloat('Gamma_sf') prm%Gamma_sf(2) = pl%get_asFloat('Gamma_sf,T',defaultVal=0.0_pReal) - endif + end if slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), & phase_lattice(ph)) if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation' - endif slipAndTwinActive + end if slipAndTwinActive slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,pl%get_as1dFloat('h_sl-tr'), & phase_lattice(ph)) if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation' - endif slipAndTransActive + end if slipAndTransActive !-------------------------------------------------------------------------------------------------- ! allocate state arrays @@ -465,7 +465,7 @@ module function plastic_dislotwin_init() result(myPlasticity) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(dislotwin)') - enddo + end do end function plastic_dislotwin_init @@ -494,11 +494,11 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) do i=1,prm%sum_N_tw homogenizedC = homogenizedC & + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i) - enddo + end do do i=1,prm%sum_N_tr homogenizedC = homogenizedC & + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i) - enddo + end do end associate @@ -566,7 +566,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) - enddo slipContribution + end do slipContribution call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) twinContibution: do i = 1, prm%sum_N_tw @@ -574,7 +574,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) - enddo twinContibution + end do twinContibution call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) transContibution: do i = 1, prm%sum_N_tr @@ -582,7 +582,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) - enddo transContibution + end do transContibution Lp = Lp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated @@ -608,10 +608,10 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) - endif significantShearBandStress - enddo + end if significantShearBandStress + end do - endif shearBandingContribution + end if shearBandingContribution end associate @@ -686,9 +686,9 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & / (d_hat-prm%d_caron(i)) - endif - endif significantSlipStress - enddo slipState + end if + end if significantSlipStress + end do slipState dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & - dot_rho_dip_formation & @@ -833,7 +833,7 @@ module subroutine plastic_dislotwin_results(ph,group) end select - enddo + end do end associate @@ -968,8 +968,8 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& end if else isFCC Ndot0=prm%dot_N_0_tw(i) - endif isFCC - enddo + end if isFCC + end do significantStress: where(tau > tol_math_check) StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r @@ -1037,8 +1037,8 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& end if else isFCC Ndot0=prm%dot_N_0_tr(i) - endif isFCC - enddo + end if isFCC + end do significantStress: where(tau > tol_math_check) StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 8b10ed5ef..bea5339c7 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -68,11 +68,11 @@ module function plastic_isotropic_init() result(myPlasticity) myPlasticity = plastic_active('isotropic') if(count(myPlasticity) == 0) return - print'(/,a)', ' <<<+- phase:mechanical:plastic:isotropic init -+>>>' - print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:mechanical:plastic:isotropic init -+>>>' + print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) - print*, 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018' - print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047' + print'(/,a)', 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018' + print'(/,a)', 'https://doi.org/10.1016/j.scriptamat.2017.09.047' phases => config_material%get('phase') allocate(param(phases%length)) @@ -140,7 +140,7 @@ module function plastic_isotropic_init() result(myPlasticity) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(isotropic)') - enddo + end do end function plastic_isotropic_init @@ -232,7 +232,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) else Li = 0.0_pReal dLi_dMi = 0.0_pReal - endif + end if end associate @@ -262,7 +262,7 @@ module subroutine isotropic_dotState(Mp,ph,en) norm_Mp = sqrt(math_tensordot(Mp,Mp)) else norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))) - endif + end if dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n @@ -273,13 +273,13 @@ module subroutine isotropic_dotState(Mp,ph,en) xi_inf_star = prm%xi_inf & + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) & / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n) - endif + end if dot%xi(en) = dot_gamma & * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) & * sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star) else dot%xi(en) = 0.0_pReal - endif + end if end associate @@ -303,7 +303,7 @@ module subroutine plastic_isotropic_results(ph,group) call results_writeDataset(stt%xi,group,trim(prm%output(o)), & 'resistance against plastic flow','Pa') end select - enddo outputsLoop + end do outputsLoop end associate end subroutine plastic_isotropic_results diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 276212576..6c0e1e0cc 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -83,8 +83,8 @@ module function plastic_kinehardening_init() result(myPlasticity) myPlasticity = plastic_active('kinehardening') if(count(myPlasticity) == 0) return - print'(/,a)', ' <<<+- phase:mechanical:plastic:kinehardening init -+>>>' - print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:mechanical:plastic:kinehardening init -+>>>' + print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) phases => config_material%get('phase') @@ -95,7 +95,7 @@ module function plastic_kinehardening_init() result(myPlasticity) do ph = 1, phases%length - if(.not. myPlasticity(ph)) cycle + if (.not. myPlasticity(ph)) cycle associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph)) @@ -125,7 +125,7 @@ module function plastic_kinehardening_init() result(myPlasticity) else prm%P_nS_pos = prm%P prm%P_nS_neg = prm%P - endif + end if prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & phase_lattice(ph)) @@ -161,7 +161,7 @@ module function plastic_kinehardening_init() result(myPlasticity) xi_0 = emptyRealArray allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray) allocate(prm%h_sl_sl(0,0)) - endif slipActive + end if slipActive !-------------------------------------------------------------------------------------------------- ! allocate state arrays @@ -217,7 +217,7 @@ module function plastic_kinehardening_init() result(myPlasticity) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(kinehardening)') - enddo + end do end function plastic_kinehardening_init @@ -258,7 +258,7 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_pos(i) * prm%P(k,l,i) * prm%P_nS_pos(m,n,i) & + ddot_gamma_dtau_neg(i) * prm%P(k,l,i) * prm%P_nS_neg(m,n,i) - enddo + end do end associate @@ -382,7 +382,7 @@ module subroutine plastic_kinehardening_results(ph,group) 'plastic shear','1',prm%systems_sl) end select - enddo + end do end associate @@ -424,7 +424,7 @@ pure subroutine kinetics(Mp,ph,en, & tau_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - stt%chi(i,en) tau_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)) - stt%chi(i,en), & 0.0_pReal, prm%nonSchmidActive) - enddo + end do where(dNeq0(tau_pos)) dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active @@ -446,14 +446,14 @@ pure subroutine kinetics(Mp,ph,en, & else where ddot_gamma_dtau_pos = 0.0_pReal end where - endif + end if if (present(ddot_gamma_dtau_neg)) then where(dNeq0(dot_gamma_neg)) ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg else where ddot_gamma_dtau_neg = 0.0_pReal end where - endif + end if end associate diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 index 544de31b1..50b4b4cbc 100644 --- a/src/phase_mechanical_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -22,16 +22,16 @@ module function plastic_none_init() result(myPlasticity) myPlasticity = plastic_active('none') - if(count(myPlasticity) == 0) return + if (count(myPlasticity) == 0) return - print'(/,a)', ' <<<+- phase:mechanical:plastic:none init -+>>>' - print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:mechanical:plastic:none init -+>>>' + print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) phases => config_material%get('phase') do ph = 1, phases%length - if(.not. myPlasticity(ph)) cycle + if (.not. myPlasticity(ph)) cycle call phase_allocateState(plasticState(ph),count(material_phaseID == ph),0,0,0) - enddo + end do end function plastic_none_init diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 1968d258e..727278e18 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -197,19 +197,19 @@ module function plastic_nonlocal_init() result(myPlasticity) myPlasticity = plastic_active('nonlocal') Ninstances = count(myPlasticity) - if(Ninstances == 0) then + if (Ninstances == 0) then call geometry_plastic_nonlocal_disable return - endif + end if - print'(/,a)', ' <<<+- phase:mechanical:plastic:nonlocal init -+>>>' - print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:mechanical:plastic:nonlocal init -+>>>' + print'(/,a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) - print*, 'C. Reuber et al., Acta Materialia 71:333–348, 2014' - print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL + print'(/,1x,a)', 'C. Reuber et al., Acta Materialia 71:333–348, 2014' + print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL - print*, 'C. Kords, Dissertation RWTH Aachen, 2014' - print*, 'http://publications.rwth-aachen.de/record/229993' + print'(/,1x,a)', 'C. Kords, Dissertation RWTH Aachen, 2014' + print'( 1x,a)', 'http://publications.rwth-aachen.de/record/229993' phases => config_material%get('phase') @@ -224,7 +224,7 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(dependentState(phases%length)) do ph = 1, phases%length - if(.not. myPlasticity(ph)) cycle + if (.not. myPlasticity(ph)) cycle associate(prm => param(ph), dot => dotState(ph), stt => state(ph), & st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph)) @@ -259,7 +259,7 @@ module function plastic_nonlocal_init() result(myPlasticity) else prm%P_nS_pos = prm%P_sl prm%P_nS_neg = prm%P_sl - endif + end if prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dFloat('h_sl-sl'), & phase_lattice(ph)) @@ -280,8 +280,8 @@ module function plastic_nonlocal_init() result(myPlasticity) if (all(dEq0 (math_cross(prm%slip_direction(1:3,s1),prm%slip_direction(1:3,s2)))) .and. & any(dNeq0(math_cross(prm%slip_normal (1:3,s1),prm%slip_normal (1:3,s2))))) & prm%colinearSystem(s1) = s2 - enddo - enddo + end do + end do ini%rho_u_ed_pos_0 = pl%get_as1dFloat('rho_u_ed_pos_0', requiredSize=size(ini%N_sl)) ini%rho_u_ed_neg_0 = pl%get_as1dFloat('rho_u_ed_neg_0', requiredSize=size(ini%N_sl)) @@ -391,7 +391,7 @@ module function plastic_nonlocal_init() result(myPlasticity) if (prm%f_ed_mult < 0.0_pReal .or. prm%f_ed_mult > 1.0_pReal) & extmsg = trim(extmsg)//' f_ed_mult' - endif slipActive + end if slipActive !-------------------------------------------------------------------------------------------------- ! allocate state arrays @@ -506,7 +506,7 @@ module function plastic_nonlocal_init() result(myPlasticity) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(nonlocal)') - enddo + end do allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,& discretization_nIPs,discretization_Nelems), source=0.0_pReal) @@ -527,24 +527,24 @@ module function plastic_nonlocal_init() result(myPlasticity) do s = 1,param(ph)%sum_N_sl l = l + 1 iRhoU(s,t,ph) = l - enddo - enddo + end do + end do l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest do t = 1,4 do s = 1,param(ph)%sum_N_sl l = l + 1 iV(s,t,ph) = l - enddo - enddo + end do + end do do t = 1,2 do s = 1,param(ph)%sum_N_sl l = l + 1 iD(s,t,ph) = l - enddo - enddo + end do + end do if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) & error stop 'state indices not properly set (nonlocal)' - enddo + end do end function plastic_nonlocal_init @@ -625,7 +625,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) / log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl) else myInteractionMatrix = prm%h_sl_sl - endif + end if dst%tau_pass(:,en) = prm%mu * prm%b_sl & * sqrt(matmul(myInteractionMatrix,sum(abs(rho),2))) @@ -680,14 +680,14 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) connection_latticeConf(1:3,n) = 0.0_pReal rho_edg_delta_neighbor(:,n) = rho_edg_delta rho_scr_delta_neighbor(:,n) = rho_scr_delta - endif + end if else ! free surface -> use central values instead connection_latticeConf(1:3,n) = 0.0_pReal rho_edg_delta_neighbor(:,n) = rho_edg_delta rho_scr_delta_neighbor(:,n) = rho_scr_delta - endif - enddo + end if + end do neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor @@ -709,12 +709,12 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) - connection_latticeConf(1:3,neighbors(2)) rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) & - neighbor_rhoExcess(c,s,neighbors(2)) - enddo + end do invConnections = math_inv33(connections) if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error') rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences)) - enddo + end do ! ... plus gradient from deads ... rhoExcessGradient(1) = rhoExcessGradient(1) + sum(rho(s,imm_edg)) / FVsize @@ -731,8 +731,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) dst%tau_back(s,en) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) & * ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) & + rhoExcessGradient_over_rho(2)) - enddo - endif + end do + end if end associate @@ -790,8 +790,8 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & else tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s)) tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s)) - endif - enddo + end if + end do tauNS = tauNS + spread(dst%tau_back(:,en),2,4) tau = tau + dst%tau_back(:,en) @@ -807,12 +807,12 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & do t = 3,4 call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph) - enddo + end do else v(:,3:4) = spread(v(:,1),2,2) dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2) dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2) - endif + end if stt%v(:,en) = pack(v,.true.) @@ -833,7 +833,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & + prm%P_sl(i,j,s) & * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) - enddo + end do end associate @@ -899,7 +899,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) do s = 1,prm%sum_N_sl tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) +dst%tau_back(s,en) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal - enddo + end do dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau)) @@ -980,7 +980,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & if (timestep <= 0.0_pReal) then plasticState(ph)%dotState = 0.0_pReal return - endif + end if associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph)) @@ -1002,7 +1002,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & do s = 1,prm%sum_N_sl tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal - enddo + end do dLower = prm%minDipoleHeight dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) @@ -1032,7 +1032,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & rhoDotMultiplication(:,1:4) = spread( & (sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) & * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 - endif isBCC + end if isBCC forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) @@ -1062,7 +1062,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & - rhoDotSingle2DipoleGlide(:,2*c-1) & - rhoDotSingle2DipoleGlide(:,2*c) - enddo + end do ! athermal annihilation @@ -1103,13 +1103,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & if (debugConstitutive%extensive) then print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip print'(a)', '<< CONST >> enforcing cutback !!!' - endif + end if #endif plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) else dot%rho(:,en) = pack(rhoDot,.true.) dot%gamma(:,en) = sum(dot_gamma,2) - endif + end if end associate @@ -1217,11 +1217,11 @@ function rhoDotFlux(timestep,ph,en,ip,el) > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & ' at a timestep of ',timestep print*, '<< CONST >> enforcing cutback !!!' - endif + end if #endif rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! enforce cutback return - endif + end if !*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!! @@ -1255,7 +1255,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) Favg = 0.5_pReal * (my_F + neighbor_F) else ! if no neighbor, take my value as average Favg = my_F - endif + end if neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below @@ -1300,10 +1300,10 @@ function rhoDotFlux(timestep,ph,en,ip,el) rhoDotFlux(:,topp) = rhoDotFlux(:,topp) & + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to opposite signed mobile dislocation type - endif - enddo - enddo - endif; endif + end if + end do + end do + end if; end if !* FLUX FROM ME TO MY NEIGHBOR @@ -1330,20 +1330,20 @@ function rhoDotFlux(timestep,ph,en,ip,el) transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor transmissivity = 0.0_pReal - endif + end if lineLength = my_rhoSgl0(s,t) * v0(s,t) & * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & * sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point - endif - enddo - enddo - endif; endif + end if + end do + end do + end if; end if - enddo neighbors - endif + end do neighbors + end if end associate @@ -1433,7 +1433,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) mis%rotate(prm%slip_normal(1:3,s2)))) & * abs(math_inner(prm%slip_direction(1:3,s1), & mis%rotate(prm%slip_direction(1:3,s2)))) - enddo neighborSlipSystems + end do neighborSlipSystems my_compatibilitySum = 0.0_pReal belowThreshold = .true. @@ -1446,15 +1446,15 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,& my_compatibility(:,:,s1,n)) my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue - enddo + end do where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal - enddo mySlipSystems - endif + end do mySlipSystems + end if - enddo neighbors + end do neighbors compatibility(:,:,:,:,i,e) = my_compatibility @@ -1532,7 +1532,7 @@ module subroutine plastic_nonlocal_results(ph,group) 'passing stress for slip','Pa', prm%systems_sl) end select - enddo + end do end associate @@ -1581,7 +1581,7 @@ subroutine stateInit(ini,phase,Nentries) s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) meanDensity = meanDensity + densityBinning * geom(phase)%V_0(e) / totalVolume stt%rhoSglMobile(s,e) = densityBinning - enddo + end do else ! homogeneous distribution with noise do e = 1, Nentries do f = 1,size(ini%N_sl,1) @@ -1594,12 +1594,12 @@ subroutine stateInit(ini,phase,Nentries) stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2) stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2) - enddo + end do stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f) stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f) - enddo - enddo - endif + end do + end do + end if end associate @@ -1688,8 +1688,8 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2.0_pReal)) dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / lambda_P - endif - enddo + end if + end do end associate @@ -1760,8 +1760,8 @@ subroutine storeGeometry(ph) do ce = 1, size(material_homogenizationEntry,1) do co = 1, homogenization_maxNconstituents if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) - enddo - enddo + end do + end do end subroutine diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 51e347eb5..6e79968f3 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -95,8 +95,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) myPlasticity = plastic_active('phenopowerlaw') if(count(myPlasticity) == 0) return - print'(/,a)', ' <<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>' - print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>' + print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) phases => config_material%get('phase') @@ -105,7 +105,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) allocate(dotState(phases%length)) do ph = 1, phases%length - if(.not. myPlasticity(ph)) cycle + if (.not. myPlasticity(ph)) cycle associate(prm => param(ph), dot => dotState(ph), stt => state(ph)) @@ -129,7 +129,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) else prm%P_nS_pos = prm%P_sl prm%P_nS_neg = prm%P_sl - endif + end if prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'),phase_lattice(ph)) xi_0_sl = pl%get_as1dFloat('xi_0_sl', requiredSize=size(N_sl)) @@ -158,7 +158,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) xi_0_sl = emptyRealArray allocate(prm%xi_inf_sl,prm%h_int,source=emptyRealArray) allocate(prm%h_sl_sl(0,0)) - endif slipActive + end if slipActive !-------------------------------------------------------------------------------------------------- ! twin related parameters @@ -192,7 +192,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) xi_0_tw = emptyRealArray allocate(prm%gamma_char,source=emptyRealArray) allocate(prm%h_tw_tw(0,0)) - endif twinActive + end if twinActive !-------------------------------------------------------------------------------------------------- ! slip-twin related parameters @@ -206,7 +206,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0 allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0 prm%h_0_tw_sl = 0.0_pReal - endif slipAndTwinActive + end if slipAndTwinActive !-------------------------------------------------------------------------------------------------- ! output pararameters @@ -263,7 +263,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(phenopowerlaw)') - enddo + end do end function plastic_phenopowerlaw_init @@ -306,7 +306,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_sl_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) & + ddot_gamma_dtau_sl_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i) - enddo slipSystems + end do slipSystems call kinetics_tw(Mp,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) twinSystems: do i = 1, prm%sum_N_tw @@ -314,7 +314,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_tw(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i) - enddo twinSystems + end do twinSystems end associate @@ -397,7 +397,7 @@ module subroutine plastic_phenopowerlaw_results(ph,group) end select - enddo + end do end associate @@ -438,7 +438,7 @@ pure subroutine kinetics_sl(Mp,ph,en, & tau_sl_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) tau_sl_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)), & 0.0_pReal, prm%nonSchmidActive) - enddo + end do where(dNeq0(tau_sl_pos)) dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active @@ -460,14 +460,14 @@ pure subroutine kinetics_sl(Mp,ph,en, & else where ddot_gamma_dtau_sl_pos = 0.0_pReal end where - endif + end if if (present(ddot_gamma_dtau_sl_neg)) then where(dNeq0(dot_gamma_sl_neg)) ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg else where ddot_gamma_dtau_sl_neg = 0.0_pReal end where - endif + end if end associate @@ -517,7 +517,7 @@ pure subroutine kinetics_tw(Mp,ph,en,& else where ddot_gamma_dtau_tw = 0.0_pReal end where - endif + end if end associate diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 9a75263ca..4ea21d039 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -86,7 +86,7 @@ module subroutine thermal_init(phases) Nmembers - print'(/,a)', ' <<<+- phase:thermal init -+>>>' + print'(/,1x,a)', '<<<+- phase:thermal init -+>>>' allocate(current(phases%length)) diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 index cbdf5bb99..898b32706 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -36,8 +36,8 @@ module function dissipation_init(source_length) result(mySources) mySources = thermal_active('dissipation',source_length) if(count(mySources) == 0) return - print'(/,a)', ' <<<+- phase:thermal:dissipation init -+>>>' - print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:thermal:dissipation init -+>>>' + print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) phases => config_material%get('phase') @@ -45,11 +45,11 @@ module function dissipation_init(source_length) result(mySources) do ph = 1, phases%length phase => phases%get(ph) - if(count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1 + if (count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1 thermal => phase%get('thermal') sources => thermal%get('source') do so = 1, sources%length - if(mySources(so,ph)) then + if (mySources(so,ph)) then associate(prm => param(ph)) src => sources%get(so) @@ -58,9 +58,9 @@ module function dissipation_init(source_length) result(mySources) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) end associate - endif - enddo - enddo + end if + end do + end do end function dissipation_init diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index 46da127f3..f7a8296a8 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -43,8 +43,8 @@ module function externalheat_init(source_length) result(mySources) mySources = thermal_active('externalheat',source_length) if(count(mySources) == 0) return - print'(/,a)', ' <<<+- phase:thermal:externalheat init -+>>>' - print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:thermal:externalheat init -+>>>' + print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) phases => config_material%get('phase') @@ -53,11 +53,11 @@ module function externalheat_init(source_length) result(mySources) do ph = 1, phases%length phase => phases%get(ph) - if(count(mySources(:,ph)) == 0) cycle + if (count(mySources(:,ph)) == 0) cycle thermal => phase%get('thermal') sources => thermal%get('source') do so = 1, sources%length - if(mySources(so,ph)) then + if (mySources(so,ph)) then source_thermal_externalheat_offset(ph) = so associate(prm => param(ph)) src => sources%get(so) @@ -70,9 +70,9 @@ module function externalheat_init(source_length) result(mySources) Nmembers = count(material_phaseID == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) end associate - endif - enddo - enddo + end if + end do + end do end function externalheat_init @@ -125,7 +125,7 @@ module function externalheat_f_T(ph,en) result(f_T) f_T = prm%f_T(interval ) * (1.0_pReal - frac_time) + & prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... ! ...or extrapolate if outside of bounds - enddo + end do end associate end function externalheat_f_T diff --git a/src/prec.f90 b/src/prec.f90 index d6d161a94..da6760403 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -69,15 +69,15 @@ contains !-------------------------------------------------------------------------------------------------- subroutine prec_init - print'(/,a)', ' <<<+- prec init -+>>>' + print'(/,1x,a)', '<<<+- prec init -+>>>' - print'(a,i3)', ' Size of integer in bit: ',bit_size(0) - print'(a,i19)', ' Maximum value: ',huge(0) - print'(/,a,i3)', ' Size of float in bit: ',storage_size(0.0_pReal) - print'(a,e10.3)', ' Maximum value: ',huge(0.0_pReal) - print'(a,e10.3)', ' Minimum value: ',PREAL_MIN - print'(a,e10.3)', ' Epsilon value: ',PREAL_EPSILON - print'(a,i3)', ' Decimal precision: ',precision(0.0_pReal) + print'(/,a,i3)', ' integer size / bit: ',bit_size(0) + print'( a,i19)', ' maximum value: ',huge(0) + print'(/,a,i3)', ' float size / bit: ',storage_size(0.0_pReal) + print'( a,e10.3)', ' maximum value: ',huge(0.0_pReal) + print'( a,e10.3)', ' minimum value: ',PREAL_MIN + print'( a,e10.3)', ' epsilon value: ',PREAL_EPSILON + print'( a,i3)', ' decimal precision: ',precision(0.0_pReal) call selfTest diff --git a/src/results.f90 b/src/results.f90 index 064829658..7c2346095 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -70,10 +70,10 @@ subroutine results_init(restart) character(len=:), allocatable :: date - print'(/,a)', ' <<<+- results init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- results init -+>>>'; flush(IO_STDOUT) - print*, 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017' - print*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL + print'(/,1x,a)', 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017' + print'( 1x,a)', 'https://doi.org/10.1007/s40192-017-0084-5' if (.not. restart) then resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w') @@ -98,7 +98,7 @@ subroutine results_init(restart) call results_closeGroup(results_addGroup('setup')) call results_addAttribute('description','input data used to run the simulation','setup') call h5gmove_f(resultsFile,'tmp','setup/previous',hdferr) - endif + end if call results_closeJobFile @@ -222,7 +222,7 @@ subroutine results_addAttribute_str(attrLabel,attrValue,path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) else call HDF5_addAttribute(resultsFile,attrLabel, attrValue) - endif + end if end subroutine results_addAttribute_str @@ -241,7 +241,7 @@ subroutine results_addAttribute_int(attrLabel,attrValue,path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) else call HDF5_addAttribute(resultsFile,attrLabel, attrValue) - endif + end if end subroutine results_addAttribute_int @@ -260,7 +260,7 @@ subroutine results_addAttribute_real(attrLabel,attrValue,path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) else call HDF5_addAttribute(resultsFile,attrLabel, attrValue) - endif + end if end subroutine results_addAttribute_real @@ -279,7 +279,7 @@ subroutine results_addAttribute_str_array(attrLabel,attrValue,path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) else call HDF5_addAttribute(resultsFile,attrLabel, attrValue) - endif + end if end subroutine results_addAttribute_str_array @@ -298,7 +298,7 @@ subroutine results_addAttribute_int_array(attrLabel,attrValue,path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) else call HDF5_addAttribute(resultsFile,attrLabel, attrValue) - endif + end if end subroutine results_addAttribute_int_array @@ -317,7 +317,7 @@ subroutine results_addAttribute_real_array(attrLabel,attrValue,path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) else call HDF5_addAttribute(resultsFile,attrLabel, attrValue) - endif + end if end subroutine results_addAttribute_real_array @@ -390,7 +390,7 @@ subroutine results_writeVectorDataset_real(dataset,group,label,description,SIuni if (present(systems)) then if (size(systems)*size(dataset,2) == 0 ) return !ToDo: maybe also implement for other results_write (not sure about scalar) - endif + end if groupHandle = results_openGroup(group) call HDF5_write(dataset,groupHandle,label) @@ -422,7 +422,7 @@ subroutine results_writeTensorDataset_real(dataset,group,label,description,SIuni transposed_ = transposed else transposed_ = .true. - endif + end if groupHandle = results_openGroup(group) if(transposed_) then @@ -430,11 +430,11 @@ subroutine results_writeTensorDataset_real(dataset,group,label,description,SIuni allocate(dataset_transposed,mold=dataset) do i=1,size(dataset_transposed,3) dataset_transposed(:,:,i) = transpose(dataset(:,:,i)) - enddo + end do call HDF5_write(dataset_transposed,groupHandle,label) else call HDF5_write(dataset,groupHandle,label) - endif + end if call executionStamp(group//'/'//label,description,SIunit) call HDF5_closeGroup(groupHandle) @@ -456,7 +456,7 @@ subroutine results_writeVectorDataset_int(dataset,group,label,description,SIunit if (present(systems)) then if (size(systems)*size(dataset,2) == 0 ) return !ToDo: maybe also implement for other results_write (not sure about scalar) - endif + end if groupHandle = results_openGroup(group) call HDF5_write(dataset,groupHandle,label) @@ -542,16 +542,16 @@ subroutine results_mapping_phase(ID,entry,label) do co = 1, size(ID,1) do ce = 1, size(ID,2) entryOffset(ID(co,ce),worldrank) = entryOffset(ID(co,ce),worldrank) +1 - enddo - enddo + end do + end do call MPI_Allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,MPI_COMM_WORLD,ierr)! get offset at each process if(ierr /= 0) error stop 'MPI error' entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2) do co = 1, size(ID,1) do ce = 1, size(ID,2) entryGlobal(co,ce) = entry(co,ce) -1 + entryOffset(ID(co,ce),worldrank) - enddo - enddo + end do + end do #endif myShape = int([size(ID,1),writeSize(worldrank)], HSIZE_T) @@ -694,13 +694,13 @@ subroutine results_mapping_homogenization(ID,entry,label) entryOffset = 0 do ce = 1, size(ID,1) entryOffset(ID(ce),worldrank) = entryOffset(ID(ce),worldrank) +1 - enddo + end do call MPI_Allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,MPI_COMM_WORLD,ierr)! get offset at each process if(ierr /= 0) error stop 'MPI error' entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2) do ce = 1, size(ID,1) entryGlobal(ce) = entry(ce) -1 + entryOffset(ID(ce),worldrank) - enddo + end do #endif myShape = int([writeSize(worldrank)], HSIZE_T) diff --git a/src/rotations.f90 b/src/rotations.f90 index c64169cab..68f9273a3 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -103,10 +103,10 @@ contains !-------------------------------------------------------------------------------------------------- subroutine rotations_init - print'(/,a)', ' <<<+- rotations init -+>>>'; flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- rotations init -+>>>'; flush(IO_STDOUT) - print*, 'D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015' - print*, 'https://doi.org/10.1088/0965-0393/23/8/083501' + print'(/,1x,a)', 'D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015' + print'( 1x,a)', 'https://doi.org/10.1088/0965-0393/23/8/083501' call selfTest From 43ae4983ed56fbfc4c1e7aa3b031816280e83dad Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 15 Nov 2021 12:58:59 -0500 Subject: [PATCH 2/5] fixed misaligned/too short "writing..." statements --- src/grid/DAMASK_grid.f90 | 4 ++-- src/grid/spectral_utilities.f90 | 2 +- src/mesh/DAMASK_mesh.f90 | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 4e79bb6a8..5b3d6084b 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -317,7 +317,7 @@ program DAMASK_grid endif writeUndeformed: if (interface_restartInc < 1) then - print'(/,1x,a)', '... writing initial configuration to file ........................' + print'(/,1x,a)', '... writing initial configuration to file .................................' flush(IO_STDOUT) call CPFEM_results(0,0.0_pReal) endif writeUndeformed @@ -444,7 +444,7 @@ program DAMASK_grid call MPI_Allreduce(interface_SIGUSR1,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) if (ierr /= 0) error stop 'MPI error' if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then - print'(/,1x,a)', '... writing results to file ......................................' + print'(/,1x,a)', '... writing results to file ...............................................' flush(IO_STDOUT) call CPFEM_results(totalIncsCounter,t) endif diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index da8c3f05e..89ce9a314 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -1094,7 +1094,7 @@ subroutine utilities_saveReferenceStiffness fileUnit,ierr if (worldrank == 0) then - print'(1x,a)', 'writing reference stiffness data required for restart to file'; flush(IO_STDOUT) + print'(/,1x,a)', '... writing reference stiffness data required for restart to file .........'; flush(IO_STDOUT) open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',& status='replace',access='stream',action='write',iostat=ierr) if (ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref') diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 372cc26ca..14872bf15 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -237,7 +237,7 @@ program DAMASK_mesh write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file end if - print'(/,1x,a)', '... writing initial configuration to file ........................' + print'(/,1x,a)', '... writing initial configuration to file .................................' flush(IO_STDOUT) call CPFEM_results(0,0.0_pReal) @@ -321,7 +321,7 @@ program DAMASK_mesh end if; flush(IO_STDOUT) if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency - print'(/,1x,a)', '... writing results to file ......................................' + print'(/,1x,a)', '... writing results to file ...............................................' call FEM_mechanical_updateCoords call CPFEM_results(totalIncsCounter,time) end if From 262ff03b9ff9c964c596aaa2ee93ef5e48cb1f34 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Wed, 17 Nov 2021 09:11:30 +0100 Subject: [PATCH 3/5] more consistent formatting of loops --- src/homogenization_mechanical_RGC.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index b02e48498..71d1542c8 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -373,7 +373,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) do i=1,3; do j=1,3; do k=1,3; do l=1,3 smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l) - end do;enddo;enddo;enddo + end do;end do;end do;end do ! projecting the material tangent dPdF into the interface ! to obtain the Jacobian matrix contribution of dPdF end if @@ -394,7 +394,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) do i=1,3; do j=1,3; do k=1,3; do l=1,3 smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l) - end do;enddo;enddo;enddo + end do;end do;end do;end do end if end do end do From b2fd78b6396012bc50d2940e1d838ac41091c487 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 17 Nov 2021 13:45:27 +0100 Subject: [PATCH 4/5] [skip ci] updated version information after successful test of v3.0.0-alpha5-128-g5a901a42f --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 3999732f4..0f5f56c96 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-109-g58410709b +v3.0.0-alpha5-128-g5a901a42f From 52d0ba4a7444f27684067c9e25c1f290ce7a7fb5 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 17 Nov 2021 19:56:38 +0100 Subject: [PATCH 5/5] [skip ci] updated version information after successful test of v3.0.0-alpha5-135-ga3f74994b --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 0f5f56c96..3b4b0236b 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-128-g5a901a42f +v3.0.0-alpha5-135-ga3f74994b