diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 9d804ec56..c0a1b0f40 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -19,7 +19,7 @@ module grid_thermal_spectral use YAML_types use config use material - + implicit none private @@ -45,11 +45,11 @@ module grid_thermal_spectral T_stagInc !< field of staggered temperature !-------------------------------------------------------------------------------------------------- -! reference diffusion tensor, mobility etc. +! reference diffusion tensor, mobility etc. integer :: totalIter = 0 !< total iteration in current increment real(pReal), dimension(3,3) :: K_ref real(pReal) :: mu_ref - + public :: & grid_thermal_spectral_init, & grid_thermal_spectral_solution, & @@ -63,8 +63,8 @@ contains !-------------------------------------------------------------------------------------------------- subroutine grid_thermal_spectral_init - PetscInt, dimension(0:worldsize-1) :: localK - integer :: i, j, k, cell + PetscInt, dimension(0:worldsize-1) :: localK + integer :: i, j, k, ce DM :: thermal_grid PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr @@ -93,11 +93,11 @@ subroutine grid_thermal_spectral_init CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) CHKERRQ(ierr) - + !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) + call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) localK = 0 localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) @@ -115,23 +115,23 @@ subroutine grid_thermal_spectral_init call DMsetUp(thermal_grid,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(thermal_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) + CHKERRQ(ierr) call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- -! init fields +! init fields call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,ierr) CHKERRQ(ierr) xend = xstart + xend - 1 yend = ystart + yend - 1 - zend = zstart + zend - 1 + zend = zstart + zend - 1 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 + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - T_current(i,j,k) = temperature(material_homogenizationAt(cell))%p(material_homogenizationMemberAt(1,cell)) + ce = ce + 1 + T_current(i,j,k) = temperature(material_homogenizationAt(ce))%p(material_homogenizationMemberAt(1,ce)) T_lastInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k) enddo; enddo; enddo @@ -143,26 +143,26 @@ subroutine grid_thermal_spectral_init end subroutine grid_thermal_spectral_init - + !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_thermal_spectral_solution(timeinc) result(solution) - + real(pReal), intent(in) :: & timeinc !< increment in time for current solution - integer :: i, j, k, cell + integer :: i, j, k, ce type(tSolutionState) :: solution PetscInt :: devNull PetscReal :: T_min, T_max, stagNorm, solnNorm - PetscErrorCode :: ierr + PetscErrorCode :: ierr SNESConvergedReason :: reason solution%converged =.false. - + !-------------------------------------------------------------------------------------------------- -! set module wide availabe data +! set module wide availabe data params%timeinc = timeinc call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) @@ -183,13 +183,13 @@ function grid_thermal_spectral_solution(timeinc) result(solution) solution%stagConverged = stagNorm < max(num%eps_thermal_atol, num%eps_thermal_rtol*solnNorm) !-------------------------------------------------------------------------------------------------- -! updating thermal state - cell = 0 +! updating thermal state + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 + ce = ce + 1 call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & - 1,cell) + 1,ce) enddo; enddo; enddo call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) @@ -198,7 +198,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution) print'(/,a)', ' ... thermal conduction converged ..................................' print'(/,a,f8.4,2x,f8.4,2x,f8.4)', ' Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm print'(/,a)', ' ===========================================================================' - flush(IO_STDOUT) + flush(IO_STDOUT) end function grid_thermal_spectral_solution @@ -207,36 +207,36 @@ end function grid_thermal_spectral_solution !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- subroutine grid_thermal_spectral_forward(cutBack) - + logical, intent(in) :: cutBack - integer :: i, j, k, cell + integer :: i, j, k, ce DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr - - if (cutBack) then + + if (cutBack) then T_current = T_lastInc T_stagInc = T_lastInc !-------------------------------------------------------------------------------------------------- -! reverting thermal field state - cell = 0 +! reverting thermal field state + ce = 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) = 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 + ce = ce + 1 call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & (T_current(i,j,k) - & T_lastInc(i,j,k))/params%timeinc, & - 1,cell) + 1,ce) enddo; enddo; enddo else T_lastInc = T_current call updateReference endif - + end subroutine grid_thermal_spectral_forward @@ -244,7 +244,7 @@ end subroutine grid_thermal_spectral_forward !> @brief forms the spectral thermal residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,f_scal,dummy,ierr) - + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & @@ -255,33 +255,33 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) f_scal PetscObject :: dummy PetscErrorCode :: ierr - integer :: i, j, k, cell + integer :: i, j, k, ce real(pReal) :: Tdot - T_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) = T_current + scalarField_real(1:grid(1),1:grid(2),1:grid3) = T_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of temperature field call utilities_FFTvectorBackward - cell = 0 + ce = 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_getConductivity(1,cell) - K_ref, & + ce = ce + 1 + vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(1,ce) - K_ref, & vectorField_real(1:3,i,j,k)) enddo; enddo; enddo call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field call utilities_FFTscalarBackward - cell = 0 + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call thermal_conduction_getSource(Tdot, T_current(i,j,k), 1, cell) + ce = ce + 1 + call thermal_conduction_getSource(Tdot, T_current(i,j,k), 1, ce) 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) - & + + thermal_conduction_getMassDensity (1,ce)* & + thermal_conduction_getSpecificHeat(1,ce)*(T_lastInc(i,j,k) - & T_current(i,j,k))& + mu_ref*T_current(i,j,k) enddo; enddo; enddo @@ -291,7 +291,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) call utilities_FFTscalarForward call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc) call utilities_FFTscalarBackward - + !-------------------------------------------------------------------------------------------------- ! constructing residual f_scal = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) @@ -304,15 +304,15 @@ end subroutine formResidual !-------------------------------------------------------------------------------------------------- subroutine updateReference - integer :: i,j,k,cell,ierr - - cell = 0 + integer :: i,j,k,ce,ierr + + ce = 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_getConductivity(1,cell) - mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* thermal_conduction_getSpecificHeat(1,cell) + ce = ce + 1 + K_ref = K_ref + thermal_conduction_getConductivity(1,ce) + mu_ref = mu_ref + thermal_conduction_getMassDensity(1,ce)* thermal_conduction_getSpecificHeat(1,ce) 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)