From 59890750099ad650e0f7e310cea7bd74126c4fc6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 10 Apr 2019 13:23:57 +0200 Subject: [PATCH 01/12] working for grid_mech_spectral_basic --- src/CPFEM2.f90 | 2 +- src/grid_mech_spectral_basic.f90 | 70 ++++++++++++++------------------ 2 files changed, 32 insertions(+), 40 deletions(-) diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 13d7f06c4..aba6006b2 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -257,7 +257,7 @@ subroutine CPFEM_age() write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file' write(rankStr,'(a1,i0)')'_',worldrank - fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a') call HDF5_write(fileHandle,material_phase, 'recordedPhase') call HDF5_write(fileHandle,crystallite_F0, 'convergedF') diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90 index 99839e50f..2ed8bc683 100644 --- a/src/grid_mech_spectral_basic.f90 +++ b/src/grid_mech_spectral_basic.f90 @@ -7,6 +7,8 @@ module grid_mech_spectral_basic #include #include + use DAMASK_interface + use HDF5_utilities use PETScdmda use PETScsnes use prec, only: & @@ -114,8 +116,9 @@ subroutine grid_mech_spectral_basic_init PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & - F ! pointer to solution data + F ! pointer to solution data PetscInt, dimension(worldsize) :: localK + integer(HID_T) :: fileHandle integer :: fileUnit character(len=1024) :: rankStr @@ -174,19 +177,13 @@ subroutine grid_mech_spectral_basic_init restart: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('F_aim') - read(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - read(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - write(rankStr,'(a1,i0)')'_',worldrank - - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') elseif (restartInc == 0) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity @@ -203,10 +200,10 @@ subroutine grid_mech_spectral_basic_init restartRead: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) + fileUnit = IO_open_jobFile_binary('C_ref') read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead @@ -321,7 +318,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: F - integer :: fileUnit + integer(HID_T) :: fileHandle character(len=32) :: rankStr call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -331,29 +328,24 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? else !-------------------------------------------------------------------------------------------------- - ! restart information for spectral solver + ! restart information for spectral solver if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - - if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim','w') - write(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - write(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) - endif - + write(6,'(/,a)') ' writing converged results for restart';flush(6) + write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + + print*, trim(getSolverJobName())//trim(rankStr)//'.hdf5';flush(6) + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_write(fileHandle,C_minMaxAvg, 'C_minMaxAvg') + call HDF5_closeFile(fileHandle) endif call CPFEM_age ! age state and kinematics From 1852f580adf220b3bead9c2b13cd1e4f2ade4dce Mon Sep 17 00:00:00 2001 From: Vitesh Shah Date: Wed, 10 Apr 2019 18:19:10 +0200 Subject: [PATCH 02/12] grid_mech_FEM works --- src/grid_mech_FEM.f90 | 70 ++++++++++++++++++------------------------- 1 file changed, 29 insertions(+), 41 deletions(-) diff --git a/src/grid_mech_FEM.f90 b/src/grid_mech_FEM.f90 index e31d93637..adcbb626b 100644 --- a/src/grid_mech_FEM.f90 +++ b/src/grid_mech_FEM.f90 @@ -7,6 +7,8 @@ module grid_mech_FEM #include #include + use DAMASK_interface + use HDF5_utilities use PETScdmda use PETScsnes use prec, only: & @@ -115,6 +117,7 @@ subroutine grid_mech_FEM_init PetscErrorCode :: ierr integer(pInt) :: rank integer :: fileUnit + integer(HID_T) :: fileHandle character(len=1024) :: rankStr real(pReal), dimension(3,3,3,3) :: devNull PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -203,24 +206,16 @@ subroutine grid_mech_FEM_init restart: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('F_aim') - read(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - read(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - write(rankStr,'(a1,i0)')'_',worldrank - - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u'//trim(rankStr)) - read(fileUnit) u_current; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr)) - read(fileUnit) u_lastInc; close (fileUnit) - + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_read(fileHandle,u_current, 'u') + call HDF5_read(fileHandle,u_lastInc, 'u_lastInc') + elseif (restartInc == 0) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) @@ -238,10 +233,9 @@ subroutine grid_mech_FEM_init restartRead: if (restartInc > 0_pInt) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) endif restartRead end subroutine grid_mech_FEM_init @@ -348,6 +342,7 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat real(pReal), dimension(3,3), intent(in) :: & rotation_BC PetscErrorCode :: ierr + integer(HID_T) :: fileHandle integer :: fileUnit character(len=32) :: rankStr PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -364,30 +359,23 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat if (restartWrite) then ! QUESTION: where is this logical properly set? write(6,'(/,a)') ' writing converged results for restart' flush(6) + + write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim','w') - write(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - write(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_write(fileHandle,u_current, 'u') + call HDF5_write(fileHandle,u_lastInc, 'u_lastInc') + call HDF5_closeFile(fileHandle) endif - write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u'//trim(rankStr),'w') - write(fileUnit) u_current; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr),'w') - write(fileUnit) u_lastInc; close (fileUnit) - endif call CPFEM_age() ! age state and kinematics call utilities_updateIPcoords(F) From 8609c959dffae0aaad950ce01db697c472fdefc7 Mon Sep 17 00:00:00 2001 From: Vitesh Shah Date: Thu, 11 Apr 2019 15:51:34 +0200 Subject: [PATCH 03/12] grid_mech_polarisation works --- src/grid_mech_spectral_polarisation.f90 | 121 +++++++++++++++--------- 1 file changed, 77 insertions(+), 44 deletions(-) diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index aff4913b1..e87b8ec46 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -7,6 +7,8 @@ module grid_mech_spectral_polarisation #include #include + use DAMASK_interface + use HDF5_utilities use PETScdmda use PETScsnes use prec, only: & @@ -124,6 +126,7 @@ subroutine grid_mech_spectral_polarisation_init F, & ! specific (sub)pointer F_tau ! specific (sub)pointer PetscInt, dimension(worldsize) :: localK + integer(HID_T) :: fileHandle integer :: fileUnit character(len=1024) :: rankStr @@ -183,23 +186,34 @@ subroutine grid_mech_spectral_polarisation_init restart: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('F_aim') - read(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - read(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr)) - read(fileUnit) F_tau; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr)) - read(fileUnit) F_tau_lastInc; close (fileUnit) + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_read(fileHandle,F_tau, 'F_tau') + call HDF5_read(fileHandle,F_tau_lastInc, 'F_tau_lastInc') + + ! fileUnit = IO_open_jobFile_binary('F_aim') + ! read(fileUnit) F_aim; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_aim_lastInc') + ! read(fileUnit) F_aim_lastInc; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_aimDot') + ! read(fileUnit) F_aimDot; close(fileUnit) + + ! write(rankStr,'(a1,i0)')'_',worldrank + + ! fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) + ! read(fileUnit) F; close (fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) + ! read(fileUnit) F_lastInc; close (fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr)) + ! read(fileUnit) F_tau; close (fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr)) + ! read(fileUnit) F_tau_lastInc; close (fileUnit) elseif (restartInc == 0) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity @@ -218,12 +232,16 @@ subroutine grid_mech_spectral_polarisation_init restartRead: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_ref') - read(fileUnit) C_minMaxAvg; close(fileUnit) + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_read(fileHandle,C_minMaxAvg,'C_ref') + call HDF5_closeFile(fileHandle) + ! fileUnit = IO_open_jobFile_binary('C_volAvg') + ! read(fileUnit) C_volAvg; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') + ! read(fileUnit) C_volAvgLastInc; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('C_ref') + ! read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead call Utilities_updateGamma(C_minMaxAvg,.true.) @@ -348,6 +366,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa real(pReal), dimension(3,3) :: F_lambda33 integer :: fileUnit + integer(HID_T) :: fileHandle character(len=32) :: rankStr call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) @@ -361,31 +380,45 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - - if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim','w') - write(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - write(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) - endif + write(6,'(/,a)') ' writing converged results for restart';flush(6) write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w') - write(fileUnit) F_tau; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w') - write(fileUnit) F_tau_lastInc; close (fileUnit) + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + + if (worldrank == 0) then + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + ! fileUnit = IO_open_jobFile_binary('C_volAvg','w') + ! write(fileUnit) C_volAvg; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') + ! write(fileUnit) C_volAvgLastInc; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_aim','w') + ! write(fileUnit) F_aim; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') + ! write(fileUnit) F_aim_lastInc; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_aimDot','w') + call HDF5_closeFile(fileHandle) + !write(fileUnit) F_aimDot; close(fileUnit) + endif + + !write(rankStr,'(a1,i0)')'_',worldrank + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_write(fileHandle,F_tau, 'F_tau') + call HDF5_write(fileHandle,F_tau_lastInc, 'F_tau_lastInc') + call HDF5_closeFile(fileHandle) + +! fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') +! write(fileUnit) F; close (fileUnit) +! fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') +! write(fileUnit) F_lastInc; close (fileUnit) +! fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w') +! write(fileUnit) F_tau; close (fileUnit) +! fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w') +! write(fileUnit) F_tau_lastInc; close (fileUnit) endif call CPFEM_age ! age state and kinematics From a5b3fa565f1299add358bddcae65e7d70203f84c Mon Sep 17 00:00:00 2001 From: Vitesh Shah Date: Thu, 11 Apr 2019 18:36:28 +0200 Subject: [PATCH 04/12] Some incorrect variable names fixed --- src/grid_mech_spectral_polarisation.f90 | 68 ++++++++++++------------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index e87b8ec46..1f64be0a0 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -232,16 +232,16 @@ subroutine grid_mech_spectral_polarisation_init restartRead: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' - call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') - call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') - call HDF5_read(fileHandle,C_minMaxAvg,'C_ref') - call HDF5_closeFile(fileHandle) - ! fileUnit = IO_open_jobFile_binary('C_volAvg') - ! read(fileUnit) C_volAvg; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - ! read(fileUnit) C_volAvgLastInc; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('C_ref') - ! read(fileUnit) C_minMaxAvg; close(fileUnit) + ! call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + ! call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + ! call HDF5_read(fileHandle,C_minMaxAvg,'C_ref') + ! call HDF5_closeFile(fileHandle) + fileUnit = IO_open_jobFile_binary('C_volAvg') + read(fileUnit) C_volAvg; close(fileUnit) + fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') + read(fileUnit) C_volAvgLastInc; close(fileUnit) + fileUnit = IO_open_jobFile_binary('C_ref') + read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead call Utilities_updateGamma(C_minMaxAvg,.true.) @@ -384,41 +384,41 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa write(rankStr,'(a1,i0)')'_',worldrank fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') - + if (worldrank == 0) then - call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') - call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + ! call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + ! call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_write(fileHandle,F_aim, 'F_aim') call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') - call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') - ! fileUnit = IO_open_jobFile_binary('C_volAvg','w') - ! write(fileUnit) C_volAvg; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - ! write(fileUnit) C_volAvgLastInc; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_aim','w') - ! write(fileUnit) F_aim; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - ! write(fileUnit) F_aim_lastInc; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_aimDot','w') - call HDF5_closeFile(fileHandle) - !write(fileUnit) F_aimDot; close(fileUnit) + call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') + fileUnit = IO_open_jobFile_binary('C_volAvg','w') + write(fileUnit) C_volAvg; close(fileUnit) + fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') + write(fileUnit) C_volAvgLastInc; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_aim','w') + ! write(fileUnit) F_aim; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') + ! write(fileUnit) F_aim_lastInc; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('F_aimDot','w') + ! call HDF5_closeFile(fileHandle) + ! write(fileUnit) F_aimDot; close(fileUnit) endif - !write(rankStr,'(a1,i0)')'_',worldrank + write(rankStr,'(a1,i0)')'_',worldrank call HDF5_write(fileHandle,F, 'F') call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') call HDF5_write(fileHandle,F_tau, 'F_tau') call HDF5_write(fileHandle,F_tau_lastInc, 'F_tau_lastInc') call HDF5_closeFile(fileHandle) -! fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') -! write(fileUnit) F; close (fileUnit) -! fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') -! write(fileUnit) F_lastInc; close (fileUnit) -! fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w') -! write(fileUnit) F_tau; close (fileUnit) -! fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w') -! write(fileUnit) F_tau_lastInc; close (fileUnit) + !fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') + !write(fileUnit) F; close (fileUnit) + !fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') + !write(fileUnit) F_lastInc; close (fileUnit) + !fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w') + !write(fileUnit) F_tau; close (fileUnit) + !fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w') + !write(fileUnit) F_tau_lastInc; close (fileUnit) endif call CPFEM_age ! age state and kinematics From 057bee92afe0d055c58630bafd5c985b256a8414 Mon Sep 17 00:00:00 2001 From: Vitesh Shah Date: Thu, 11 Apr 2019 18:43:20 +0200 Subject: [PATCH 05/12] adding C_ref to HDF5 file doesnt work --- src/grid_mech_spectral_polarisation.f90 | 28 ++++++++++++------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index 1f64be0a0..c549cd591 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -232,14 +232,14 @@ subroutine grid_mech_spectral_polarisation_init restartRead: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' - ! call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') - ! call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') - ! call HDF5_read(fileHandle,C_minMaxAvg,'C_ref') - ! call HDF5_closeFile(fileHandle) - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + !call HDF5_read(fileHandle,C_minMaxAvg,'C_ref') + call HDF5_closeFile(fileHandle) + ! fileUnit = IO_open_jobFile_binary('C_volAvg') + ! read(fileUnit) C_volAvg; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') + ! read(fileUnit) C_volAvgLastInc; close(fileUnit) fileUnit = IO_open_jobFile_binary('C_ref') read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead @@ -386,15 +386,15 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') if (worldrank == 0) then - ! call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') - ! call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_write(fileHandle,F_aim, 'F_aim') call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('C_volAvg','w') + ! write(fileUnit) C_volAvg; close(fileUnit) + ! fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') + ! write(fileUnit) C_volAvgLastInc; close(fileUnit) ! fileUnit = IO_open_jobFile_binary('F_aim','w') ! write(fileUnit) F_aim; close(fileUnit) ! fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') From a59c5c30c90e134ea58cefdf2f429ab0f81cd16a Mon Sep 17 00:00:00 2001 From: Vitesh Shah Date: Thu, 11 Apr 2019 18:48:28 +0200 Subject: [PATCH 06/12] some cleanup --- src/grid_mech_spectral_polarisation.f90 | 41 ------------------------- 1 file changed, 41 deletions(-) diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index c549cd591..12183db0e 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -197,23 +197,6 @@ subroutine grid_mech_spectral_polarisation_init call HDF5_read(fileHandle,F_tau, 'F_tau') call HDF5_read(fileHandle,F_tau_lastInc, 'F_tau_lastInc') - ! fileUnit = IO_open_jobFile_binary('F_aim') - ! read(fileUnit) F_aim; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - ! read(fileUnit) F_aim_lastInc; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_aimDot') - ! read(fileUnit) F_aimDot; close(fileUnit) - - ! write(rankStr,'(a1,i0)')'_',worldrank - - ! fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - ! read(fileUnit) F; close (fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - ! read(fileUnit) F_lastInc; close (fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr)) - ! read(fileUnit) F_tau; close (fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr)) - ! read(fileUnit) F_tau_lastInc; close (fileUnit) elseif (restartInc == 0) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity @@ -234,12 +217,7 @@ subroutine grid_mech_spectral_polarisation_init write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') - !call HDF5_read(fileHandle,C_minMaxAvg,'C_ref') call HDF5_closeFile(fileHandle) - ! fileUnit = IO_open_jobFile_binary('C_volAvg') - ! read(fileUnit) C_volAvg; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - ! read(fileUnit) C_volAvgLastInc; close(fileUnit) fileUnit = IO_open_jobFile_binary('C_ref') read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead @@ -391,17 +369,6 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa call HDF5_write(fileHandle,F_aim, 'F_aim') call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') - ! fileUnit = IO_open_jobFile_binary('C_volAvg','w') - ! write(fileUnit) C_volAvg; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - ! write(fileUnit) C_volAvgLastInc; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_aim','w') - ! write(fileUnit) F_aim; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - ! write(fileUnit) F_aim_lastInc; close(fileUnit) - ! fileUnit = IO_open_jobFile_binary('F_aimDot','w') - ! call HDF5_closeFile(fileHandle) - ! write(fileUnit) F_aimDot; close(fileUnit) endif write(rankStr,'(a1,i0)')'_',worldrank @@ -411,14 +378,6 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa call HDF5_write(fileHandle,F_tau_lastInc, 'F_tau_lastInc') call HDF5_closeFile(fileHandle) - !fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - !write(fileUnit) F; close (fileUnit) - !fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - !write(fileUnit) F_lastInc; close (fileUnit) - !fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w') - !write(fileUnit) F_tau; close (fileUnit) - !fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w') - !write(fileUnit) F_tau_lastInc; close (fileUnit) endif call CPFEM_age ! age state and kinematics From 18f9deef1a7f59b9ebc26cacaafcc2368a91a528 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Apr 2019 15:27:58 +0200 Subject: [PATCH 07/12] fixed deadlock --- src/grid_mech_spectral_basic.f90 | 1 - src/grid_mech_spectral_polarisation.f90 | 19 +++++++------------ 2 files changed, 7 insertions(+), 13 deletions(-) diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90 index 2ed8bc683..d07c4c4d9 100644 --- a/src/grid_mech_spectral_basic.f90 +++ b/src/grid_mech_spectral_basic.f90 @@ -335,7 +335,6 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi write(rankStr,'(a1,i0)')'_',worldrank fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') - print*, trim(getSolverJobName())//trim(rankStr)//'.hdf5';flush(6) call HDF5_write(fileHandle,F_aim, 'F_aim') call HDF5_write(fileHandle,F_aim_lastInc,'F_aim_lastInc') call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index 12183db0e..568451971 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -107,8 +107,7 @@ subroutine grid_mech_spectral_polarisation_init use spectral_utilities, only: & utilities_constitutiveResponse, & utilities_updateGamma, & - utilities_updateIPcoords, & - wgt + utilities_updateIPcoords use mesh, only: & grid, & grid3 @@ -196,7 +195,6 @@ subroutine grid_mech_spectral_polarisation_init call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') call HDF5_read(fileHandle,F_tau, 'F_tau') call HDF5_read(fileHandle,F_tau_lastInc, 'F_tau_lastInc') - elseif (restartInc == 0) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity @@ -343,7 +341,6 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa integer :: i, j, k real(pReal), dimension(3,3) :: F_lambda33 - integer :: fileUnit integer(HID_T) :: fileHandle character(len=32) :: rankStr @@ -363,15 +360,13 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa write(rankStr,'(a1,i0)')'_',worldrank fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') - if (worldrank == 0) then - call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') - call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') - call HDF5_write(fileHandle,F_aim, 'F_aim') - call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') - call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') - endif + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') + - write(rankStr,'(a1,i0)')'_',worldrank call HDF5_write(fileHandle,F, 'F') call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') call HDF5_write(fileHandle,F_tau, 'F_tau') From 4793f964f8cc309f8126f6c4e948164012575bc0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Apr 2019 15:53:46 +0200 Subject: [PATCH 08/12] unified style --- src/grid_mech_spectral_basic.f90 | 33 +- src/grid_mech_spectral_polarisation.f90 | 830 ++++++++++++------------ 2 files changed, 430 insertions(+), 433 deletions(-) diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90 index d07c4c4d9..06f898224 100644 --- a/src/grid_mech_spectral_basic.f90 +++ b/src/grid_mech_spectral_basic.f90 @@ -27,8 +27,7 @@ module grid_mech_spectral_basic type(tSolutionParams), private :: params type, private :: tNumerics - logical :: & - update_gamma !< update gamma operator with current stiffness + logical :: update_gamma !< update gamma operator with current stiffness end type tNumerics type(tNumerics) :: num ! numerics parameters. Better name? @@ -42,8 +41,8 @@ module grid_mech_spectral_basic !-------------------------------------------------------------------------------------------------- ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: & - F_lastInc, & - Fdot + F_lastInc, & !< field of previous compatible deformation gradients + Fdot !< field of assumed rate of compatible deformation gradient !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. @@ -101,8 +100,7 @@ subroutine grid_mech_spectral_basic_init use spectral_utilities, only: & utilities_constitutiveResponse, & utilities_updateGamma, & - utilities_updateIPcoords, & - wgt + utilities_updateIPcoords use mesh, only: & grid, & grid3 @@ -179,6 +177,7 @@ subroutine grid_mech_spectral_basic_init write(rankStr,'(a1,i0)')'_',worldrank fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + call HDF5_read(fileHandle,F_aim, 'F_aim') call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') @@ -226,7 +225,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ terminallyIll implicit none - !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & @@ -248,8 +246,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (num%update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite) + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg,restartWrite) !-------------------------------------------------------------------------------------------------- ! set module wide available data @@ -329,7 +327,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi else !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? + if (restartWrite) then write(6,'(/,a)') ' writing converged results for restart';flush(6) write(rankStr,'(a1,i0)')'_',worldrank @@ -344,6 +342,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_write(fileHandle,C_minMaxAvg, 'C_minMaxAvg') + call HDF5_closeFile(fileHandle) endif @@ -390,7 +389,7 @@ end subroutine grid_mech_spectral_basic_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) +subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr) use numerics, only: & itmax, & itmin, & @@ -403,11 +402,11 @@ subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) implicit none SNES :: snes_local - PetscInt :: PETScIter - PetscReal :: & - xnorm, & ! not used - snorm, & ! not used - fnorm ! not used + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + devNull3 SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr @@ -506,7 +505,7 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) + call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & F,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index 568451971..1d470bed3 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -7,80 +7,79 @@ module grid_mech_spectral_polarisation #include #include - use DAMASK_interface - use HDF5_utilities - use PETScdmda - use PETScsnes - use prec, only: & - pReal - use math, only: & - math_I3 - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private + use DAMASK_interface + use HDF5_utilities + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private !-------------------------------------------------------------------------------------------------- ! derived types type(tSolutionParams), private :: params type, private :: tNumerics - logical :: & - update_gamma !< update gamma operator with current stiffness + logical :: update_gamma !< update gamma operator with current stiffness end type tNumerics type(tNumerics) :: num ! numerics parameters. Better name? !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: da - SNES, private :: snes - Vec, private :: solution_vec + DM, private :: da + SNES, private :: snes + Vec, private :: solution_vec !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: & - F_lastInc, & !< field of previous compatible deformation gradients - F_tau_lastInc, & !< field of previous incompatible deformation gradient - Fdot, & !< field of assumed rate of compatible deformation gradient - F_tauDot !< field of assumed rate of incopatible deformation gradient + real(pReal), private, dimension(:,:,:,:,:), allocatable :: & + F_lastInc, & !< field of previous compatible deformation gradients + F_tau_lastInc, & !< field of previous incompatible deformation gradient + Fdot, & !< field of assumed rate of compatible deformation gradient + F_tauDot !< field of assumed rate of incopatible deformation gradient !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & - F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient - F_aim = math_I3, & !< current prescribed deformation gradient - F_aim_lastInc = math_I3, & !< previous average deformation gradient - F_av = 0.0_pReal, & !< average incompatible def grad field - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - - character(len=1024), private :: incInfo !< time and increment information - real(pReal), private, dimension(3,3,3,3) :: & - C_volAvg = 0.0_pReal, & !< current volume average stiffness - C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness - C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness - C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness - S = 0.0_pReal, & !< current compliance (filled up with zeros) - C_scale = 0.0_pReal, & - S_scale = 0.0_pReal - - real(pReal), private :: & - err_BC, & !< deviation from stress BC - err_curl, & !< RMS of curl of F - err_div !< RMS of div of P + real(pReal), private, dimension(3,3) :: & + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient + F_aim = math_I3, & !< current prescribed deformation gradient + F_aim_lastInc = math_I3, & !< previous average deformation gradient + F_av = 0.0_pReal, & !< average incompatible def grad field + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - integer, private :: & - totalIter = 0 !< total iteration in current increment + character(len=1024), private :: incInfo !< time and increment information + real(pReal), private, dimension(3,3,3,3) :: & + C_volAvg = 0.0_pReal, & !< current volume average stiffness + C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness + C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness + S = 0.0_pReal, & !< current compliance (filled up with zeros) + C_scale = 0.0_pReal, & + S_scale = 0.0_pReal - public :: & - grid_mech_spectral_polarisation_init, & - grid_mech_spectral_polarisation_solution, & - grid_mech_spectral_polarisation_forward - private :: & - converged, & - formResidual + real(pReal), private :: & + err_BC, & !< deviation from stress BC + err_curl, & !< RMS of curl of F + err_div !< RMS of div of P + + integer, private :: & + totalIter = 0 !< total iteration in current increment + + public :: & + grid_mech_spectral_polarisation_init, & + grid_mech_spectral_polarisation_solution, & + grid_mech_spectral_polarisation_forward + private :: & + converged, & + formResidual contains @@ -88,141 +87,142 @@ contains !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine grid_mech_spectral_polarisation_init - use IO, only: & - IO_intOut, & - IO_error, & - IO_open_jobFile_binary - use FEsolving, only: & - restartInc - use config, only :& - config_numerics - use numerics, only: & - worldrank, & - worldsize, & - petsc_options - use homogenization, only: & - materialpoint_F0 - use DAMASK_interface, only: & - getSolverJobName - use spectral_utilities, only: & - utilities_constitutiveResponse, & - utilities_updateGamma, & - utilities_updateIPcoords - use mesh, only: & - grid, & - grid3 - use math, only: & - math_invSym3333 - - implicit none - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P - real(pReal), dimension(3,3) :: & - temp33_Real = 0.0_pReal - - PetscErrorCode :: ierr - PetscScalar, pointer, dimension(:,:,:,:) :: & - FandF_tau, & ! overall pointer to solution data - F, & ! specific (sub)pointer - F_tau ! specific (sub)pointer - PetscInt, dimension(worldsize) :: localK - integer(HID_T) :: fileHandle - integer :: fileUnit - character(len=1024) :: rankStr + use IO, only: & + IO_intOut, & + IO_error, & + IO_open_jobFile_binary + use FEsolving, only: & + restartInc + use config, only :& + config_numerics + use numerics, only: & + worldrank, & + worldsize, & + petsc_options + use homogenization, only: & + materialpoint_F0 + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + utilities_constitutiveResponse, & + utilities_updateGamma, & + utilities_updateIPcoords + use mesh, only: & + grid, & + grid3 + use math, only: & + math_invSym3333 + + implicit none + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal - write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>' - - write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' - write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' + PetscErrorCode :: ierr + PetscScalar, pointer, dimension(:,:,:,:) :: & + FandF_tau, & ! overall pointer to solution data + F, & ! specific (sub)pointer + F_tau ! specific (sub)pointer + PetscInt, dimension(worldsize) :: localK + integer(HID_T) :: fileHandle + integer :: fileUnit + character(len=1024) :: rankStr + + write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>' - num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0 + write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' + write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' + + num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0 !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) - CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) - CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 - call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) - call DMDACreate3d(PETSC_COMM_WORLD, & - DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary - DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - grid(1),grid(2),grid(3), & ! global grid - 1 , 1, worldsize, & - 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) - [grid(1)],[grid(2)],localK, & ! local grid - da,ierr) ! handle, error - CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da - call DMsetFromOptions(da,ierr); CHKERRQ(ierr) - call DMsetUp(da,ierr); CHKERRQ(ierr) - call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) - call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" - CHKERRQ(ierr) - call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments + call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) + localK = 0 + localK(worldrank+1) = grid3 + call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary + DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & ! global grid + 1 , 1, worldsize, & + 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) + [grid(1)],[grid(2)],localK, & ! local grid + da,ierr) ! handle, error + CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da + call DMsetFromOptions(da,ierr); CHKERRQ(ierr) + call DMsetUp(da,ierr); CHKERRQ(ierr) + call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) + call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" + CHKERRQ(ierr) + call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data - F => FandF_tau( 0: 8,:,:,:) - F_tau => FandF_tau( 9:17,:,:,:) + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data + F => FandF_tau( 0: 8,:,:,:) + F_tau => FandF_tau( 9:17,:,:,:) + + restart: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' + + write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_read(fileHandle,F_tau, 'F_tau') + call HDF5_read(fileHandle,F_tau_lastInc,'F_tau_lastInc') + + elseif (restartInc == 0) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_tau = 2.0_pReal*F + F_tau_lastInc = 2.0_pReal*F_lastInc + endif restart + + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) + call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + reshape(F,shape(F_lastInc)), & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer + + restartRead: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) - restart: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' - - write(rankStr,'(a1,i0)')'_',worldrank - fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') - - call HDF5_read(fileHandle,F_aim, 'F_aim') - call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') - call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') - call HDF5_read(fileHandle,F, 'F') - call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') - call HDF5_read(fileHandle,F_tau, 'F_tau') - call HDF5_read(fileHandle,F_tau_lastInc, 'F_tau_lastInc') - - elseif (restartInc == 0) then restart - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) - F_tau = 2.0_pReal*F - F_tau_lastInc = 2.0_pReal*F_lastInc - endif restart - - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 - reshape(F,shape(F_lastInc)), & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition - call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer - - restartRead: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' - call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') - call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') - call HDF5_closeFile(fileHandle) - fileUnit = IO_open_jobFile_binary('C_ref') - read(fileUnit) C_minMaxAvg; close(fileUnit) - endif restartRead - - call Utilities_updateGamma(C_minMaxAvg,.true.) - C_scale = C_minMaxAvg - S_scale = math_invSym3333(C_minMaxAvg) + fileUnit = IO_open_jobFile_binary('C_ref') + read(fileUnit) C_minMaxAvg; close(fileUnit) + endif restartRead + + call Utilities_updateGamma(C_minMaxAvg,.true.) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) end subroutine grid_mech_spectral_polarisation_init @@ -231,66 +231,65 @@ end subroutine grid_mech_spectral_polarisation_init !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) - use math, only: & - math_invSym3333 - use spectral_utilities, only: & - tBoundaryCondition, & - utilities_maskedCompliance, & - utilities_updateGamma - use FEsolving, only: & - restartWrite, & - terminallyIll + use math, only: & + math_invSym3333 + use spectral_utilities, only: & + tBoundaryCondition, & + utilities_maskedCompliance, & + utilities_updateGamma + use FEsolving, only: & + restartWrite, & + terminallyIll implicit none - !-------------------------------------------------------------------------------------------------- ! input data for solution - character(len=*), intent(in) :: & - incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC - type(tSolutionState) :: & - solution + character(len=*), intent(in) :: & + incInfoIn + real(pReal), intent(in) :: & + timeinc, & !< time increment of current solution + timeinc_old !< time increment of last successful increment + type(tBoundaryCondition), intent(in) :: & + stress_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(tSolutionState) :: & + solution !-------------------------------------------------------------------------------------------------- ! PETSc Data - PetscErrorCode :: ierr - SNESConvergedReason :: reason - - incInfo = incInfoIn + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (num%update_gamma) then - call utilities_updateGamma(C_minMaxAvg,restartWrite) - C_scale = C_minMaxAvg - S_scale = math_invSym3333(C_minMaxAvg) - endif + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + if (num%update_gamma) then + call utilities_updateGamma(C_minMaxAvg,restartWrite) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) + endif !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old + params%stress_mask = stress_BC%maskFloat + params%stress_BC = stress_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) + call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) - - solution%converged = reason > 0 - solution%iterationsNeeded = totalIter - solution%termIll = terminallyIll - terminallyIll = .false. + call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) + + solution%converged = reason > 0 + solution%iterationsNeeded = totalIter + solution%termIll = terminallyIll + terminallyIll = .false. end function grid_mech_spectral_polarisation_solution @@ -324,13 +323,12 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa use FEsolving, only: & restartWrite - implicit none logical, intent(in) :: & guess real(pReal), intent(in) :: & timeinc_old, & timeinc, & - loadCaseTime !< remaining time of current load case + loadCaseTime !< remaining time of current load case type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC @@ -349,30 +347,29 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa F_tau => FandF_tau( 9:17,:,:,:) if (cutBack) then - C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? - C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? + C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? + C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? else !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? + if (restartWrite) then write(6,'(/,a)') ' writing converged results for restart';flush(6) write(rankStr,'(a1,i0)')'_',worldrank fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_write(fileHandle,F_tau, 'F_tau') + call HDF5_write(fileHandle,F_tau_lastInc, 'F_tau_lastInc') + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') - call HDF5_write(fileHandle,F_aim, 'F_aim') - call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') - call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') - - call HDF5_write(fileHandle,F, 'F') - call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') - call HDF5_write(fileHandle,F_tau, 'F_tau') - call HDF5_write(fileHandle,F_tau_lastInc, 'F_tau_lastInc') call HDF5_closeFile(fileHandle) - endif call CPFEM_age ! age state and kinematics @@ -386,13 +383,13 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F F_aimDot = & F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed F_aimDot = & F_aimDot + deformation_BC%maskFloat * deformation_BC%values - elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed F_aimDot = & F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime endif @@ -404,33 +401,33 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa F_tauDot = utilities_calculateRate(guess, & F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & math_rotate_backward33(F_aimDot,rotation_BC)) - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward - F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau forward - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward + F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau forward + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent endif !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc - F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average math_rotate_backward33(F_aim,rotation_BC)),& [9,grid(1),grid(2),grid3]) - if (guess) then - F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & - [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition - else - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) - F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) - F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, & - math_mul3333xx33(C_scale,& - matmul(transpose(F_lambda33),& - F_lambda33)-math_I3))*0.5_pReal)& - + math_I3 - F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) - enddo; enddo; enddo - endif - - call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) + if (guess) then + F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & + [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + else + do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) + F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, & + math_mul3333xx33(C_scale,& + matmul(transpose(F_lambda33),& + F_lambda33)-math_I3))*0.5_pReal)& + + math_I3 + F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) + enddo; enddo; enddo + endif + + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) end subroutine grid_mech_spectral_polarisation_forward @@ -438,61 +435,61 @@ end subroutine grid_mech_spectral_polarisation_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) - use numerics, only: & - itmax, & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_curl_tolRel, & - err_curl_tolAbs, & - err_stress_tolRel, & - err_stress_tolAbs - use FEsolving, only: & - terminallyIll - - implicit none - SNES :: snes_local - PetscInt :: PETScIter - PetscReal :: & - xnorm, & ! not used - snorm, & ! not used - fnorm ! not used - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal) :: & - curlTol, & - divTol, & - BCTol - - curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) - divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) - BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) - - if ((totalIter >= itmin .and. & - all([ err_div /divTol, & - err_curl/curlTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then - reason = 1 - elseif (totalIter >= itmax) then - reason = -1 - else - reason = 0 - endif +subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_curl_tolRel, & + err_curl_tolAbs, & + err_stress_tolRel, & + err_stress_tolAbs + use FEsolving, only: & + terminallyIll + + implicit none + SNES :: snes_local + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + devNull3 + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal) :: & + curlTol, & + divTol, & + BCTol + + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) + + if ((totalIter >= itmin .and. & + all([ err_div /divTol, & + err_curl/curlTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then + reason = -1 + else + reason = 0 + endif !-------------------------------------------------------------------------------------------------- ! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' - write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' - write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine converged @@ -502,144 +499,145 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, FandF_tau, & residuum, dummy,ierr) - use numerics, only: & - itmax, & - itmin, & - polarAlpha, & - polarBeta - use mesh, only: & - grid, & - grid3 - use math, only: & - math_rotate_forward33, & - math_rotate_backward33, & - math_mul3333xx33, & - math_invSym3333 - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation - use spectral_utilities, only: & - wgt, & - tensorField_real, & - utilities_FFTtensorForward, & - utilities_fourierGammaConvolution, & - utilities_FFTtensorBackward, & - utilities_constitutiveResponse, & - utilities_divergenceRMS, & - utilities_curlRMS - use IO, only: & - IO_intOut - use homogenization, only: & - materialpoint_dPdF - use FEsolving, only: & - terminallyIll + use numerics, only: & + itmax, & + itmin, & + polarAlpha, & + polarBeta + use mesh, only: & + grid, & + grid3 + use math, only: & + math_rotate_forward33, & + math_rotate_backward33, & + math_mul3333xx33, & + math_invSym3333 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use spectral_utilities, only: & + wgt, & + tensorField_real, & + utilities_FFTtensorForward, & + utilities_fourierGammaConvolution, & + utilities_FFTtensorBackward, & + utilities_constitutiveResponse, & + utilities_divergenceRMS, & + utilities_curlRMS + use IO, only: & + IO_intOut + use homogenization, only: & + materialpoint_dPdF + use FEsolving, only: & + terminallyIll - implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) - PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), & - target, intent(in) :: FandF_tau - PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& - target, intent(out) :: residuum !< residuum field - PetscScalar, pointer, dimension(:,:,:,:,:) :: & - F, & - F_tau, & - residual_F, & - residual_F_tau - PetscInt :: & - PETScIter, & - nfuncs - PetscObject :: dummy - PetscErrorCode :: ierr - integer :: & - i, j, k, e + implicit none + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) + PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), & + target, intent(in) :: FandF_tau + PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& + target, intent(out) :: residuum !< residuum field + PetscScalar, pointer, dimension(:,:,:,:,:) :: & + F, & + F_tau, & + residual_F, & + residual_F_tau + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + integer :: & + i, j, k, e - F => FandF_tau(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) - F_tau => FandF_tau(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) - residual_F => residuum(1:3,1:3,1,& - X_RANGE, Y_RANGE, Z_RANGE) - residual_F_tau => residuum(1:3,1:3,2,& - X_RANGE, Y_RANGE, Z_RANGE) + F => FandF_tau(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_tau => FandF_tau(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => residuum(1:3,1:3,1,& + X_RANGE, Y_RANGE, Z_RANGE) + residual_F_tau => residuum(1:3,1:3,2,& + X_RANGE, Y_RANGE, Z_RANGE) - F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt - call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt + call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) - if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment !-------------------------------------------------------------------------------------------------- ! begin of new iteration - newIteration: if (totalIter <= PETScIter) then - totalIter = totalIter + 1 - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & - trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax - if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) - flush(6) - endif newIteration + newIteration: if (totalIter <= PETScIter) then + totalIter = totalIter + 1 + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim =', transpose(F_aim) + flush(6) + endif newIteration !-------------------------------------------------------------------------------------------------- ! - tensorField_real = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) - tensorField_real(1:3,1:3,i,j,k) = & - polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& - polarAlpha*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 + tensorField_real = 0.0_pReal + do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + tensorField_real(1:3,1:3,i,j,k) = & + polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + polarAlpha*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 !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space - call utilities_FFTtensorForward - call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) - call utilities_FFTtensorBackward + call utilities_FFTtensorForward + call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) - P_av,C_volAvg,C_minMaxAvg, & - F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) + P_av,C_volAvg,C_minMaxAvg, & + F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & - -math_rotate_forward33(F_av,params%rotation_BC)) + & - params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc + err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & + -math_rotate_forward33(F_av,params%rotation_BC)) + & + params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc ! calculate divergence - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise - call utilities_FFTtensorForward - err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise + call utilities_FFTtensorForward + err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress + !-------------------------------------------------------------------------------------------------- ! constructing residual - e = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) - e = e + 1 - residual_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & - residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & - math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & - + residual_F_tau(1:3,1:3,i,j,k) - enddo; enddo; enddo + e = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + e = e + 1 + residual_F(1:3,1:3,i,j,k) = & + math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & + residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + + residual_F_tau(1:3,1:3,i,j,k) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! calculating curl - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F - call utilities_FFTtensorForward - err_curl = Utilities_curlRMS() + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + call utilities_FFTtensorForward + err_curl = Utilities_curlRMS() end subroutine formResidual From 19dfefc2f48ed34255e1595a7814441af4dea7c3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Apr 2019 16:16:14 +0200 Subject: [PATCH 09/12] fixed deadlock in grid_mech_FEM --- src/grid_mech_FEM.f90 | 58 ++++++++++++------------- src/grid_mech_spectral_basic.f90 | 9 ++-- src/grid_mech_spectral_polarisation.f90 | 7 +-- 3 files changed, 33 insertions(+), 41 deletions(-) diff --git a/src/grid_mech_FEM.f90 b/src/grid_mech_FEM.f90 index adcbb626b..8ec42d918 100644 --- a/src/grid_mech_FEM.f90 +++ b/src/grid_mech_FEM.f90 @@ -208,6 +208,7 @@ subroutine grid_mech_FEM_init write(rankStr,'(a1,i0)')'_',worldrank fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + call HDF5_read(fileHandle,F_aim, 'F_aim') call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') @@ -221,8 +222,8 @@ subroutine grid_mech_FEM_init F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restart materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateIPcoords(F) - call Utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_updateIPcoords(F) + call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F 0.0_pReal, & ! time increment math_I3) ! no rotation of boundary condition @@ -255,7 +256,6 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation terminallyIll implicit none - !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & @@ -268,7 +268,6 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation real(pReal), dimension(3,3), intent(in) :: rotation_BC type(tSolutionState) :: & solution - !-------------------------------------------------------------------------------------------------- ! PETSc Data PetscErrorCode :: ierr @@ -278,7 +277,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) !-------------------------------------------------------------------------------------------------- ! set module wide available data params%stress_mask = stress_BC%maskFloat @@ -357,27 +356,26 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) + write(6,'(/,a)') ' writing converged results for restart';flush(6) write(rankStr,'(a1,i0)')'_',worldrank fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') - if (worldrank == 0) then - call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') - call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') - call HDF5_write(fileHandle,F_aim, 'F_aim') - call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') - call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') - call HDF5_write(fileHandle,F, 'F') - call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') - call HDF5_write(fileHandle,u_current, 'u') - call HDF5_write(fileHandle,u_lastInc, 'u_lastInc') - call HDF5_closeFile(fileHandle) - endif + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_write(fileHandle,u_current, 'u') + call HDF5_write(fileHandle,u_lastInc, 'u_lastInc') + + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + + call HDF5_closeFile(fileHandle) endif - call CPFEM_age() ! age state and kinematics + call CPFEM_age ! age state and kinematics call utilities_updateIPcoords(F) C_volAvgLastInc = C_volAvg @@ -429,9 +427,9 @@ end subroutine grid_mech_FEM_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) -use mesh -use spectral_utilities +subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr) + use mesh + use spectral_utilities use numerics, only: & itmax, & itmin, & @@ -444,11 +442,11 @@ use spectral_utilities implicit none SNES :: snes_local - PetscInt :: PETScIter - PetscReal :: & - xnorm, & ! not used - snorm, & ! not used - fnorm + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + devNull3 SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr @@ -461,7 +459,6 @@ use spectral_utilities divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) - if ((totalIter >= itmin .and. & all([ err_div/divTol, & err_BC /BCTol ] < 1.0_pReal)) & @@ -489,7 +486,8 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- !> @brief forms the residual vector !-------------------------------------------------------------------------------------------------- -subroutine formResidual(da_local,x_local,f_local,dummy,ierr) +subroutine formResidual(da_local,x_local, & + f_local,dummy,ierr) use numerics, only: & itmax, & itmin diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90 index 06f898224..91dcd03b4 100644 --- a/src/grid_mech_spectral_basic.f90 +++ b/src/grid_mech_spectral_basic.f90 @@ -117,7 +117,6 @@ subroutine grid_mech_spectral_basic_init F ! pointer to solution data PetscInt, dimension(worldsize) :: localK integer(HID_T) :: fileHandle - integer :: fileUnit character(len=1024) :: rankStr write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>' @@ -164,7 +163,7 @@ subroutine grid_mech_spectral_basic_init call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "converged" + call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" CHKERRQ(ierr) call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments @@ -203,11 +202,9 @@ subroutine grid_mech_spectral_basic_init call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_closeFile(fileHandle) - fileUnit = IO_open_jobFile_binary('C_ref') - read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead - call Utilities_updateGamma(C_minMaxAvg,.true.) + call utilities_updateGamma(C_minMaxAvg,.true.) end subroutine grid_mech_spectral_basic_init @@ -442,7 +439,7 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the basic residual vector +!> @brief forms the residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, F, & residuum, dummy, ierr) diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index 1d470bed3..a7bf75209 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -125,7 +125,6 @@ subroutine grid_mech_spectral_polarisation_init F_tau ! specific (sub)pointer PetscInt, dimension(worldsize) :: localK integer(HID_T) :: fileHandle - integer :: fileUnit character(len=1024) :: rankStr write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>' @@ -216,11 +215,9 @@ subroutine grid_mech_spectral_polarisation_init call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_closeFile(fileHandle) - fileUnit = IO_open_jobFile_binary('C_ref') - read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead - call Utilities_updateGamma(C_minMaxAvg,.true.) + call utilities_updateGamma(C_minMaxAvg,.true.) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) @@ -495,7 +492,7 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the polarisation residual vector +!> @brief forms the residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, FandF_tau, & residuum, dummy,ierr) From 8f58f1348ea91d6b55190ee607b1a4c9dc34f357 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Apr 2019 16:19:30 +0200 Subject: [PATCH 10/12] _pInt not needed --- src/grid_mech_FEM.f90 | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/src/grid_mech_FEM.f90 b/src/grid_mech_FEM.f90 index 8ec42d918..70332a040 100644 --- a/src/grid_mech_FEM.f90 +++ b/src/grid_mech_FEM.f90 @@ -12,7 +12,6 @@ module grid_mech_FEM use PETScdmda use PETScsnes use prec, only: & - pInt, & pReal use math, only: & math_I3 @@ -61,8 +60,8 @@ module grid_mech_FEM real(pReal), private :: & err_BC !< deviation from stress BC - integer(pInt), private :: & - totalIter = 0_pInt !< total iteration in current increment + integer, private :: & + totalIter = 0 !< total iteration in current increment public :: & grid_mech_FEM_init, & @@ -115,8 +114,7 @@ subroutine grid_mech_FEM_init 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) PetscErrorCode :: ierr - integer(pInt) :: rank - integer :: fileUnit + integer :: rank integer(HID_T) :: fileHandle character(len=1024) :: rankStr real(pReal), dimension(3,3,3,3) :: devNull @@ -232,7 +230,7 @@ subroutine grid_mech_FEM_init call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) CHKERRQ(ierr) - restartRead: if (restartInc > 0_pInt) then + restartRead: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') @@ -342,7 +340,6 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat rotation_BC PetscErrorCode :: ierr integer(HID_T) :: fileHandle - integer :: fileUnit character(len=32) :: rankStr PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc @@ -530,12 +527,12 @@ subroutine formResidual(da_local,x_local, & call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) - if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment !-------------------------------------------------------------------------------------------------- ! begin of new iteration newIteration: if (totalIter <= PETScIter) then - totalIter = totalIter + 1_pInt + totalIter = totalIter + 1 write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & From d0a4cbf8d5804d4e6c093970e924f206fbfabb79 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Apr 2019 16:30:54 +0200 Subject: [PATCH 11/12] one implicit none is enough --- src/grid_damage_spectral.f90 | 4 - src/grid_mech_FEM.f90 | 962 ++++++++++++------------ src/grid_mech_spectral_basic.f90 | 5 - src/grid_mech_spectral_polarisation.f90 | 4 - src/grid_thermal_spectral.f90 | 4 - 5 files changed, 478 insertions(+), 501 deletions(-) diff --git a/src/grid_damage_spectral.f90 b/src/grid_damage_spectral.f90 index c5e9a254b..3ce37c5ff 100644 --- a/src/grid_damage_spectral.f90 +++ b/src/grid_damage_spectral.f90 @@ -64,7 +64,6 @@ subroutine grid_damage_spectral_init worldsize, & petsc_options - implicit none PetscInt, dimension(worldsize) :: localK integer :: i, j, k, cell DM :: damage_grid @@ -164,7 +163,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result( use damage_nonlocal, only: & damage_nonlocal_putNonLocalDamage - implicit none real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment @@ -236,7 +234,6 @@ subroutine grid_damage_spectral_forward damage_nonlocal_getDiffusion33, & damage_nonlocal_getMobility - implicit none integer :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal @@ -301,7 +298,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) damage_nonlocal_getDiffusion33, & damage_nonlocal_getMobility - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & diff --git a/src/grid_mech_FEM.f90 b/src/grid_mech_FEM.f90 index 70332a040..2feec76df 100644 --- a/src/grid_mech_FEM.f90 +++ b/src/grid_mech_FEM.f90 @@ -7,66 +7,66 @@ module grid_mech_FEM #include #include - use DAMASK_interface - use HDF5_utilities - use PETScdmda - use PETScsnes - use prec, only: & - pReal - use math, only: & - math_I3 - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private + use DAMASK_interface + use HDF5_utilities + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private !-------------------------------------------------------------------------------------------------- ! derived types - type(tSolutionParams), private :: params + type(tSolutionParams), private :: params !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: mech_grid - SNES, private :: mech_snes - Vec, private :: solution_current, solution_lastInc, solution_rate + DM, private :: mech_grid + SNES, private :: mech_snes + Vec, private :: solution_current, solution_lastInc, solution_rate !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc - real(pReal), private :: detJ - real(pReal), private, dimension(3) :: delta - real(pReal), private, dimension(3,8) :: BMat - real(pReal), private, dimension(8,8) :: HGMat - PetscInt, private :: xstart,ystart,zstart,xend,yend,zend + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc + real(pReal), private :: detJ + real(pReal), private, dimension(3) :: delta + real(pReal), private, dimension(3,8) :: BMat + real(pReal), private, dimension(8,8) :: HGMat + PetscInt, private :: xstart,ystart,zstart,xend,yend,zend !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & - F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient - F_aim = math_I3, & !< current prescribed deformation gradient - F_aim_lastIter = math_I3, & - F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - - character(len=1024), private :: incInfo !< time and increment information - - real(pReal), private, dimension(3,3,3,3) :: & - C_volAvg = 0.0_pReal, & !< current volume average stiffness - C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness - S = 0.0_pReal !< current compliance (filled up with zeros) - - real(pReal), private :: & - err_BC !< deviation from stress BC - - integer, private :: & - totalIter = 0 !< total iteration in current increment - - public :: & - grid_mech_FEM_init, & - grid_mech_FEM_solution, & - grid_mech_FEM_forward + real(pReal), private, dimension(3,3) :: & + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient + F_aim = math_I3, & !< current prescribed deformation gradient + F_aim_lastIter = math_I3, & + F_aim_lastInc = math_I3, & !< previous average deformation gradient + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress + + character(len=1024), private :: incInfo !< time and increment information + + real(pReal), private, dimension(3,3,3,3) :: & + C_volAvg = 0.0_pReal, & !< current volume average stiffness + C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness + S = 0.0_pReal !< current compliance (filled up with zeros) + + real(pReal), private :: & + err_BC !< deviation from stress BC + + integer, private :: & + totalIter = 0 !< total iteration in current increment + + public :: & + grid_mech_FEM_init, & + grid_mech_FEM_solution, & + grid_mech_FEM_forward contains @@ -78,164 +78,163 @@ subroutine grid_mech_FEM_init IO_intOut, & IO_error, & IO_open_jobFile_binary - use FEsolving, only: & - restartInc - use numerics, only: & - worldrank, & - worldsize, & - petsc_options - use homogenization, only: & - materialpoint_F0 - use DAMASK_interface, only: & - getSolverJobName - use spectral_utilities, only: & - utilities_constitutiveResponse, & - utilities_updateIPcoords, & - wgt - use mesh, only: & - geomSize, & - grid, & - grid3 - use math, only: & - math_invSym3333 - - implicit none - real(pReal) :: HGCoeff = 0e-2_pReal - PetscInt, dimension(:), allocatable :: localK - real(pReal), dimension(3,3) :: & - temp33_Real = 0.0_pReal - real(pReal), dimension(4,8) :: & - HGcomp = reshape([ 1.0_pReal, 1.0_pReal, 1.0_pReal,-1.0_pReal, & - 1.0_pReal,-1.0_pReal,-1.0_pReal, 1.0_pReal, & - -1.0_pReal, 1.0_pReal,-1.0_pReal, 1.0_pReal, & - -1.0_pReal,-1.0_pReal, 1.0_pReal,-1.0_pReal, & - -1.0_pReal,-1.0_pReal, 1.0_pReal, 1.0_pReal, & - -1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, & - 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & - 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) - PetscErrorCode :: ierr - integer :: rank - integer(HID_T) :: fileHandle - character(len=1024) :: rankStr - real(pReal), dimension(3,3,3,3) :: devNull - PetscScalar, pointer, dimension(:,:,:,:) :: & - u_current,u_lastInc - - write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>' + use FEsolving, only: & + restartInc + use numerics, only: & + worldrank, & + worldsize, & + petsc_options + use homogenization, only: & + materialpoint_F0 + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + utilities_constitutiveResponse, & + utilities_updateIPcoords, & + wgt + use mesh, only: & + geomSize, & + grid, & + grid3 + use math, only: & + math_invSym3333 + + real(pReal) :: HGCoeff = 0e-2_pReal + PetscInt, dimension(:), allocatable :: localK + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal + real(pReal), dimension(4,8) :: & + HGcomp = reshape([ 1.0_pReal, 1.0_pReal, 1.0_pReal,-1.0_pReal, & + 1.0_pReal,-1.0_pReal,-1.0_pReal, 1.0_pReal, & + -1.0_pReal, 1.0_pReal,-1.0_pReal, 1.0_pReal, & + -1.0_pReal,-1.0_pReal, 1.0_pReal,-1.0_pReal, & + -1.0_pReal,-1.0_pReal, 1.0_pReal, 1.0_pReal, & + -1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, & + 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & + 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) + PetscErrorCode :: ierr + integer :: rank + integer(HID_T) :: fileHandle + character(len=1024) :: rankStr + real(pReal), dimension(3,3,3,3) :: devNull + PetscScalar, pointer, dimension(:,:,:,:) :: & + u_current,u_lastInc + + write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>' !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & - &-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) - CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) - CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & + &-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 - do rank = 1, worldsize - call MPI_Bcast(localK(rank),1,MPI_INTEGER,rank-1,PETSC_COMM_WORLD,ierr) - enddo - call DMDACreate3d(PETSC_COMM_WORLD, & - DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, & - DMDA_STENCIL_BOX, & - grid(1),grid(2),grid(3), & - 1, 1, worldsize, & - 3, 1, & - [grid(1)],[grid(2)],localK, & - mech_grid,ierr) - CHKERRQ(ierr) - call DMDASetUniformCoordinates(mech_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr) - CHKERRQ(ierr) - call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(ierr) - call DMsetFromOptions(mech_grid,ierr); CHKERRQ(ierr) - call DMsetUp(mech_grid,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(mech_grid,solution_current,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(mech_grid,solution_lastInc,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(mech_grid,solution_rate ,ierr); CHKERRQ(ierr) - call DMSNESSetFunctionLocal(mech_grid,formResidual,PETSC_NULL_SNES,ierr) - CHKERRQ(ierr) - call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) - CHKERRQ(ierr) - call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) - CHKERRQ(ierr) ! specify custom convergence check function "_converged" - call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures - call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments + call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 + do rank = 1, worldsize + call MPI_Bcast(localK(rank),1,MPI_INTEGER,rank-1,PETSC_COMM_WORLD,ierr) + enddo + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, & + DMDA_STENCIL_BOX, & + grid(1),grid(2),grid(3), & + 1, 1, worldsize, & + 3, 1, & + [grid(1)],[grid(2)],localK, & + mech_grid,ierr) + CHKERRQ(ierr) + call DMDASetUniformCoordinates(mech_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr) + CHKERRQ(ierr) + call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(ierr) + call DMsetFromOptions(mech_grid,ierr); CHKERRQ(ierr) + call DMsetUp(mech_grid,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(mech_grid,solution_current,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(mech_grid,solution_lastInc,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(mech_grid,solution_rate ,ierr); CHKERRQ(ierr) + call DMSNESSetFunctionLocal(mech_grid,formResidual,PETSC_NULL_SNES,ierr) + CHKERRQ(ierr) + call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) + CHKERRQ(ierr) + call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" + CHKERRQ(ierr) + call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures + call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments !-------------------------------------------------------------------------------------------------- ! init fields - call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr) - call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr) - call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) - - call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent - CHKERRQ(ierr) - xend = xstart+xend-1 - yend = ystart+yend-1 - zend = zstart+zend-1 - delta = geomSize/real(grid,pReal) ! grid spacing - detJ = product(delta) ! cell volume - - BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & - 1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & - -1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & - 1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & - -1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & - 1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & - -1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), & - 1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix - - HGMat = matmul(transpose(HGcomp),HGcomp) & - * HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix + call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr) + call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr) + call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) + + call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent + CHKERRQ(ierr) + xend = xstart+xend-1 + yend = ystart+yend-1 + zend = zstart+zend-1 + delta = geomSize/real(grid,pReal) ! grid spacing + detJ = product(delta) ! cell volume + + BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & + 1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & + -1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & + 1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & + -1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & + 1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & + -1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), & + 1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix + + HGMat = matmul(transpose(HGcomp),HGcomp) & + * HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix !-------------------------------------------------------------------------------------------------- ! init fields - restart: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file' - - write(rankStr,'(a1,i0)')'_',worldrank - fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') - - call HDF5_read(fileHandle,F_aim, 'F_aim') - call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') - call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') - call HDF5_read(fileHandle,F, 'F') - call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') - call HDF5_read(fileHandle,u_current, 'u') - call HDF5_read(fileHandle,u_lastInc, 'u_lastInc') - - elseif (restartInc == 0) then restart - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) - endif restart - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call utilities_updateIPcoords(F) - call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 - F, & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition - call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) - CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) - CHKERRQ(ierr) - - restartRead: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' - call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') - call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') - call HDF5_closeFile(fileHandle) - endif restartRead + restart: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file' + + write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_read(fileHandle,u_current, 'u') + call HDF5_read(fileHandle,u_lastInc, 'u_lastInc') + + elseif (restartInc == 0) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) + endif restart + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + call utilities_updateIPcoords(F) + call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + F, & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition + call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) + CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) + CHKERRQ(ierr) + + restartRead: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) + endif restartRead end subroutine grid_mech_FEM_init @@ -244,58 +243,57 @@ end subroutine grid_mech_FEM_init !> @brief solution for the FEM scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) - use IO, only: & - IO_error - use spectral_utilities, only: & - tBoundaryCondition, & - utilities_maskedCompliance - use FEsolving, only: & - restartWrite, & - terminallyIll + use IO, only: & + IO_error + use spectral_utilities, only: & + tBoundaryCondition, & + utilities_maskedCompliance + use FEsolving, only: & + restartWrite, & + terminallyIll - implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution - character(len=*), intent(in) :: & - incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC - type(tSolutionState) :: & - solution + character(len=*), intent(in) :: & + incInfoIn + real(pReal), intent(in) :: & + timeinc, & !< time increment of current solution + timeinc_old !< time increment of last successful increment + type(tBoundaryCondition), intent(in) :: & + stress_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(tSolutionState) :: & + solution !-------------------------------------------------------------------------------------------------- ! PETSc Data - PetscErrorCode :: ierr - SNESConvergedReason :: reason - - incInfo = incInfoIn + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old + params%stress_mask = stress_BC%maskFloat + params%stress_BC = stress_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr);CHKERRQ(ierr) + call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr);CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr) - - solution%converged = reason > 0 - solution%iterationsNeeded = totalIter - solution%termIll = terminallyIll - terminallyIll = .false. + call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr) + + solution%converged = reason > 0 + solution%iterationsNeeded = totalIter + solution%termIll = terminallyIll + terminallyIll = .false. end function grid_mech_FEM_solution @@ -326,13 +324,12 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat use FEsolving, only: & restartWrite - implicit none logical, intent(in) :: & guess real(pReal), intent(in) :: & timeinc_old, & timeinc, & - loadCaseTime !< remaining time of current load case + loadCaseTime !< remaining time of current load case type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC @@ -344,15 +341,15 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc - call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) if (cutBack) then C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? else !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? + if (restartWrite) then write(6,'(/,a)') ' writing converged results for restart';flush(6) write(rankStr,'(a1,i0)')'_',worldrank @@ -372,7 +369,7 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat call HDF5_closeFile(fileHandle) endif - call CPFEM_age ! age state and kinematics + call CPFEM_age ! age state and kinematics call utilities_updateIPcoords(F) C_volAvgLastInc = C_volAvg @@ -382,13 +379,13 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F F_aimDot = & F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed F_aimDot = & F_aimDot + deformation_BC%maskFloat * deformation_BC%values - elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed F_aimDot = & F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime endif @@ -403,20 +400,20 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat endif call VecCopy(solution_current,solution_lastInc,ierr); CHKERRQ(ierr) - F_lastInc = F ! winding F forward - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + F_lastInc = F ! winding F forward + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent endif !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients - F_aim = F_aim_lastInc + F_aimDot * timeinc - call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) + F_aim = F_aim_lastInc + F_aimDot * timeinc + call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) - CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) - CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) + CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) + CHKERRQ(ierr) end subroutine grid_mech_FEM_forward @@ -424,58 +421,57 @@ end subroutine grid_mech_FEM_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr) - use mesh - use spectral_utilities - use numerics, only: & - itmax, & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_stress_tolRel, & - err_stress_tolAbs - use FEsolving, only: & - terminallyIll +subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,ierr) + use mesh + use spectral_utilities + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_stress_tolRel, & + err_stress_tolAbs + use FEsolving, only: & + terminallyIll - implicit none - SNES :: snes_local - PetscInt, intent(in) :: PETScIter - PetscReal, intent(in) :: & - devNull1, & - devNull2, & - devNull3 - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal) :: & - err_div, & - divTol, & - BCTol + SNES :: snes_local + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + fnorm + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal) :: & + err_div, & + divTol, & + BCTol - err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ - divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) - BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) + err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ + divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) - if ((totalIter >= itmin .and. & - all([ err_div/divTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then - reason = 1 - elseif (totalIter >= itmax) then - reason = -1 - else - reason = 0 - endif + if ((totalIter >= itmin .and. & + all([ err_div/divTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then + reason = -1 + else + reason = 0 + endif !-------------------------------------------------------------------------------------------------- ! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' - write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine converged @@ -485,135 +481,134 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- subroutine formResidual(da_local,x_local, & f_local,dummy,ierr) - use numerics, only: & - itmax, & - itmin - use numerics, only: & - worldrank - use mesh, only: & - grid - use math, only: & - math_rotate_backward33, & - math_mul3333xx33 - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation - use spectral_utilities, only: & - utilities_constitutiveResponse - use IO, only: & - IO_intOut - use FEsolving, only: & - terminallyIll - use homogenization, only: & - materialpoint_dPdF + use numerics, only: & + itmax, & + itmin + use numerics, only: & + worldrank + use mesh, only: & + grid + use math, only: & + math_rotate_backward33, & + math_mul3333xx33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use spectral_utilities, only: & + utilities_constitutiveResponse + use IO, only: & + IO_intOut + use FEsolving, only: & + terminallyIll + use homogenization, only: & + materialpoint_dPdF - implicit none - DM :: da_local - Vec :: x_local, f_local - PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal - PetscScalar, dimension(8,3) :: x_elem, f_elem - PetscInt :: i, ii, j, jj, k, kk, ctr, ele - real(pReal), dimension(3,3) :: & - deltaF_aim - PetscInt :: & - PETScIter, & - nfuncs - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal), dimension(3,3,3,3) :: devNull + DM :: da_local + Vec :: x_local, f_local + PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal + PetscScalar, dimension(8,3) :: x_elem, f_elem + PetscInt :: i, ii, j, jj, k, kk, ctr, ele + real(pReal), dimension(3,3) :: & + deltaF_aim + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal), dimension(3,3,3,3) :: devNull - call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) + call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) - if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment !-------------------------------------------------------------------------------------------------- ! begin of new iteration - newIteration: if (totalIter <= PETScIter) then - totalIter = totalIter + 1 - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & - trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax - if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) - flush(6) - endif newIteration + newIteration: if (totalIter <= PETScIter) then + totalIter = totalIter + 1 + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim =', transpose(F_aim) + flush(6) + endif newIteration !-------------------------------------------------------------------------------------------------- ! get deformation gradient - call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) - do k = zstart, zend; do j = ystart, yend; do i = xstart, xend - ctr = 0 - do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 - ctr = ctr + 1 - x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) - enddo; enddo; enddo - ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 - F(1:3,1:3,ii,jj,kk) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem)) - enddo; enddo; enddo - call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) + enddo; enddo; enddo + ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 + F(1:3,1:3,ii,jj,kk) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem)) + enddo; enddo; enddo + call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(P_current,& - P_av,C_volAvg,devNull, & - F,params%timeinc,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + call Utilities_constitutiveResponse(P_current,& + P_av,C_volAvg,devNull, & + F,params%timeinc,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim_lastIter = F_aim - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) - F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + F_aim_lastIter = F_aim + deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) + F_aim = F_aim - deltaF_aim + err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- ! constructing residual - call VecSet(f_local,0.0_pReal,ierr);CHKERRQ(ierr) - call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) - call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) - ele = 0 - do k = zstart, zend; do j = ystart, yend; do i = xstart, xend - ctr = 0 - do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 - ctr = ctr + 1 - x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) - enddo; enddo; enddo - ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 - ele = ele + 1 - f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & - matmul(HGMat,x_elem)*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - ctr = 0 - do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 - ctr = ctr + 1 - f_scal(0:2,i+ii,j+jj,k+kk) = f_scal(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3) - enddo; enddo; enddo - enddo; enddo; enddo - call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + call VecSet(f_local,0.0_pReal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + ele = 0 + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) + enddo; enddo; enddo + ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 + ele = ele + 1 + f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & + matmul(HGMat,x_elem)*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + f_scal(0:2,i+ii,j+jj,k+kk) = f_scal(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3) + enddo; enddo; enddo + enddo; enddo; enddo + call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! applying boundary conditions - call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) - if (zstart == 0) then - f_scal(0:2,xstart,ystart,zstart) = 0.0 - f_scal(0:2,xend+1,ystart,zstart) = 0.0 - f_scal(0:2,xstart,yend+1,zstart) = 0.0 - f_scal(0:2,xend+1,yend+1,zstart) = 0.0 - endif - if (zend + 1 == grid(3)) then - f_scal(0:2,xstart,ystart,zend+1) = 0.0 - f_scal(0:2,xend+1,ystart,zend+1) = 0.0 - f_scal(0:2,xstart,yend+1,zend+1) = 0.0 - f_scal(0:2,xend+1,yend+1,zend+1) = 0.0 - endif - call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + if (zstart == 0) then + f_scal(0:2,xstart,ystart,zstart) = 0.0 + f_scal(0:2,xend+1,ystart,zstart) = 0.0 + f_scal(0:2,xstart,yend+1,zstart) = 0.0 + f_scal(0:2,xend+1,yend+1,zstart) = 0.0 + endif + if (zend + 1 == grid(3)) then + f_scal(0:2,xstart,ystart,zend+1) = 0.0 + f_scal(0:2,xend+1,ystart,zend+1) = 0.0 + f_scal(0:2,xstart,yend+1,zend+1) = 0.0 + f_scal(0:2,xend+1,yend+1,zend+1) = 0.0 + endif + call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) end subroutine formResidual @@ -622,97 +617,96 @@ end subroutine formResidual !> @brief forms the FEM stiffness matrix !-------------------------------------------------------------------------------------------------- subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr) - use mesh, only: & - mesh_ipCoordinates - use homogenization, only: & - materialpoint_dPdF - - implicit none - - DM :: da_local - Vec :: x_local, coordinates - Mat :: Jac_pre, Jac - MatStencil,dimension(4,24) :: row, col - PetscScalar,pointer,dimension(:,:,:,:) :: x_scal - PetscScalar,dimension(24,24) :: K_ele - PetscScalar,dimension(9,24) :: BMatFull - PetscInt :: i, ii, j, jj, k, kk, ctr, ele - PetscInt,dimension(3) :: rows - PetscScalar :: diag - PetscObject :: dummy - MatNullSpace :: matnull - PetscErrorCode :: ierr + use mesh, only: & + mesh_ipCoordinates + use homogenization, only: & + materialpoint_dPdF - BMatFull = 0.0 - BMatFull(1:3,1 :8 ) = BMat - BMatFull(4:6,9 :16) = BMat - BMatFull(7:9,17:24) = BMat - call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr) - call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr) - call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) - ele = 0 - do k = zstart, zend; do j = ystart, yend; do i = xstart, xend - ctr = 0 - do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 - ctr = ctr + 1 - col(MatStencil_i,ctr ) = i+ii - col(MatStencil_j,ctr ) = j+jj - col(MatStencil_k,ctr ) = k+kk - col(MatStencil_c,ctr ) = 0 - col(MatStencil_i,ctr+8 ) = i+ii - col(MatStencil_j,ctr+8 ) = j+jj - col(MatStencil_k,ctr+8 ) = k+kk - col(MatStencil_c,ctr+8 ) = 1 - col(MatStencil_i,ctr+16) = i+ii - col(MatStencil_j,ctr+16) = j+jj - col(MatStencil_k,ctr+16) = k+kk - col(MatStencil_c,ctr+16) = 2 - enddo; enddo; enddo - row = col - ele = ele + 1 - K_ele = 0.0 - K_ele(1 :8 ,1 :8 ) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele(9 :16,9 :16) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele(17:24,17:24) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele = K_ele + & - matmul(transpose(BMatFull), & - matmul(reshape(reshape(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,ele), & - shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ - call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) - CHKERRQ(ierr) - enddo; enddo; enddo - call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) - call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) - call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) - call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + + DM :: da_local + Vec :: x_local, coordinates + Mat :: Jac_pre, Jac + MatStencil,dimension(4,24) :: row, col + PetscScalar,pointer,dimension(:,:,:,:) :: x_scal + PetscScalar,dimension(24,24) :: K_ele + PetscScalar,dimension(9,24) :: BMatFull + PetscInt :: i, ii, j, jj, k, kk, ctr, ele + PetscInt,dimension(3) :: rows + PetscScalar :: diag + PetscObject :: dummy + MatNullSpace :: matnull + PetscErrorCode :: ierr + + BMatFull = 0.0 + BMatFull(1:3,1 :8 ) = BMat + BMatFull(4:6,9 :16) = BMat + BMatFull(7:9,17:24) = BMat + call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr) + call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr) + call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) + ele = 0 + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + col(MatStencil_i,ctr ) = i+ii + col(MatStencil_j,ctr ) = j+jj + col(MatStencil_k,ctr ) = k+kk + col(MatStencil_c,ctr ) = 0 + col(MatStencil_i,ctr+8 ) = i+ii + col(MatStencil_j,ctr+8 ) = j+jj + col(MatStencil_k,ctr+8 ) = k+kk + col(MatStencil_c,ctr+8 ) = 1 + col(MatStencil_i,ctr+16) = i+ii + col(MatStencil_j,ctr+16) = j+jj + col(MatStencil_k,ctr+16) = k+kk + col(MatStencil_c,ctr+16) = 2 + enddo; enddo; enddo + row = col + ele = ele + 1 + K_ele = 0.0 + K_ele(1 :8 ,1 :8 ) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(9 :16,9 :16) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(17:24,17:24) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele = K_ele + & + matmul(transpose(BMatFull), & + matmul(reshape(reshape(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,ele), & + shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ + call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) + CHKERRQ(ierr) + enddo; enddo; enddo + call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! applying boundary conditions - rows = [0, 1, 2] - diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + & - C_volAvg(2,2,2,2)/delta(2)**2.0_pReal + & - C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ - call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr) - CHKERRQ(ierr) - call DMGetGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) - call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) - ele = 0 - do k = zstart, zend; do j = ystart, yend; do i = xstart, xend - ele = ele + 1 - x_scal(0:2,i,j,k) = mesh_ipCoordinates(1:3,1,ele) - enddo; enddo; enddo - call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates) - call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes - call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) - call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) - call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) - call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr) + rows = [0, 1, 2] + diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + & + C_volAvg(2,2,2,2)/delta(2)**2.0_pReal + & + C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ + call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr) + CHKERRQ(ierr) + call DMGetGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) + ele = 0 + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ele = ele + 1 + x_scal(0:2,i,j,k) = mesh_ipCoordinates(1:3,1,ele) + enddo; enddo; enddo + call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates) + call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes + call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) + call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) + call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) + call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr) end subroutine formJacobian diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90 index 91dcd03b4..f698fdc6d 100644 --- a/src/grid_mech_spectral_basic.f90 +++ b/src/grid_mech_spectral_basic.f90 @@ -107,7 +107,6 @@ subroutine grid_mech_spectral_basic_init use math, only: & math_invSym3333 - implicit none real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal @@ -221,7 +220,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ restartWrite, & terminallyIll - implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & @@ -298,7 +296,6 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi use FEsolving, only: & restartWrite - implicit none logical, intent(in) :: & guess real(pReal), intent(in) :: & @@ -397,7 +394,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm use FEsolving, only: & terminallyIll - implicit none SNES :: snes_local PetscInt, intent(in) :: PETScIter PetscReal, intent(in) :: & @@ -468,7 +464,6 @@ subroutine formResidual(in, F, & use FEsolving, only: & terminallyIll - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), & intent(in) :: F !< deformation gradient field diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index a7bf75209..55df4455f 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -113,7 +113,6 @@ subroutine grid_mech_spectral_polarisation_init use math, only: & math_invSym3333 - implicit none real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal @@ -238,7 +237,6 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, restartWrite, & terminallyIll - implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & @@ -445,7 +443,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm use FEsolving, only: & terminallyIll - implicit none SNES :: snes_local PetscInt, intent(in) :: PETScIter PetscReal, intent(in) :: & @@ -529,7 +526,6 @@ subroutine formResidual(in, FandF_tau, & use FEsolving, only: & terminallyIll - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), & target, intent(in) :: FandF_tau diff --git a/src/grid_thermal_spectral.f90 b/src/grid_thermal_spectral.f90 index 31740feb9..e899fd89a 100644 --- a/src/grid_thermal_spectral.f90 +++ b/src/grid_thermal_spectral.f90 @@ -69,7 +69,6 @@ subroutine grid_thermal_spectral_init worldsize, & petsc_options - implicit none PetscInt, dimension(worldsize) :: localK integer :: i, j, k, cell DM :: thermal_grid @@ -167,7 +166,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result use thermal_conduction, only: & thermal_conduction_putTemperatureAndItsRate - implicit none real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment @@ -242,7 +240,6 @@ subroutine grid_thermal_spectral_forward thermal_conduction_getMassDensity, & thermal_conduction_getSpecificHeat - implicit none integer :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal @@ -311,7 +308,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) thermal_conduction_getMassDensity, & thermal_conduction_getSpecificHeat - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & From ed8270d432ab971b09d8f4838fcf700470492718 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Apr 2019 14:53:36 +0200 Subject: [PATCH 12/12] need to read reference stiffness --- src/grid_mech_spectral_basic.f90 | 3 +++ src/grid_mech_spectral_polarisation.f90 | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90 index f698fdc6d..2daebefbd 100644 --- a/src/grid_mech_spectral_basic.f90 +++ b/src/grid_mech_spectral_basic.f90 @@ -116,6 +116,7 @@ subroutine grid_mech_spectral_basic_init F ! pointer to solution data PetscInt, dimension(worldsize) :: localK integer(HID_T) :: fileHandle + integer :: fileUnit character(len=1024) :: rankStr write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>' @@ -201,6 +202,8 @@ subroutine grid_mech_spectral_basic_init call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_closeFile(fileHandle) + fileUnit = IO_open_jobFile_binary('C_ref') + read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead call utilities_updateGamma(C_minMaxAvg,.true.) diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index 55df4455f..3bd30a360 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -124,6 +124,7 @@ subroutine grid_mech_spectral_polarisation_init F_tau ! specific (sub)pointer PetscInt, dimension(worldsize) :: localK integer(HID_T) :: fileHandle + integer :: fileUnit character(len=1024) :: rankStr write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>' @@ -214,6 +215,8 @@ subroutine grid_mech_spectral_polarisation_init call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_closeFile(fileHandle) + fileUnit = IO_open_jobFile_binary('C_ref') + read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead call utilities_updateGamma(C_minMaxAvg,.true.)