diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index 8909f10e3..85c856076 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -355,18 +355,14 @@ program DAMASK_spectral_Driver if (appendToOutFile) then ! after restart, append to existing results file open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& '.spectralOut',form='UNFORMATTED', position='APPEND', status='OLD') - ! '.spectralOut',access='STREAM', position='APPEND', status='OLD') open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& '.sta',form='FORMATTED', position='APPEND', status='OLD') else ! open new files ... open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& '.spectralOut',form='UNFORMATTED',status='REPLACE') - ! '.spectralOut',access='STREAM',status='REPLACE') write(resUnit) 'load', trim(loadCaseFile) ! ... and write header write(resUnit) 'workingdir', trim(getSolverWorkingDirectoryName()) write(resUnit) 'geometry', trim(geometryFile) - !write(resUnit) 'grid', grid - !write(resUnit) 'size', geomSize write(resUnit) 'resolution', grid write(resUnit) 'dimension', geomSize write(resUnit) 'materialpoint_sizeResults', materialpoint_sizeResults diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index d56168f6f..0b8e69505 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -163,9 +163,6 @@ subroutine AL_init(temperature) real(pReal), dimension(:,:,:,:,:), allocatable :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal - real(pReal), dimension(3,3,3,3) :: & - temp3333_Real = 0.0_pReal, & - temp3333_Real2 = 0.0_pReal PetscErrorCode :: ierr PetscObject :: dummy @@ -214,7 +211,8 @@ subroutine AL_init(temperature) F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) F_lambda = F F_lambda_lastInc = F_lastInc - elseif (restartInc > 1_pInt) then + + elseif (restartInc > 1_pInt) then ! read in F to calculate coordinates and initialize CPFEM general if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading values of increment ', restartInc - 1_pInt, ' from file' @@ -227,7 +225,7 @@ subroutine AL_init(temperature) close (777) call IO_read_realFile(777,'F_lastInc2', trim(getSolverJobName()),size(F_lastInc2)) read (777,rec=1) F_lastInc2 - close (777) + close (777) call IO_read_realFile(777,'F_lambda',trim(getSolverJobName()),size(F_lambda)) read (777,rec=1) F_lambda close (777) @@ -244,28 +242,35 @@ subroutine AL_init(temperature) call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) read (777,rec=1) f_aimDot close (777) + endif + + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) + call Utilities_constitutiveResponse(F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]),& + temperature,0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) + nullify(F) + nullify(F_lambda) + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc + + if (restartInc > 1_pInt) then ! using old values from files + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading more values of increment', restartInc - 1_pInt, 'from file' + flush(6) + call IO_read_realFile(777,'F_lambda',trim(getSolverJobName()),size(F_lambda)) + read (777,rec=1) F_lambda call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) read (777,rec=1) C_volAvg close (777) call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) read (777,rec=1) C_volAvgLastInc close (777) - call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real)) + call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) read (777,rec=1) C_minMaxAvg close (777) endif - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& - F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) - call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real,temp3333_Real2,& - temp33_Real,.false.,math_I3) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) -!-------------------------------------------------------------------------------------------------- -! reference stiffness - if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness - C_minMaxAvg = temp3333_Real2 - C_volAvg = temp3333_Real - endif + call Utilities_updateGamma(C_minMaxAvg,.True.) C_scale = C_minMaxAvg @@ -419,11 +424,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) integer(pInt) :: & i, j, k, e - F => x_scal(1:3,1:3,1,& + F => x_scal(1:3,1:3,1,& XG_RANGE,YG_RANGE,ZG_RANGE) - F_lambda => x_scal(1:3,1:3,2,& + F_lambda => x_scal(1:3,1:3,2,& XG_RANGE,YG_RANGE,ZG_RANGE) - residual_F => f_scal(1:3,1:3,1,& + residual_F => f_scal(1:3,1:3,1,& X_RANGE,Y_RANGE,Z_RANGE) residual_F_lambda => f_scal(1:3,1:3,2,& X_RANGE,Y_RANGE,Z_RANGE) @@ -434,7 +439,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - if (totalIter <= PETScIter) then ! new iteration + if(totalIter <= PETScIter) then ! new iteration !-------------------------------------------------------------------------------------------------- ! report begin of new iteration totalIter = totalIter + 1_pInt @@ -507,7 +512,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) e = 0_pInt do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) e = e + 1_pInt - residual_F(1:3,1:3,i,j,k) = math_mul3333xx33(math_invSym3333(materialpoint_dPdF(:,:,:,:,1,e) + C_scale), & + 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) - & math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3))) & @@ -641,7 +646,7 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ ! update coordinates and rate and forward last inc call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) F => xx_psc(0:8,:,:,:) - F_lambda => xx_psc(9:17,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) if (restartWrite) then write(6,'(/,a)') ' writing converged results for restart' flush(6) @@ -675,6 +680,7 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) write (777,rec=1) C_volAvgLastInc close(777) + endif mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) @@ -682,13 +688,13 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ if (cutBack) then F_aim = F_aim_lastInc F_lambda = reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid(3)]) - F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)]) + F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)]) C_volAvg = C_volAvgLastInc else ForwardData = .True. + C_volAvgLastInc = C_volAvg !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - C_volAvgLastInc = C_volAvg if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed @@ -707,8 +713,8 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)])) - F_lastInc2 = F_lastInc - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) + F_lastInc2 = F_lastInc + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)]) endif @@ -720,13 +726,13 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition if (.not. guess) then ! large strain forwarding do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) - F_lambda33 = reshape(F_lambda(:,i,j,k),[3,3]) + F_lambda33 = reshape(F_lambda(1:9,i,j,k),[3,3]) F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & math_mul3333xx33(C_scale,& math_mul33x33(math_transpose33(F_lambda33),& F_lambda33) -math_I3))*0.5_pReal)& + math_I3 - F_lambda(:,i,j,k) = reshape(F_lambda33,[9]) + F_lambda(1:9,i,j,k) = reshape(F_lambda33,[9]) enddo; enddo; enddo endif call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index d070df3d9..3d6ca14b6 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -119,7 +119,7 @@ subroutine basic_init(temperature) elseif (restartInc > 1_pInt) then ! using old values from file if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading deformation gradients of increment', restartInc - 1_pInt, 'from file' + 'reading values of increment ', restartInc - 1_pInt, ' from file' flush(6) call IO_read_realFile(777,'F',trim(getSolverJobName()),size(F)) read (777,rec=1) F @@ -127,31 +127,30 @@ subroutine basic_init(temperature) call IO_read_realFile(777,'F_lastInc',trim(getSolverJobName()),size(F_lastInc)) read (777,rec=1) F_lastInc close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) + read (777,rec=1) f_aimDot + close (777) + F_aim = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt ! average of F + F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc endif mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,F),[3,1,product(grid)]) - call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,C_minmaxAvg,& + call Utilities_constitutiveResponse(F_lastInc,F,temperature,0.0_pReal,P,C,C_minMaxAvg,& temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness - if (restartInc > 1_pInt) then ! using old values from file - F_aim = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt ! average of F - F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc - + if (restartInc > 1_pInt) then ! using old values from files if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading more values of increment', restartInc - 1_pInt, 'from file' flush(6) - call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) - read (777,rec=1) f_aimDot - close (777) call IO_read_realFile(777,'C',trim(getSolverJobName()),size(C)) read (777,rec=1) C close (777) call IO_read_realFile(777,'C_lastInc',trim(getSolverJobName()),size(C_lastInc)) read (777,rec=1) C_lastInc close (777) - call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minmaxAvg)) - read (777,rec=1) C_minmaxAvg + call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) + read (777,rec=1) C_minMaxAvg close (777) endif diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index a1a0b66dd..ba727785f 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -198,7 +198,7 @@ subroutine basicPETSc_init(temperature) elseif (restartInc > 1_pInt) then ! using old values from file if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading deformation gradients of increment', restartInc - 1_pInt, 'from file' + 'reading values of increment ', restartInc - 1_pInt, ' from file' flush(6) call IO_read_realFile(777,'F',trim(getSolverJobName()),size(F)) read (777,rec=1) F @@ -209,35 +209,33 @@ subroutine basicPETSc_init(temperature) call IO_read_realFile(777,'F_lastInc2',trim(getSolverJobName()),size(F_lastInc2)) read (777,rec=1) F_lastInc2 close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) + read (777,rec=1) f_aimDot + close (777) + F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F + F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc endif mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) - call Utilities_constitutiveResponse(& + call Utilities_constitutiveResponse(F_lastInc, & reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),& - reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),& - temperature,0.0_pReal,P,C_volAvg,C_minmaxAvg,temp33_Real,.false.,math_I3) - call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back into PETSc - - if (restartInc > 1_pInt) then ! using old values from file - F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F - F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc + temperature,0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) + call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc + if (restartInc > 1_pInt) then ! using old values from files if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading more values of increment', restartInc - 1_pInt, 'from file' flush(6) - call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) - read (777,rec=1) f_aimDot - close (777) call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) read (777,rec=1) C_volAvg close (777) call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) read (777,rec=1) C_volAvgLastInc close (777) - call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minmaxAvg)) - read (777,rec=1) C_minmaxAvg + call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) + read (777,rec=1) C_minMaxAvg close (777) endif diff --git a/code/DAMASK_spectral_solverPolarisation.f90 b/code/DAMASK_spectral_solverPolarisation.f90 index 6782f2120..cec0ea945 100644 --- a/code/DAMASK_spectral_solverPolarisation.f90 +++ b/code/DAMASK_spectral_solverPolarisation.f90 @@ -163,9 +163,6 @@ subroutine Polarisation_init(temperature) real(pReal), dimension(:,:,:,:,:), allocatable :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal - real(pReal), dimension(3,3,3,3) :: & - temp3333_Real = 0.0_pReal, & - temp3333_Real2 = 0.0_pReal PetscErrorCode :: ierr PetscObject :: dummy @@ -211,9 +208,9 @@ subroutine Polarisation_init(temperature) if (restartInc == 1_pInt) then ! no deformation (no restart) F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity F_lastInc2 = F_lastInc - F_tau_lastInc = 2.0_pReal*F_lastInc F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) F_tau = 2.0_pReal* F + F_tau_lastInc = 2.0_pReal*F_lastInc elseif (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & @@ -228,8 +225,6 @@ subroutine Polarisation_init(temperature) call IO_read_realFile(777,'F_lastInc2',trim(getSolverJobName()),size(F_lastInc2)) read (777,rec=1) F_lastInc2 close (777) - F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F - F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc call IO_read_realFile(777,'F_tau',trim(getSolverJobName()),size(F_tau)) read (777,rec=1) F_tau close (777) @@ -246,32 +241,35 @@ subroutine Polarisation_init(temperature) call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) read (777,rec=1) f_aimDot close (777) + endif + + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) + call Utilities_constitutiveResponse(F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]),& + temperature,0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) + nullify(F) + nullify(F_tau) + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc + + if (restartInc > 1_pInt) then ! using old values from files + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading more values of increment', restartInc - 1_pInt, 'from file' + flush(6) call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) read (777,rec=1) C_volAvg close (777) call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) read (777,rec=1) C_volAvgLastInc close (777) - call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real)) + call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) read (777,rec=1) C_minMaxAvg close (777) endif - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& - F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) - call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real,temp3333_Real2,& - temp33_Real,.false.,math_I3) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) -!-------------------------------------------------------------------------------------------------- -! reference stiffness - if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness - C_minMaxAvg = temp3333_Real2 - C_volAvg = temp3333_Real - endif - - call Utilities_updateGamma(temp3333_Real2,.True.) - C_scale = temp3333_Real2 - S_scale = math_invSym3333(temp3333_Real2) + call Utilities_updateGamma(C_minMaxAvg,.True.) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) end subroutine Polarisation_init @@ -508,7 +506,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) e = e + 1_pInt residual_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(materialpoint_dPdF(:,:,:,:,1,e) + C_scale), & + 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) - math_mul33x33(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) @@ -720,13 +718,13 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition if (.not. guess) then ! large strain forwarding do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) - F_lambda33 = reshape(F_tau(:,i,j,k)-F(:,i,j,k),[3,3]) + F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(:,i,j,k),[3,3]) F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & math_mul3333xx33(C_scale,& math_mul33x33(math_transpose33(F_lambda33),& F_lambda33) -math_I3))*0.5_pReal)& + math_I3 - F_tau(:,i,j,k) = reshape(F_lambda33,[9])+F(:,i,j,k) + F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(:,i,j,k) enddo; enddo; enddo endif call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 5e3e047d5..724bb4106 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -862,6 +862,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& if (forwardData) then ! aging results calcMode = ior(calcMode, CPFEM_AGERESULTS) collectMode = ior(collectMode, CPFEM_BACKUPJACOBIAN) + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)]) endif if (cutBack) then ! restore saved variables collectMode = ior(collectMode , CPFEM_RESTOREJACOBIAN) @@ -872,8 +873,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& call CPFEM_general(collectMode,usePingPong,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), & ! collect mode handles Jacobian backup / restoration crystallite_temperature(1,1),timeinc,1_pInt,1_pInt) - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)]) - materialpoint_F = reshape(F, [3,3,1,product(grid)]) + materialpoint_F = reshape(F,[3,3,1,product(grid)]) call debug_reset()