223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'more-flexible-L' into 'development' 

more flexibility for the L in the load case

See merge request damask/DAMASK!420
This commit is contained in:
Philip Eisenlohr 2021-08-09 21:27:13 +00:00
parent f1aa4e9d2c 574cfd7034
commit f75235f6a9
8 changed files with 90 additions and 99 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/IO.f90
View File

@ -503,8 +503,6 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
! errors related to the grid solver ! errors related to the grid solver
case (831) case (831)
msg = 'mask consistency violated in grid load case' msg = 'mask consistency violated in grid load case'
case (832)
msg = 'ill-defined L (line partly defined) in grid load case'
case (833) case (833)
msg = 'non-positive ratio for geometric progression' msg = 'non-positive ratio for geometric progression'
case (834) case (834)

67
src/grid/DAMASK_grid.f90
View File

@ -53,12 +53,11 @@ program DAMASK_grid
integer, parameter :: & integer, parameter :: &
subStepFactor = 2 !< for each substep, divide the last time increment by 2.0 subStepFactor = 2 !< for each substep, divide the last time increment by 2.0
real(pReal) :: & real(pReal) :: &
T_0 = 300.0_pReal, & t = 0.0_pReal, & !< elapsed time
time = 0.0_pReal, & !< elapsed time t_0 = 0.0_pReal, & !< begin of interval
time0 = 0.0_pReal, & !< begin of interval Delta_t = 1.0_pReal, & !< current time interval
timeinc = 1.0_pReal, & !< current time interval Delta_t_prev = 0.0_pReal, & !< previous time interval
timeIncOld = 0.0_pReal, & !< previous time interval t_remaining = 0.0_pReal !< remaining time of current load case
remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case
logical :: & logical :: &
guess, & !< guess along former trajectory guess, & !< guess along former trajectory
stagIterate, & stagIterate, &
@ -230,12 +229,6 @@ program DAMASK_grid
reportAndCheck: if (worldrank == 0) then reportAndCheck: if (worldrank == 0) then
print'(/,a,i0)', ' load case: ', l print'(/,a,i0)', ' load case: ', l
print*, ' estimate_rate:', loadCases(l)%estimate_rate 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 if (loadCases(l)%deformation%myType == 'F') then
print*, ' F:' print*, ' F:'
else else
@ -243,14 +236,14 @@ program DAMASK_grid
endif endif
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if (loadCases(l)%deformation%mask(i,j)) then 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 ' write(IO_STDOUT,'(2x,12a)',advance='no') ' x '
else
write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j)
endif endif
enddo; write(IO_STDOUT,'(/)',advance='no') enddo; write(IO_STDOUT,'(/)',advance='no')
enddo enddo
if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831 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 errorID = 838 ! no rotation is allowed by stress BC
if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:' if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:'
@ -259,9 +252,9 @@ program DAMASK_grid
if (loadCases(l)%stress%myType /= '') then if (loadCases(l)%stress%myType /= '') then
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if (loadCases(l)%stress%mask(i,j)) then 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 ' 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 endif
enddo; write(IO_STDOUT,'(/)',advance='no') enddo; write(IO_STDOUT,'(/)',advance='no')
enddo enddo
@ -303,7 +296,7 @@ program DAMASK_grid
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict) initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict)
thermal => initial_conditions%get('thermal',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) case(FIELD_DAMAGE_ID)
call grid_damage_spectral_init call grid_damage_spectral_init
@ -330,7 +323,7 @@ program DAMASK_grid
endif writeUndeformed endif writeUndeformed
loadCaseLooping: do l = 1, size(loadCases) 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 guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1, loadCases(l)%N incLooping: do inc = 1, loadCases(l)%N
@ -338,31 +331,31 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forwarding time ! 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 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 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) / (1.0_pReal-loadCases(l)%r**loadCases(l)%N)
endif 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? 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 guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping else skipping
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(l)%t+time0 - time t_remaining = loadCases(l)%t + t_0 - t
time = time + timeinc ! forward target time t = t + Delta_t ! forward target time
stepFraction = stepFraction + 1 ! count step stepFraction = stepFraction + 1 ! count step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
print'(/,a)', ' ###########################################################################' print'(/,a)', ' ###########################################################################'
print'(1x,a,es12.5,6(a,i0))', & print'(1x,a,es12.5,6(a,i0))', &
'Time', time, & 'Time', t, &
's: Increment ', inc,'/',loadCases(l)%N,& 's: Increment ', inc,'/',loadCases(l)%N,&
'-', stepFraction,'/',subStepFactor**cutBackLevel,& '-', stepFraction,'/',subStepFactor**cutBackLevel,&
' of load case ', l,'/',size(loadCases) ' of load case ', l,'/',size(loadCases)
@ -377,7 +370,7 @@ program DAMASK_grid
select case(ID(field)) select case(ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call mechanical_forward (& call mechanical_forward (&
cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & cutBack,guess,Delta_t,Delta_t_prev,t_remaining, &
deformation_BC = loadCases(l)%deformation, & deformation_BC = loadCases(l)%deformation, &
stress_BC = loadCases(l)%stress, & stress_BC = loadCases(l)%stress, &
rotation_BC = loadCases(l)%rot) rotation_BC = loadCases(l)%rot)
@ -398,9 +391,9 @@ program DAMASK_grid
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
solres(field) = mechanical_solution(incInfo) solres(field) = mechanical_solution(incInfo)
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
solres(field) = grid_thermal_spectral_solution(timeinc) solres(field) = grid_thermal_spectral_solution(Delta_t)
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
solres(field) = grid_damage_spectral_solution(timeinc) solres(field) = grid_damage_spectral_solution(Delta_t)
end select end select
if (.not. solres(field)%converged) exit ! no solution found if (.not. solres(field)%converged) exit ! no solution found
@ -418,11 +411,11 @@ program DAMASK_grid
if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged
.and. .not. solres(1)%termIll) then ! and acceptable solution found .and. .not. solres(1)%termIll) then ! and acceptable solution found
call mechanical_updateCoords call mechanical_updateCoords
timeIncOld = timeinc Delta_t_prev = Delta_t
cutBack = .false. cutBack = .false.
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
if (worldrank == 0) then if (worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, & write(statUnit,*) totalIncsCounter, t, cutBackLevel, &
solres(1)%converged, solres(1)%iterationsNeeded solres(1)%converged, solres(1)%iterationsNeeded
flush(statUnit) flush(statUnit)
endif endif
@ -430,8 +423,8 @@ program DAMASK_grid
cutBack = .true. cutBack = .true.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1 cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time t = t - Delta_t
timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep Delta_t = Delta_t/real(subStepFactor,pReal) ! cut timestep
print'(/,a)', ' cutting back ' print'(/,a)', ' cutting back '
else ! no more options to continue else ! no more options to continue
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
@ -453,7 +446,7 @@ program DAMASK_grid
if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then
print'(1/,a)', ' ... writing results to file ......................................' print'(1/,a)', ' ... writing results to file ......................................'
flush(IO_STDOUT) flush(IO_STDOUT)
call CPFEM_results(totalIncsCounter,time) call CPFEM_results(totalIncsCounter,t)
endif endif
if(signal) call interface_setSIGUSR1(.false.) if(signal) call interface_setSIGUSR1(.false.)
call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) 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 do i = 1,3
row => tensor%get(i) row => tensor%get(i)
do j = 1,3 do j = 1,3
mask(i,j) = row%get_asString(j) /= 'x' mask(i,j) = row%get_asString(j) == 'x'
if (mask(i,j)) values(i,j) = row%get_asFloat(j) if (.not. mask(i,j)) values(i,j) = row%get_asFloat(j)
enddo enddo
enddo enddo

10
src/grid/grid_damage_spectral.f90
View File

@ -160,10 +160,10 @@ end subroutine grid_damage_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral damage scheme with internal iterations !> @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) :: & real(pReal), intent(in) :: &
timeinc !< increment in time for current solution Delta_t !< increment in time for current solution
integer :: i, j, k, ce integer :: i, j, k, ce
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
@ -176,7 +176,7 @@ function grid_damage_spectral_solution(timeinc) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! 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 SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,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 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 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)) & + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
+ mu_ref*phi_current(i,j,k) + mu_ref*phi_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -292,7 +292,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator ! convolution of damage field with green operator
call utilities_FFTscalarForward 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 call utilities_FFTscalarBackward
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) & where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) &

20
src/grid/grid_mech_FEM.f90
View File

@ -324,7 +324,7 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution)
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. 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 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 else
C_volAvgLastInc = C_volAvg 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 F_aim_lastInc = F_aim
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! 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 & 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 elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot & 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 elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & 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 endif
if (guess) then 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 ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if (stress_BC%myType=='P') P_aim = P_aim & if (stress_BC%myType=='P') P_aim = P_aim &
+ merge((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 & 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) 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 ! set module wide available data
params%stress_mask = stress_BC%mask params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%Delta_t = Delta_t
end subroutine grid_mechanical_FEM_forward end subroutine grid_mechanical_FEM_forward
@ -568,13 +568,13 @@ subroutine formResidual(da_local,x_local, &
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(P_current,& call utilities_constitutiveResponse(P_current,&
P_av,C_volAvg,devNull, & 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) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc 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 ! constructing residual

20
src/grid/grid_mech_spectral_basic.f90
View File

@ -277,7 +277,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. 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 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_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg 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 F_aim_lastInc = F_aim
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! 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 & 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 elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot & 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 elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & 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 endif
Fdot = utilities_calculateRate(guess, & 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 ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if (stress_BC%myType=='P') P_aim = P_aim & if (stress_BC%myType=='P') P_aim = P_aim &
+ merge((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 & 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 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]) 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 ! set module wide available data
params%stress_mask = stress_BC%mask params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%Delta_t = Delta_t
end subroutine grid_mechanical_spectral_basic_forward end subroutine grid_mechanical_spectral_basic_forward
@ -492,14 +492,14 @@ subroutine formResidual(in, F, &
! evaluate constitutive response ! 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, & 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) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
F_aim = F_aim - deltaF_aim 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 ! updated deformation gradient using fix point algorithm of basic scheme

22
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -309,7 +309,7 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. 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 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_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg 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 F_aim_lastInc = F_aim
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! 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 & 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 elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot & 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 elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & 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 endif
Fdot = utilities_calculateRate(guess, & 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 ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if(stress_BC%myType=='P') P_aim = P_aim & if(stress_BC%myType=='P') P_aim = P_aim &
+ merge((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 & 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 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.)),& 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 ! set module wide available data
params%stress_mask = stress_BC%mask params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%Delta_t = Delta_t
end subroutine grid_mechanical_spectral_polarisation_forward end subroutine grid_mechanical_spectral_polarisation_forward
@ -592,14 +592,14 @@ subroutine formResidual(in, FandF_tau, &
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, & 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) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc 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, & err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), &
math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),& P_av-P_aim, &
params%stress_mask))) params%stress_mask)))
! calculate divergence ! calculate divergence
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal

14
src/grid/grid_thermal_spectral.f90
View File

@ -155,10 +155,10 @@ end subroutine grid_thermal_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral thermal scheme with internal iterations !> @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) :: & real(pReal), intent(in) :: &
timeinc !< increment in time for current solution Delta_t !< increment in time for current solution
integer :: i, j, k, ce integer :: i, j, k, ce
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
@ -171,7 +171,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! 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 SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(thermal_snes,reason,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 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 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 enddo; enddo; enddo
call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
@ -233,7 +233,7 @@ subroutine grid_thermal_spectral_forward(cutBack)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 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 enddo; enddo; enddo
else else
T_lastInc = T_current T_lastInc = T_current
@ -279,7 +279,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 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)) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
+ mu_ref*T_current(i,j,k) + mu_ref*T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -287,7 +287,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convolution of temperature field with green operator ! convolution of temperature field with green operator
call utilities_FFTscalarForward 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 call utilities_FFTscalarBackward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

28
src/grid/spectral_utilities.f90
View File

@ -88,7 +88,7 @@ module spectral_utilities
type, public :: tBoundaryCondition !< set of parameters defining a boundary condition type, public :: tBoundaryCondition !< set of parameters defining a boundary condition
real(pReal), dimension(3,3) :: values = 0.0_pReal 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 character(len=:), allocatable :: myType
end type tBoundaryCondition end type tBoundaryCondition
@ -96,7 +96,7 @@ module spectral_utilities
real(pReal), dimension(3,3) :: stress_BC real(pReal), dimension(3,3) :: stress_BC
logical, dimension(3,3) :: stress_mask logical, dimension(3,3) :: stress_mask
type(rotation) :: rotation_BC type(rotation) :: rotation_BC
real(pReal) :: timeinc real(pReal) :: Delta_t
end type tSolutionParams end type tSolutionParams
type :: tNumerics type :: tNumerics
@ -684,7 +684,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
logical :: errmatinv logical :: errmatinv
character(len=pStringLen):: formatString character(len=pStringLen):: formatString
mask_stressVector = reshape(transpose(mask_stress), [9]) mask_stressVector = .not. reshape(transpose(mask_stress), [9])
size_reduced = count(mask_stressVector) size_reduced = count(mask_stressVector)
if(size_reduced > 0) then if(size_reduced > 0) then
temp99_real = math_3333to99(rot_BC%rotate(C)) 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,& 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,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) :: P_av !< average PK stress
real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< 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), 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 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 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) & 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) & 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 = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt 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) pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
avRate !< homogeneous addon avRate !< homogeneous addon
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
dt !< timeinc between field0 and field dt !< Delta_t between field0 and field
logical, intent(in) :: & logical, intent(in) :: &
heterogeneous !< calculate field of rates heterogeneous !< calculate field of rates
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & 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, !> @brief forwards a field with a pointwise given rate, if aim is given,
!> ensures that the average matches the aim !> 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) :: & 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) :: & real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
field_lastInc, & !< initial field field_lastInc, & !< initial field
rate !< rate by which to forward rate !< rate by which to forward
@ -913,7 +913,7 @@ function utilities_forwardField(timeinc,field_lastInc,rate,aim)
real(pReal), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim real(pReal), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
PetscErrorCode :: ierr PetscErrorCode :: ierr
utilities_forwardField = field_lastInc + rate*timeinc utilities_forwardField = field_lastInc + rate*Delta_t
if (present(aim)) then !< correct to match average if (present(aim)) then !< correct to match average
fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt 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) call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)