diff --git a/src/IO.f90 b/src/IO.f90 index 29107b59c..cd7c09c75 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -503,8 +503,6 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) ! errors related to the grid solver case (831) msg = 'mask consistency violated in grid load case' - case (832) - msg = 'ill-defined L (line partly defined) in grid load case' case (833) msg = 'non-positive ratio for geometric progression' case (834) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index b6e45a75b..7ffa6b80a 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -53,12 +53,11 @@ program DAMASK_grid integer, parameter :: & subStepFactor = 2 !< for each substep, divide the last time increment by 2.0 real(pReal) :: & - T_0 = 300.0_pReal, & - time = 0.0_pReal, & !< elapsed time - time0 = 0.0_pReal, & !< begin of interval - timeinc = 1.0_pReal, & !< current time interval - timeIncOld = 0.0_pReal, & !< previous time interval - remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case + t = 0.0_pReal, & !< elapsed time + t_0 = 0.0_pReal, & !< begin of interval + Delta_t = 1.0_pReal, & !< current time interval + Delta_t_prev = 0.0_pReal, & !< previous time interval + t_remaining = 0.0_pReal !< remaining time of current load case logical :: & guess, & !< guess along former trajectory stagIterate, & @@ -230,12 +229,6 @@ program DAMASK_grid reportAndCheck: if (worldrank == 0) then print'(/,a,i0)', ' load case: ', l print*, ' estimate_rate:', loadCases(l)%estimate_rate - if (loadCases(l)%deformation%myType == 'L') then - do j = 1, 3 - if (any(loadCases(l)%deformation%mask(j,1:3) .eqv. .true.) .and. & - any(loadCases(l)%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined - enddo - endif if (loadCases(l)%deformation%myType == 'F') then print*, ' F:' else @@ -243,14 +236,14 @@ program DAMASK_grid endif do i = 1, 3; do j = 1, 3 if (loadCases(l)%deformation%mask(i,j)) then - write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j) - else write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' + else + write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j) endif enddo; write(IO_STDOUT,'(/)',advance='no') enddo if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831 - if (any(loadCases(l)%stress%mask .and. transpose(loadCases(l)%stress%mask) .and. (math_I3<1))) & + if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) & errorID = 838 ! no rotation is allowed by stress BC if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:' @@ -259,9 +252,9 @@ program DAMASK_grid if (loadCases(l)%stress%myType /= '') then do i = 1, 3; do j = 1, 3 if (loadCases(l)%stress%mask(i,j)) then - write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal - else write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' + else + write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal endif enddo; write(IO_STDOUT,'(/)',advance='no') enddo @@ -303,7 +296,7 @@ program DAMASK_grid case(FIELD_THERMAL_ID) initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict) thermal => initial_conditions%get('thermal',defaultVal=emptyDict) - call grid_thermal_spectral_init(thermal%get_asFloat('T',defaultVal = T_0)) + call grid_thermal_spectral_init(thermal%get_asFloat('T')) case(FIELD_DAMAGE_ID) call grid_damage_spectral_init @@ -330,7 +323,7 @@ program DAMASK_grid endif writeUndeformed loadCaseLooping: do l = 1, size(loadCases) - time0 = time ! load case start time + t_0 = t ! load case start time guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc incLooping: do inc = 1, loadCases(l)%N @@ -338,31 +331,31 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! forwarding time - timeIncOld = timeinc ! last timeinc that brought former inc to an end + Delta_t_prev = Delta_t ! last time intervall that brought former inc to an end if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then ! linear scale - timeinc = loadCases(l)%t/real(loadCases(l)%N,pReal) + Delta_t = loadCases(l)%t/real(loadCases(l)%N,pReal) else - timeinc = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) & + Delta_t = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) & / (1.0_pReal-loadCases(l)%r**loadCases(l)%N) endif - timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step + Delta_t = Delta_t * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc? - time = time + timeinc ! just advance time, skip already performed calculation + t = t + Delta_t ! just advance time, skip already performed calculation guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference else skipping stepFraction = 0 ! fraction scaled by stepFactor**cutLevel subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) - remainingLoadCaseTime = loadCases(l)%t+time0 - time - time = time + timeinc ! forward target time + t_remaining = loadCases(l)%t + t_0 - t + t = t + Delta_t ! forward target time stepFraction = stepFraction + 1 ! count step !-------------------------------------------------------------------------------------------------- ! report begin of new step print'(/,a)', ' ###########################################################################' print'(1x,a,es12.5,6(a,i0))', & - 'Time', time, & + 'Time', t, & 's: Increment ', inc,'/',loadCases(l)%N,& '-', stepFraction,'/',subStepFactor**cutBackLevel,& ' of load case ', l,'/',size(loadCases) @@ -377,10 +370,10 @@ program DAMASK_grid select case(ID(field)) case(FIELD_MECH_ID) call mechanical_forward (& - cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & - deformation_BC = loadCases(l)%deformation, & - stress_BC = loadCases(l)%stress, & - rotation_BC = loadCases(l)%rot) + cutBack,guess,Delta_t,Delta_t_prev,t_remaining, & + deformation_BC = loadCases(l)%deformation, & + stress_BC = loadCases(l)%stress, & + rotation_BC = loadCases(l)%rot) case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) @@ -398,9 +391,9 @@ program DAMASK_grid case(FIELD_MECH_ID) solres(field) = mechanical_solution(incInfo) case(FIELD_THERMAL_ID) - solres(field) = grid_thermal_spectral_solution(timeinc) + solres(field) = grid_thermal_spectral_solution(Delta_t) case(FIELD_DAMAGE_ID) - solres(field) = grid_damage_spectral_solution(timeinc) + solres(field) = grid_damage_spectral_solution(Delta_t) end select if (.not. solres(field)%converged) exit ! no solution found @@ -418,11 +411,11 @@ program DAMASK_grid if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged .and. .not. solres(1)%termIll) then ! and acceptable solution found call mechanical_updateCoords - timeIncOld = timeinc + Delta_t_prev = Delta_t cutBack = .false. guess = .true. ! start guessing after first converged (sub)inc if (worldrank == 0) then - write(statUnit,*) totalIncsCounter, time, cutBackLevel, & + write(statUnit,*) totalIncsCounter, t, cutBackLevel, & solres(1)%converged, solres(1)%iterationsNeeded flush(statUnit) endif @@ -430,8 +423,8 @@ program DAMASK_grid cutBack = .true. stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator cutBackLevel = cutBackLevel + 1 - time = time - timeinc ! rewind time - timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep + t = t - Delta_t + Delta_t = Delta_t/real(subStepFactor,pReal) ! cut timestep print'(/,a)', ' cutting back ' else ! no more options to continue if (worldrank == 0) close(statUnit) @@ -453,7 +446,7 @@ program DAMASK_grid if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then print'(1/,a)', ' ... writing results to file ......................................' flush(IO_STDOUT) - call CPFEM_results(totalIncsCounter,time) + call CPFEM_results(totalIncsCounter,t) endif if(signal) call interface_setSIGUSR1(.false.) call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) @@ -497,8 +490,8 @@ subroutine getMaskedTensor(values,mask,tensor) do i = 1,3 row => tensor%get(i) do j = 1,3 - mask(i,j) = row%get_asString(j) /= 'x' - if (mask(i,j)) values(i,j) = row%get_asFloat(j) + mask(i,j) = row%get_asString(j) == 'x' + if (.not. mask(i,j)) values(i,j) = row%get_asFloat(j) enddo enddo diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 162d665cb..f43a6cb18 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -160,10 +160,10 @@ end subroutine grid_damage_spectral_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral damage scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_damage_spectral_solution(timeinc) result(solution) +function grid_damage_spectral_solution(Delta_t) result(solution) real(pReal), intent(in) :: & - timeinc !< increment in time for current solution + Delta_t !< increment in time for current solution integer :: i, j, k, ce type(tSolutionState) :: solution PetscInt :: devNull @@ -176,7 +176,7 @@ function grid_damage_spectral_solution(timeinc) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - params%timeinc = timeinc + params%Delta_t = Delta_t call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) @@ -284,7 +284,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + mu_ref*phi_current(i,j,k) enddo; enddo; enddo @@ -292,7 +292,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! convolution of damage field with green operator call utilities_FFTscalarForward - call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc) + call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) call utilities_FFTscalarBackward where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) & diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 111e9fece..9cccd2fc0 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -324,7 +324,7 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution) solution%iterationsNeeded = totalIter solution%termIll = terminallyIll terminallyIll = .false. - P_aim = merge(P_aim,P_av,params%stress_mask) + P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_FEM_solution @@ -363,20 +363,20 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai else C_volAvgLastInc = C_volAvg - F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components + F_aimDot = merge(merge(.0_pReal,(F_aim-F_aim_lastInc)/Delta_t_old,stress_BC%mask),.0_pReal,guess) ! estimate deformation rate for prescribed stress components F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F F_aimDot = F_aimDot & - + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) + + matmul(merge(.0_pReal,deformation_BC%values,deformation_BC%mask),F_aim_lastInc) elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed F_aimDot = F_aimDot & - + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) + + merge(.0_pReal,deformation_BC%values,deformation_BC%mask) elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed F_aimDot = F_aimDot & - + merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask) + + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask) endif if (guess) then @@ -397,9 +397,9 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * Delta_t if (stress_BC%myType=='P') P_aim = P_aim & - + merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t if (stress_BC%myType=='dot_P') P_aim = P_aim & - + merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t call VecAXPY(solution_current,Delta_t,solution_rate,ierr); CHKERRQ(ierr) @@ -412,7 +412,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai ! set module wide available data params%stress_mask = stress_BC%mask params%rotation_BC = rotation_BC - params%timeinc = Delta_t + params%Delta_t = Delta_t end subroutine grid_mechanical_FEM_forward @@ -568,13 +568,13 @@ subroutine formResidual(da_local,x_local, & ! evaluate constitutive response call utilities_constitutiveResponse(P_current,& P_av,C_volAvg,devNull, & - F,params%timeinc,params%rotation_BC) + F,params%Delta_t,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc - err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) + err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index c4fce61a7..6f382d639 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -277,7 +277,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution) solution%iterationsNeeded = totalIter solution%termIll = terminallyIll terminallyIll = .false. - P_aim = merge(P_aim,P_av,params%stress_mask) + P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_spectral_basic_solution @@ -314,20 +314,20 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components + F_aimDot = merge(merge(.0_pReal,(F_aim-F_aim_lastInc)/Delta_t_old,stress_BC%mask),.0_pReal,guess) ! estimate deformation rate for prescribed stress components F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F F_aimDot = F_aimDot & - + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) + + matmul(merge(.0_pReal,deformation_BC%values,deformation_BC%mask),F_aim_lastInc) elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed F_aimDot = F_aimDot & - + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) + + merge(.0_pReal,deformation_BC%values,deformation_BC%mask) elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed F_aimDot = F_aimDot & - + merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask) + + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -342,9 +342,9 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * Delta_t if (stress_BC%myType=='P') P_aim = P_aim & - + merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t if (stress_BC%myType=='dot_P') P_aim = P_aim & - + merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) @@ -354,7 +354,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ ! set module wide available data params%stress_mask = stress_BC%mask params%rotation_BC = rotation_BC - params%timeinc = Delta_t + params%Delta_t = Delta_t end subroutine grid_mechanical_spectral_basic_forward @@ -492,14 +492,14 @@ subroutine formResidual(in, F, & ! evaluate constitutive response 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) + F,params%Delta_t,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) + err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index a52554e3c..7a4ae7595 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -309,7 +309,7 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti solution%iterationsNeeded = totalIter solution%termIll = terminallyIll terminallyIll = .false. - P_aim = merge(P_aim,P_av,params%stress_mask) + P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_spectral_polarisation_solution @@ -350,20 +350,20 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components + F_aimDot = merge(merge(.0_pReal,(F_aim-F_aim_lastInc)/Delta_t_old,stress_BC%mask),.0_pReal,guess) ! estimate deformation rate for prescribed stress components F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F F_aimDot = F_aimDot & - + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) + + matmul(merge(.0_pReal,deformation_BC%values,deformation_BC%mask),F_aim_lastInc) elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed F_aimDot = F_aimDot & - + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) + + merge(.0_pReal,deformation_BC%values,deformation_BC%mask) elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed F_aimDot = F_aimDot & - + merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask) + + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -382,9 +382,9 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * Delta_t if(stress_BC%myType=='P') P_aim = P_aim & - + merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t if(stress_BC%myType=='dot_P') P_aim = P_aim & - + merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),& @@ -408,7 +408,7 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D ! set module wide available data params%stress_mask = stress_BC%mask params%rotation_BC = rotation_BC - params%timeinc = Delta_t + params%Delta_t = Delta_t end subroutine grid_mechanical_spectral_polarisation_forward @@ -592,14 +592,14 @@ subroutine formResidual(in, FandF_tau, & ! 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/num%beta,params%timeinc,params%rotation_BC) + F - residual_F_tau/num%beta,params%Delta_t,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc - err_BC = maxval(abs(merge(P_av-P_aim, & - math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),& + err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), & + P_av-P_aim, & params%stress_mask))) ! calculate divergence tensorField_real = 0.0_pReal diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 47c49e76f..41fca2133 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -155,10 +155,10 @@ end subroutine grid_thermal_spectral_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_thermal_spectral_solution(timeinc) result(solution) +function grid_thermal_spectral_solution(Delta_t) result(solution) real(pReal), intent(in) :: & - timeinc !< increment in time for current solution + Delta_t !< increment in time for current solution integer :: i, j, k, ce type(tSolutionState) :: solution PetscInt :: devNull @@ -171,7 +171,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - params%timeinc = timeinc + params%Delta_t = Delta_t call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) @@ -195,7 +195,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce) + call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) enddo; enddo; enddo call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) @@ -233,7 +233,7 @@ subroutine grid_thermal_spectral_forward(cutBack) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce) + call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) enddo; enddo; enddo else T_lastInc = T_current @@ -279,7 +279,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & + scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + mu_ref*T_current(i,j,k) enddo; enddo; enddo @@ -287,7 +287,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! convolution of temperature field with green operator call utilities_FFTscalarForward - call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc) + call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) call utilities_FFTscalarBackward !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index bd37e2654..e33ee94d0 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -88,7 +88,7 @@ module spectral_utilities type, public :: tBoundaryCondition !< set of parameters defining a boundary condition real(pReal), dimension(3,3) :: values = 0.0_pReal - logical, dimension(3,3) :: mask = .false. + logical, dimension(3,3) :: mask = .true. character(len=:), allocatable :: myType end type tBoundaryCondition @@ -96,7 +96,7 @@ module spectral_utilities real(pReal), dimension(3,3) :: stress_BC logical, dimension(3,3) :: stress_mask type(rotation) :: rotation_BC - real(pReal) :: timeinc + real(pReal) :: Delta_t end type tSolutionParams type :: tNumerics @@ -684,7 +684,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) logical :: errmatinv character(len=pStringLen):: formatString - mask_stressVector = reshape(transpose(mask_stress), [9]) + mask_stressVector = .not. reshape(transpose(mask_stress), [9]) size_reduced = count(mask_stressVector) if(size_reduced > 0) then temp99_real = math_3333to99(rot_BC%rotate(C)) @@ -791,16 +791,16 @@ end subroutine utilities_fourierTensorDivergence !-------------------------------------------------------------------------------------------------- -!> @brief calculate constitutive response from homogenization_F0 to F during timeinc +!> @brief calculate constitutive response from homogenization_F0 to F during Delta_t !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& - F,timeinc,rotation_BC) + F,Delta_t,rotation_BC) real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target - real(pReal), intent(in) :: timeinc !< loading time + real(pReal), intent(in) :: Delta_t !< loading time type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame @@ -815,11 +815,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field - call homogenization_mechanical_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field + call homogenization_mechanical_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field if (.not. terminallyIll) & - call homogenization_thermal_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + call homogenization_thermal_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) if (.not. terminallyIll) & - call homogenization_mechanical_response2(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt @@ -870,14 +870,14 @@ end subroutine utilities_constitutiveResponse !-------------------------------------------------------------------------------------------------- -!> @brief calculates forward rate, either guessing or just add delta/timeinc +!> @brief calculates forward rate, either guessing or just add delta/Delta_t !-------------------------------------------------------------------------------------------------- pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) real(pReal), intent(in), dimension(3,3) :: & avRate !< homogeneous addon real(pReal), intent(in) :: & - dt !< timeinc between field0 and field + dt !< Delta_t between field0 and field logical, intent(in) :: & heterogeneous !< calculate field of rates real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & @@ -899,10 +899,10 @@ end function utilities_calculateRate !> @brief forwards a field with a pointwise given rate, if aim is given, !> ensures that the average matches the aim !-------------------------------------------------------------------------------------------------- -function utilities_forwardField(timeinc,field_lastInc,rate,aim) +function utilities_forwardField(Delta_t,field_lastInc,rate,aim) real(pReal), intent(in) :: & - timeinc !< timeinc of current step + Delta_t !< Delta_t of current step real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & field_lastInc, & !< initial field rate !< rate by which to forward @@ -913,7 +913,7 @@ function utilities_forwardField(timeinc,field_lastInc,rate,aim) real(pReal), dimension(3,3) :: fieldDiff !< - aim PetscErrorCode :: ierr - utilities_forwardField = field_lastInc + rate*timeinc + utilities_forwardField = field_lastInc + rate*Delta_t if (present(aim)) then !< correct to match average fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)