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

use nomenclature from the DAMASK paper

This commit is contained in:
Martin Diehl 2020-01-30 23:07:45 +01:00
parent 348a91d503
commit 269d65005b
2 changed files with 128 additions and 135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

119
src/grid/grid_damage_spectral.f90
View File

@ -29,15 +29,15 @@ module grid_damage_spectral
Vec, private :: solution_vec
PetscInt, private :: xstart, xend, ystart, yend, zstart, zend
real(pReal), private, dimension(:,:,:), allocatable :: &
damage_current, & !< field of current damage
damage_lastInc, & !< field of previous damage
damage_stagInc !< field of staggered damage
phi_current, & !< field of current damage
phi_lastInc, & !< field of previous damage
phi_stagInc !< field of staggered damage
!--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc.
integer, private :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3), private :: D_ref
real(pReal), private :: mobility_ref
real(pReal), dimension(3,3), private :: K_ref
real(pReal), private :: mu_ref
public :: &
grid_damage_spectral_init, &
@ -55,7 +55,6 @@ contains
subroutine grid_damage_spectral_init
PetscInt, dimension(0:worldsize-1) :: localK
integer :: i, j, k, cell
DM :: damage_grid
Vec :: uBound, lBound
PetscErrorCode :: ierr
@ -116,25 +115,14 @@ subroutine grid_damage_spectral_init
xend = xstart + xend - 1
yend = ystart + yend - 1
zend = zstart + zend - 1
allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! damage reference diffusion update
cell = 0
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell)
mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
call updateReference
end subroutine grid_damage_spectral_init
@ -150,7 +138,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
integer :: i, j, k, cell
type(tSolutionState) :: solution
PetscInt :: devNull
PetscReal :: minDamage, maxDamage, stagNorm, solnNorm
PetscReal :: phi_min, phi_max, stagNorm, solnNorm
PetscErrorCode :: ierr
SNESConvergedReason :: reason
@ -172,11 +160,11 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
solution%converged = .true.
solution%iterationsNeeded = totalIter
endif
stagNorm = maxval(abs(damage_current - damage_stagInc))
solnNorm = maxval(abs(damage_current))
stagNorm = maxval(abs(phi_current - phi_stagInc))
solnNorm = maxval(abs(phi_current))
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
damage_stagInc = damage_current
phi_stagInc = phi_current
solution%stagConverged = stagNorm < max(err_damage_tolAbs, err_damage_tolRel*solnNorm)
!--------------------------------------------------------------------------------------------------
@ -184,15 +172,15 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
enddo; enddo; enddo
call VecMin(solution_vec,devNull,minDamage,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,devNull,maxDamage,ierr); CHKERRQ(ierr)
call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr)
if (solution%converged) &
write(6,'(/,a)') ' ... nonlocal damage converged .....................................'
write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',&
minDamage, maxDamage, stagNorm
phi_min, phi_max, stagNorm
write(6,'(/,a)') ' ==========================================================================='
flush(6)
@ -211,35 +199,22 @@ subroutine grid_damage_spectral_forward(cutBack)
PetscErrorCode :: ierr
if (cutBack) then
damage_current = damage_lastInc
damage_stagInc = damage_lastInc
phi_current = phi_lastInc
phi_stagInc = phi_lastInc
!--------------------------------------------------------------------------------------------------
! reverting damage field state
cell = 0
call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current
x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
enddo; enddo; enddo
else
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
damage_lastInc = damage_current
cell = 0
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell)
mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
phi_lastInc = phi_current
call updateReference
endif
end subroutine grid_damage_spectral_forward
@ -263,18 +238,18 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
integer :: i, j, k, cell
real(pReal) :: phiDot, dPhiDot_dPhi, mobility
damage_current = x_scal
phi_current = x_scal
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current
scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current
call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion33(1,cell) - K_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward
@ -283,29 +258,51 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell)
call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k), 1, cell)
mobility = damage_nonlocal_getMobility(1,cell)
scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + &
params%timeinc*phiDot + &
mobility*damage_lastInc(i,j,k) - &
mobility*damage_current(i,j,k) + &
mobility_ref*damage_current(i,j,k)
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) &
+ mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
+ mu_ref*phi_current(i,j,k)
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator
call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc)
call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc)
call utilities_FFTscalarBackward
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current
f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current
end subroutine formResidual
!--------------------------------------------------------------------------------------------------
!> @brief update reference viscosity and conductivity
!--------------------------------------------------------------------------------------------------
subroutine updateReference
integer :: i,j,k,cell,ierr
cell = 0
K_ref = 0.0_pReal
mu_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
K_ref = K_ref + damage_nonlocal_getDiffusion33(1,cell)
mu_ref = mu_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo
K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mu_ref = mu_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mu_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine updateReference
end module grid_damage_spectral

144
src/grid/grid_thermal_spectral.f90
View File

@ -30,15 +30,15 @@ module grid_thermal_spectral
Vec, private :: solution_vec
PetscInt, private :: xstart, xend, ystart, yend, zstart, zend
real(pReal), private, dimension(:,:,:), allocatable :: &
temperature_current, & !< field of current temperature
temperature_lastInc, & !< field of previous temperature
temperature_stagInc !< field of staggered temperature
T_current, & !< field of current temperature
T_lastInc, & !< field of previous temperature
T_stagInc !< field of staggered temperature
!--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc.
integer, private :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3), private :: D_ref
real(pReal), private :: mobility_ref
real(pReal), dimension(3,3), private :: K_ref
real(pReal), private :: mu_ref
public :: &
grid_thermal_spectral_init, &
@ -104,36 +104,22 @@ subroutine grid_thermal_spectral_init
xend = xstart + xend - 1
yend = ystart + yend - 1
zend = zstart + zend - 1
allocate(temperature_current(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(temperature_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(temperature_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
temperature_current(i,j,k) = temperature(material_homogenizationAt(cell))% &
T_current(i,j,k) = temperature(material_homogenizationAt(cell))% &
p(thermalMapping(material_homogenizationAt(cell))%p(1,cell))
temperature_lastInc(i,j,k) = temperature_current(i,j,k)
temperature_stagInc(i,j,k) = temperature_current(i,j,k)
T_lastInc(i,j,k) = T_current(i,j,k)
T_stagInc(i,j,k) = T_current(i,j,k)
enddo; enddo; enddo
call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current
x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current
call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! thermal reference diffusion update
cell = 0
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
D_ref = D_ref + thermal_conduction_getConductivity33(1,cell)
mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
call updateReference
end subroutine grid_thermal_spectral_init
@ -148,8 +134,8 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
timeinc_old !< increment in time of last increment
integer :: i, j, k, cell
type(tSolutionState) :: solution
PetscInt :: position
PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm
PetscInt :: devNull
PetscReal :: T_min, T_max, stagNorm, solnNorm
PetscErrorCode :: ierr
SNESConvergedReason :: reason
@ -171,11 +157,11 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
solution%converged = .true.
solution%iterationsNeeded = totalIter
endif
stagNorm = maxval(abs(temperature_current - temperature_stagInc))
solnNorm = maxval(abs(temperature_current))
stagNorm = maxval(abs(T_current - T_stagInc))
solnNorm = maxval(abs(T_current))
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
temperature_stagInc = temperature_current
T_stagInc = T_current
solution%stagConverged = stagNorm < max(err_thermal_tolAbs, err_thermal_tolRel*solnNorm)
!--------------------------------------------------------------------------------------------------
@ -183,17 +169,17 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), &
(temperature_current(i,j,k)-temperature_lastInc(i,j,k))/params%timeinc, &
call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), &
(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, &
1,cell)
enddo; enddo; enddo
call VecMin(solution_vec,position,minTemperature,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,position,maxTemperature,ierr); CHKERRQ(ierr)
call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr)
if (solution%converged) &
write(6,'(/,a)') ' ... thermal conduction converged ..................................'
write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',&
minTemperature, maxTemperature, stagNorm
T_min, T_max, stagNorm
write(6,'(/,a)') ' ==========================================================================='
flush(6)
@ -212,40 +198,26 @@ subroutine grid_thermal_spectral_forward(cutBack)
PetscErrorCode :: ierr
if (cutBack) then
temperature_current = temperature_lastInc
temperature_stagInc = temperature_lastInc
T_current = T_lastInc
T_stagInc = T_lastInc
!--------------------------------------------------------------------------------------------------
! reverting thermal field state
cell = 0
call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current
x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), &
(temperature_current(i,j,k) - &
temperature_lastInc(i,j,k))/params%timeinc, &
call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), &
(T_current(i,j,k) - &
T_lastInc(i,j,k))/params%timeinc, &
1,cell)
enddo; enddo; enddo
else
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
temperature_lastInc = temperature_current
cell = 0
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
D_ref = D_ref + thermal_conduction_getConductivity33(1,cell)
mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)
enddo; enddo; enddo
D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
T_lastInc = T_current
call updateReference
endif
end subroutine grid_thermal_spectral_forward
@ -269,45 +241,69 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
integer :: i, j, k, cell
real(pReal) :: Tdot, dTdot_dT
temperature_current = x_scal
T_current = x_scal
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current
scalarField_real(1:grid(1),1:grid(2),1:grid3) = T_current
call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_fourierScalarGradient !< calculate gradient of temperature field
call utilities_FFTvectorBackward
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity33(1,cell) - K_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field
call utilities_FFTscalarBackward
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell)
scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + &
params%timeinc*Tdot + &
thermal_conduction_getMassDensity (1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)*(temperature_lastInc(i,j,k) - &
temperature_current(i,j,k)) + &
mobility_ref*temperature_current(i,j,k)
call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T_current(i,j,k), 1, cell)
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
+ thermal_conduction_getMassDensity (1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)*(T_lastInc(i,j,k) - &
T_current(i,j,k))&
+ mu_ref*T_current(i,j,k)
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator
! convolution of temperature field with green operator
call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc)
call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc)
call utilities_FFTscalarBackward
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3)
f_scal = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3)
end subroutine formResidual
!--------------------------------------------------------------------------------------------------
!> @brief update reference viscosity and conductivity
!--------------------------------------------------------------------------------------------------
subroutine updateReference
integer :: i,j,k,cell,ierr
cell = 0
K_ref = 0.0_pReal
mu_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
K_ref = K_ref + thermal_conduction_getConductivity33(1,cell)
mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)
enddo; enddo; enddo
K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mu_ref = mu_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mu_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine updateReference
end module grid_thermal_spectral