diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index d7bbbdefe..b5f283b30 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -43,8 +43,8 @@ program DAMASK_spectral_Driver debug_levelBasic use math ! need to include the whole module for FFTW use mesh, only: & - gridGlobal, & - geomSizeGlobal + grid, & + geomSize use CPFEM, only: & CPFEM_initAll use FEsolving, only: & @@ -390,8 +390,8 @@ program DAMASK_spectral_Driver write(resUnit) 'load:', trim(loadCaseFile) ! ... and write header write(resUnit) 'workingdir:', trim(getSolverWorkingDirectoryName()) write(resUnit) 'geometry:', trim(geometryFile) - write(resUnit) 'grid:', gridGlobal - write(resUnit) 'size:', geomSizeGlobal + write(resUnit) 'grid:', grid + write(resUnit) 'size:', geomSize write(resUnit) 'materialpoint_sizeResults:', materialpoint_sizeResults write(resUnit) 'loadcases:', size(loadCases) write(resUnit) 'frequencies:', loadCases%outputfrequency ! one entry per LoadCase diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index 6eea15bd9..17e4a5748 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -12,10 +12,10 @@ module DAMASK_spectral_solverAL pReal use math, only: & math_I3 - use DAMASK_spectral_utilities, only: & + use DAMASK_spectral_Utilities, only: & tSolutionState, & tSolutionParams - + implicit none private #include @@ -121,8 +121,8 @@ subroutine AL_init Utilities_updateGamma, & Utilities_updateIPcoords use mesh, only: & - gridLocal, & - gridGlobal + grid, & + grid3 use math, only: & math_invSym3333 @@ -145,36 +145,35 @@ subroutine AL_init #include "compilation_info.f90" endif - allocate (P (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) + allocate (P (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) - allocate (Fdot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) - allocate (F_lambda_lastInc(3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) - allocate (F_lambdaDot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) + allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_lambda_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_lambdaDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! PETSc Init call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 do proc = 1, worldsize call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) enddo call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid + grid(1),grid(2),grid(3), & ! global grid 1 , 1, worldsize, & 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) - gridLocal (1),gridLocal (2),localK, & ! local grid + grid(1),grid(2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,dummy,ierr) CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr) CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) @@ -184,13 +183,7 @@ subroutine AL_init call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data F => xx_psc(0:8,:,:,:) F_lambda => xx_psc(9:17,:,:,:) - if (restartInc == 1_pInt) then ! no deformation (no restart) - F_lastInc = spread(spread(spread(math_I3,3,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) ! initialize to identity - F = reshape(F_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)]) - F_lambda = F - F_lambda_lastInc = F_lastInc - - elseif (restartInc > 1_pInt) then ! read in F to calculate coordinates and initialize CPFEM general + restart: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading values of increment ', restartInc - 1_pInt, ' from file' @@ -218,16 +211,22 @@ subroutine AL_init call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) read (777,rec=1) f_aimDot close (777) - endif + elseif (restartInc == 1_pInt) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_lambda = F + F_lambda_lastInc = F_lastInc + endif restart + - call utilities_updateIPcoords(F) + call Utilities_updateIPcoords(F) call Utilities_constitutiveResponse(F_lastInc,F,& 0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) nullify(F) nullify(F_lambda) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc - if (restartInc > 1_pInt) then ! using old values from files + readRestart: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading more values of increment', restartInc - 1_pInt, 'from file' @@ -241,7 +240,7 @@ subroutine AL_init call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) read (777,rec=1) C_minMaxAvg close (777) - endif + endif readRestart call Utilities_updateGamma(C_minMaxAvg,.True.) C_scale = C_minMaxAvg @@ -299,7 +298,7 @@ type(tSolutionState) function & C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) endif - + AL_solution%converged =.false. !-------------------------------------------------------------------------------------------------- @@ -339,7 +338,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) polarBeta, & worldrank use mesh, only: & - gridLocal + grid3, & + grid use IO, only: & IO_intOut use math, only: & @@ -350,7 +350,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) math_mul33x33 use DAMASK_spectral_Utilities, only: & wgt, & - tensorField_realMPI, & + tensorField_real, & utilities_FFTtensorForward, & utilities_fourierGammaConvolution, & Utilities_inverseLaplace, & @@ -408,7 +408,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - if(totalIter <= PETScIter) then ! new iteration + newIteration: if(totalIter <= PETScIter) then !-------------------------------------------------------------------------------------------------- ! report begin of new iteration totalIter = totalIter + 1_pInt @@ -420,15 +420,15 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & math_transpose33(F_aim) + flush(6) endif - flush(6) - endif + endif newIteration !-------------------------------------------------------------------------------------------------- ! - tensorField_realMPI = 0.0_pReal - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) - tensorField_realMPI(1:3,1:3,i,j,k) = & + tensorField_real = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + tensorField_real(1:3,1:3,i,j,k) = & polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3)) @@ -443,7 +443,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_lambda = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) + residual_F_lambda = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -455,8 +455,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculate divergence - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F call utilities_FFTtensorForward() err_div = Utilities_divergenceRMS() call utilities_FFTtensorBackward() @@ -464,7 +464,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! constructing residual e = 0_pInt - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) e = e + 1_pInt residual_F(1:3,1:3,i,j,k) = math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & residual_F(1:3,1:3,i,j,k) - & @@ -475,8 +475,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculating curl - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F call utilities_FFTtensorForward() err_curl = Utilities_curlRMS() call utilities_FFTtensorBackward() @@ -570,7 +570,8 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ use numerics, only: & worldrank use mesh, only: & - gridLocal + grid3, & + grid use DAMASK_spectral_Utilities, only: & Utilities_calculateRate, & Utilities_forwardField, & @@ -639,14 +640,14 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ write (777,rec=1) C_volAvgLastInc close(777) endif + endif - endif call utilities_updateIPcoords(F) if (cutBack) then F_aim = F_aim_lastInc - F_lambda = reshape(F_lambda_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)]) - F = reshape(F_lastInc, [9,gridLocal(1),gridLocal(2),gridLocal(3)]) + F_lambda = reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid3]) + F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) C_volAvg = C_volAvgLastInc else ForwardData = .True. @@ -667,22 +668,25 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ ! update coordinates and rate and forward last inc call utilities_updateIPcoords(F) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_lastInc,reshape(F,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)])) + timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3])) F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,gridLocal(1), & - gridLocal(2),gridLocal(3)])) - F_lastInc = reshape(F, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)]) - F_lambda_lastInc = reshape(F_lambda,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]) + timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1), & + grid(2),grid3])) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) + F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid3]) endif F_aim = F_aim + f_aimDot * timeinc + +!-------------------------------------------------------------------------------------------------- +! update local deformation gradient F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim math_rotate_backward33(F_aim,rotation_BC)), & - [9,gridLocal(1),gridLocal(2),gridLocal(3)]) + [9,grid(1),grid(2),grid3]) F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & - [9,gridLocal(1),gridLocal(2),gridLocal(3)]) ! does not have any average value as boundary condition + [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition if (.not. guess) then ! large strain forwarding - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) F_lambda33 = reshape(F_lambda(1:9,i,j,k),[3,3]) F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & math_mul3333xx33(C_scale,& diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 5651ff00e..5b7e0b14c 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -94,9 +94,9 @@ subroutine basicPETSc_init IO_read_realFile, & IO_timeStamp use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRestart + debug_level, & + debug_spectral, & + debug_spectralRestart use FEsolving, only: & restartInc use numerics, only: & @@ -110,8 +110,8 @@ subroutine basicPETSc_init utilities_updateIPcoords, & wgt use mesh, only: & - gridLocal, & - gridGlobal + grid, & + grid3 use math, only: & math_invSym3333 @@ -133,27 +133,27 @@ subroutine basicPETSc_init #include "compilation_info.f90" endif mainProcess - allocate (P (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) + allocate (P (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) - allocate (Fdot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) + allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 do proc = 1, worldsize call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) enddo call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid + grid(1),grid(2),grid(3), & ! global grid 1, 1, worldsize, & 9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) - gridLocal (1),gridLocal (2),localK, & ! local grid + grid (1),grid (2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) @@ -169,10 +169,7 @@ subroutine basicPETSc_init ! init fields call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with - if (restartInc == 1_pInt) then ! no deformation (no restart) - F_lastInc = spread(spread(spread(math_I3,3,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) ! initialize to identity - F = reshape(F_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)]) - elseif (restartInc > 1_pInt) then ! using old values from file + restart: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading values of increment ', restartInc - 1_pInt, ' from file' @@ -189,7 +186,10 @@ subroutine basicPETSc_init close (777) F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc - endif + elseif (restartInc == 1_pInt) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + endif restart call Utilities_updateIPcoords(F) call Utilities_constitutiveResponse(F_lastInc, F, & @@ -202,7 +202,7 @@ subroutine basicPETSc_init call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc - if (restartInc > 1_pInt) then ! using old values from files + restartRead: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading more values of increment', restartInc - 1_pInt, 'from file' @@ -216,7 +216,7 @@ subroutine basicPETSc_init call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) read (777,rec=1) C_minMaxAvg close (777) - endif + endif restartRead call Utilities_updateGamma(C_minmaxAvg,.True.) @@ -302,7 +302,8 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) use numerics, only: & worldrank use mesh, only: & - gridLocal + grid, & + grid3 use math, only: & math_rotate_backward33, & math_transpose33, & @@ -312,9 +313,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) debug_spectral, & debug_spectralRotation use DAMASK_spectral_Utilities, only: & - wgt, & - tensorField_realMPI, & - tensorField_fourierMPI, & + tensorField_real, & utilities_FFTtensorForward, & utilities_FFTtensorBackward, & utilities_fourierGammaConvolution, & @@ -345,7 +344,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - if (totalIter <= PETScIter) then ! new iteration + newIteration: if (totalIter <= PETScIter) then !-------------------------------------------------------------------------------------------------- ! report begin of new iteration totalIter = totalIter + 1_pInt @@ -353,13 +352,13 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & - math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & + math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & math_transpose33(F_aim) - endif flush(6) - endif + endif + endif newIteration !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -376,8 +375,8 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = f_scal + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = f_scal call utilities_FFTtensorForward() err_div = Utilities_divergenceRMS() call utilities_fourierGammaConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC)) @@ -385,7 +384,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) + f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) end subroutine BasicPETSc_formResidual @@ -455,7 +454,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r math_mul33x33 ,& math_rotate_backward33 use mesh, only: & - gridLocal + grid, & + grid3 use DAMASK_spectral_Utilities, only: & Utilities_calculateRate, & Utilities_forwardField, & @@ -513,9 +513,10 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r endif call utilities_updateIPcoords(F) + if (cutBack) then F_aim = F_aim_lastInc - F = reshape(F_lastInc, [9,gridLocal(1),gridLocal(2),gridLocal(3)]) + F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) C_volAvg = C_volAvgLastInc else ForwardData = .True. @@ -536,15 +537,15 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r ! update coordinates and rate and forward last inc call utilities_updateIPcoords(F) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_lastInc,reshape(F,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)])) - F_lastInc = reshape(F, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)]) + timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3])) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) endif F_aim = F_aim + f_aimDot * timeinc !-------------------------------------------------------------------------------------------------- ! update local deformation gradient F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim - math_rotate_backward33(F_aim,rotation_BC)),[9,gridLocal(1),gridLocal(2),gridLocal(3)]) + math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) end subroutine BasicPETSc_forward diff --git a/code/DAMASK_spectral_solverPolarisation.f90 b/code/DAMASK_spectral_solverPolarisation.f90 index 4ae57842c..ca4e1a556 100644 --- a/code/DAMASK_spectral_solverPolarisation.f90 +++ b/code/DAMASK_spectral_solverPolarisation.f90 @@ -12,10 +12,10 @@ module DAMASK_spectral_solverPolarisation pReal use math, only: & math_I3 - use DAMASK_spectral_utilities, only: & + use DAMASK_spectral_Utilities, only: & tSolutionState, & tSolutionParams - + implicit none private #include @@ -105,7 +105,7 @@ subroutine Polarisation_init IO_intOut, & IO_read_realFile, & IO_timeStamp - use debug, only : & + use debug, only: & debug_level, & debug_spectral, & debug_spectralRestart @@ -121,8 +121,8 @@ subroutine Polarisation_init Utilities_updateGamma, & Utilities_updateIPcoords use mesh, only: & - gridLocal, & - gridGlobal + grid, & + grid3 use math, only: & math_invSym3333 @@ -145,29 +145,29 @@ subroutine Polarisation_init #include "compilation_info.f90" endif - allocate (P (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) + allocate (P (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) - allocate (Fdot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) - allocate (F_tau_lastInc(3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) - allocate (F_tauDot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal) + allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! PETSc Init call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 do proc = 1, worldsize call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) enddo call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid + grid(1),grid(2),grid(3), & ! global grid 1 , 1, worldsize, & 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) - gridLocal (1),gridLocal (2),localK, & ! local grid + grid (1),grid (2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) @@ -183,12 +183,7 @@ subroutine Polarisation_init call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data F => xx_psc(0:8,:,:,:) F_tau => xx_psc(9:17,:,:,:) - if (restartInc == 1_pInt) then ! no deformation (no restart) - F_lastInc = spread(spread(spread(math_I3,3,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) ! initialize to identity - F = reshape(F_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)]) - F_tau = 2.0_pReal* F - F_tau_lastInc = 2.0_pReal*F_lastInc - elseif (restartInc > 1_pInt) then + restart: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading values of increment', restartInc - 1_pInt, 'from file' @@ -216,7 +211,12 @@ subroutine Polarisation_init call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) read (777,rec=1) f_aimDot close (777) - endif + elseif (restartInc == 1_pInt) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_tau = 2.0_pReal* F + F_tau_lastInc = 2.0_pReal*F_lastInc + endif restart call Utilities_updateIPcoords(F) call Utilities_constitutiveResponse(F_lastInc,F,& @@ -225,7 +225,7 @@ subroutine Polarisation_init nullify(F_tau) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc - if (restartInc > 1_pInt) then ! using old values from files + readRestart: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading more values of increment', restartInc - 1_pInt, 'from file' @@ -239,7 +239,7 @@ subroutine Polarisation_init call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) read (777,rec=1) C_minMaxAvg close (777) - endif + endif readRestart call Utilities_updateGamma(C_minMaxAvg,.True.) C_scale = C_minMaxAvg @@ -300,7 +300,7 @@ type(tSolutionState) function & Polarisation_solution%converged =.false. !-------------------------------------------------------------------------------------------------- -! set module wide availabe data +! set module wide availabe data mask_stress = P_BC%maskFloat params%P_BC = P_BC%values params%rotation_BC = rotation_BC @@ -335,10 +335,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) polarAlpha, & polarBeta, & worldrank + use mesh, only: & + grid3, & + grid use IO, only: & IO_intOut - use mesh, only: & - gridLocal use math, only: & math_rotate_backward33, & math_transpose33, & @@ -347,7 +348,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) math_mul33x33 use DAMASK_spectral_Utilities, only: & wgt, & - tensorField_realMPI, & + tensorField_real, & utilities_FFTtensorForward, & utilities_fourierGammaConvolution, & Utilities_inverseLaplace, & @@ -405,7 +406,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - if (totalIter <= PETScIter) then ! new iteration + newIteration: if(totalIter <= PETScIter) then !-------------------------------------------------------------------------------------------------- ! report begin of new iteration totalIter = totalIter + 1_pInt @@ -413,19 +414,19 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & - math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & + math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & math_transpose33(F_aim) - endif - flush(6) - endif + flush(6) + endif + endif newIteration !-------------------------------------------------------------------------------------------------- ! - tensorField_realMPI = 0.0_pReal - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) - tensorField_realMPI(1:3,1:3,i,j,k) = & + tensorField_real = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + tensorField_real(1:3,1:3,i,j,k) = & polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) @@ -439,7 +440,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_tau = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) + residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -451,8 +452,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculate divergence - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F call utilities_FFTtensorForward() err_div = Utilities_divergenceRMS() call utilities_FFTtensorBackward() @@ -460,7 +461,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! constructing residual e = 0_pInt - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) e = e + 1_pInt residual_F(1:3,1:3,i,j,k) = & math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & @@ -471,8 +472,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculating curl - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F call utilities_FFTtensorForward() err_curl = Utilities_curlRMS() call utilities_FFTtensorBackward() @@ -491,8 +492,8 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, err_div_tolAbs, & err_curl_tolRel, & err_curl_tolAbs, & - err_stress_tolabs, & - err_stress_tolrel, & + err_stress_tolAbs, & + err_stress_tolRel, & worldrank use math, only: & math_mul3333xx33 @@ -563,20 +564,21 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC math_mul3333xx33, & math_transpose33, & math_rotate_backward33 + use numerics, only: & + worldrank + use mesh, only: & + grid3, & + grid use DAMASK_spectral_Utilities, only: & Utilities_calculateRate, & Utilities_forwardField, & Utilities_updateIPcoords, & tBoundaryCondition, & cutBack - use mesh, only: & - gridLocal use IO, only: & IO_write_JobRealFile use FEsolving, only: & restartWrite - use numerics, only: & - worldrank implicit none real(pReal), intent(in) :: & @@ -636,11 +638,10 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC endif call utilities_updateIPcoords(F) - if (cutBack) then F_aim = F_aim_lastInc - F_tau= reshape(F_tau_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)]) - F = reshape(F_lastInc, [9,gridLocal(1),gridLocal(2),gridLocal(3)]) + F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid3]) + F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) C_volAvg = C_volAvgLastInc else ForwardData = .True. @@ -662,24 +663,25 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC call utilities_updateIPcoords(F) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & timeinc_old,guess,F_lastInc, & - reshape(F,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)])) + reshape(F,[3,3,grid(1),grid(2),grid3])) F_tauDot = Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), & timeinc_old,guess,F_tau_lastInc, & - reshape(F_tau,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)])) - F_lastInc = reshape(F, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)]) - F_tau_lastInc = reshape(F_tau,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]) + reshape(F_tau,[3,3,grid(1),grid(2),grid3])) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) + F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) endif + F_aim = F_aim + f_aimDot * timeinc !-------------------------------------------------------------------------------------------------- ! update local deformation gradient - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim math_rotate_backward33(F_aim,rotation_BC)), & - [9,gridLocal(1),gridLocal(2),gridLocal(3)]) + [9,grid(1),grid(2),grid3]) F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & ! does not have any average value as boundary condition - [9,gridLocal(1),gridLocal(2),gridLocal(3)]) + [9,grid(1),grid(2),grid3]) if (.not. guess) then ! large strain forwarding - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & math_mul3333xx33(C_scale,& @@ -699,7 +701,7 @@ end subroutine Polarisation_forward subroutine Polarisation_destroy() use DAMASK_spectral_Utilities, only: & Utilities_destroy - + implicit none PetscErrorCode :: ierr diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index a399736f0..b31d835a2 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -42,12 +42,12 @@ module DAMASK_spectral_utilities !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW integer(pInt), public :: grid1Red !< grid(1)/2 - real (C_DOUBLE), public, dimension(:,:,:,:,:), pointer :: tensorField_realMPI !< real representation (some stress or deformation) of field_fourier - complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:,:), pointer :: tensorField_fourierMPI !< field on which the Fourier transform operates - real(C_DOUBLE), public, dimension(:,:,:,:), pointer :: vectorField_realMPI !< vector field real representation for fftw - complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:), pointer :: vectorField_fourierMPI !< vector field fourier representation for fftw - real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_realMPI !< scalar field real representation for fftw - complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:), pointer :: scalarField_fourierMPI !< scalar field fourier representation for fftw + real (C_DOUBLE), public, dimension(:,:,:,:,:), pointer :: tensorField_real !< real representation (some stress or deformation) of field_fourier + complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:,:), pointer :: tensorField_fourier !< field on which the Fourier transform operates + real(C_DOUBLE), public, dimension(:,:,:,:), pointer :: vectorField_real !< vector field real representation for fftw + complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field fourier representation for fftw + real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for fftw + complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field fourier representation for fftw real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method real(pReal), private, dimension(:,:,:,:), allocatable :: xi !< wave vector field for divergence and for gamma operator real(pReal), private, dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness @@ -56,22 +56,17 @@ module DAMASK_spectral_utilities !-------------------------------------------------------------------------------------------------- ! plans for FFTW type(C_PTR), private :: & - planTensorForth, & !< FFTW MPI plan P(x) to P(k) - planTensorBack, & !< FFTW MPI plan F(k) to F(x) - planVectorForth, & !< FFTW MPI plan P(x) to P(k) - planVectorBack, & !< FFTW MPI plan F(k) to F(x) - planScalarForth, & !< FFTW MPI plan P(x) to P(k) - planScalarBack, & !< FFTW MPI plan F(k) to F(x) - planDebugForth, & !< FFTW MPI plan for scalar field - planDebugBack, & !< FFTW MPI plan for scalar field inverse - planDiv !< FFTW MPI plan for FFTW in case of debugging divergence calculation + planTensorForth, & !< FFTW MPI plan P(x) to P(k) + planTensorBack, & !< FFTW MPI plan F(k) to F(x) + planVectorForth, & !< FFTW MPI plan v(x) to v(k) + planVectorBack, & !< FFTW MPI plan v(k) to v(x) + planScalarForth, & !< FFTW MPI plan s(x) to s(k) + planScalarBack !< FFTW MPI plan s(k) to s(x) !-------------------------------------------------------------------------------------------------- ! variables controlling debugging logical, private :: & debugGeneral, & !< general debugging of spectral solver - debugDivergence, & !< debugging of divergence calculation (comparison to function used for post processing) - debugFFTW, & !< doing additional FFT on scalar field and compare to results of strided 3D FFT debugRotation, & !< also printing out results in lab frame debugPETSc !< use some in debug defined options for more verbose PETSc solution @@ -190,10 +185,10 @@ subroutine utilities_init() #endif use math use mesh, only: & - gridGlobal, & - gridLocal, & - gridOffset, & - geomSizeGlobal + grid, & + grid3, & + grid3Offset, & + geomSize implicit none #ifdef PETSc @@ -206,9 +201,9 @@ subroutine utilities_init() integer(pInt) :: i, j, k integer(pInt), dimension(3) :: k_s type(C_PTR) :: & - tensorFieldMPI, & !< field cotaining data for FFTW in real and fourier space (in place) - vectorFieldMPI, & !< field cotaining data for FFTW in real space when debugging FFTW (no in place) - scalarFieldMPI !< field cotaining data for FFTW in real space when debugging FFTW (no in place) + tensorField, & !< field cotaining data for FFTW in real and fourier space (in place) + vectorField, & !< field cotaining data for FFTW in real space when debugging FFTW (no in place) + scalarField !< field cotaining data for FFTW in real space when debugging FFTW (no in place) integer(C_INTPTR_T) :: gridFFTW(3), alloc_local, local_K, local_K_offset integer(C_INTPTR_T), parameter :: & scalarSize = 1_C_INTPTR_T, & @@ -225,8 +220,6 @@ subroutine utilities_init() !-------------------------------------------------------------------------------------------------- ! set debugging parameters debugGeneral = iand(debug_level(debug_SPECTRAL),debug_LEVELBASIC) /= 0 - debugDivergence = iand(debug_level(debug_SPECTRAL),debug_SPECTRALDIVERGENCE) /= 0 - debugFFTW = iand(debug_level(debug_SPECTRAL),debug_SPECTRALFFTW) /= 0 debugRotation = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION) /= 0 debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0 @@ -240,117 +233,96 @@ subroutine utilities_init() call PetscOptionsInsertString(trim(petsc_defaultOptions),ierr); CHKERRQ(ierr) call PetscOptionsInsertString(trim(petsc_options),ierr); CHKERRQ(ierr) - grid1Red = gridLocal(1)/2_pInt + 1_pInt - wgt = 1.0/real(product(gridGlobal),pReal) + grid1Red = grid(1)/2_pInt + 1_pInt + wgt = 1.0/real(product(grid),pReal) if (worldrank == 0) then - write(6,'(a,3(i12 ))') ' grid a b c: ', gridGlobal - write(6,'(a,3(es12.5))') ' size x y z: ', geomSizeGlobal + write(6,'(a,3(i12 ))') ' grid a b c: ', grid + write(6,'(a,3(es12.5))') ' size x y z: ', geomSize endif !-------------------------------------------------------------------------------------------------- -! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso- -! lution-independent divergence +! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and +! resolution-independent divergence if (divergence_correction == 1_pInt) then do j = 1_pInt, 3_pInt - if (j /= minloc(geomSizeGlobal,1) .and. j /= maxloc(geomSizeGlobal,1)) & - scaledGeomSize = geomSizeGlobal/geomSizeGlobal(j) + if (j /= minloc(geomSize,1) .and. j /= maxloc(geomSize,1)) & + scaledGeomSize = geomSize/geomSize(j) enddo elseif (divergence_correction == 2_pInt) then do j = 1_pInt, 3_pInt - if (j /= minloc(geomSizeGlobal/gridGlobal,1) .and. j /= maxloc(geomSizeGlobal/gridGlobal,1)) & - scaledGeomSize = geomSizeGlobal/geomSizeGlobal(j)*gridGlobal(j) + if (j /= minloc(geomSize/grid,1) .and. j /= maxloc(geomSize/grid,1)) & + scaledGeomSize = geomSize/geomSize(j)*grid(j) enddo else - scaledGeomSize = geomSizeGlobal + scaledGeomSize = geomSize endif !-------------------------------------------------------------------------------------------------- ! MPI allocation - gridFFTW = int(gridGlobal,C_INTPTR_T) + gridFFTW = int(grid,C_INTPTR_T) alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, & MPI_COMM_WORLD, local_K, local_K_offset) + allocate (xi(3,grid1Red,grid(2),grid3),source = 0.0_pReal) ! frequencies, only half the size for first dimension - tensorFieldMPI = fftw_alloc_complex(tensorSize*alloc_local) - call c_f_pointer(tensorFieldMPI, tensorField_realMPI, [3_C_INTPTR_T,3_C_INTPTR_T, & + tensorField = fftw_alloc_complex(tensorSize*alloc_local) + call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, & 2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K]) ! place a pointer for a real tensor representation - call c_f_pointer(tensorFieldMPI, tensorField_fourierMPI, [3_C_INTPTR_T,3_C_INTPTR_T, & + call c_f_pointer(tensorField, tensorField_fourier, [3_C_INTPTR_T,3_C_INTPTR_T, & gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T , gridFFTW(2),local_K]) ! place a pointer for a fourier tensor representation - vectorFieldMPI = fftw_alloc_complex(vecSize*alloc_local) - call c_f_pointer(vectorFieldMPI, vectorField_realMPI, [3_C_INTPTR_T,& + vectorField = fftw_alloc_complex(vecSize*alloc_local) + call c_f_pointer(vectorField, vectorField_real, [3_C_INTPTR_T,& 2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K]) ! place a pointer for a real vector representation - call c_f_pointer(vectorFieldMPI, vectorField_fourierMPI,[3_C_INTPTR_T,& + call c_f_pointer(vectorField, vectorField_fourier,[3_C_INTPTR_T,& gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T, gridFFTW(2),local_K]) ! place a pointer for a fourier vector representation - scalarFieldMPI = fftw_alloc_complex(scalarSize*alloc_local) ! allocate data for real representation (no in place transform) - call c_f_pointer(scalarFieldMPI, scalarField_realMPI, & + scalarField = fftw_alloc_complex(scalarSize*alloc_local) ! allocate data for real representation (no in place transform) + call c_f_pointer(scalarField, scalarField_real, & [2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1),gridFFTW(2),local_K]) ! place a pointer for a real scalar representation - call c_f_pointer(scalarFieldMPI, scalarField_fourierMPI, & + call c_f_pointer(scalarField, scalarField_fourier, & [ gridFFTW(1)/2_C_INTPTR_T + 1 ,gridFFTW(2),local_K]) ! place a pointer for a fourier scarlar representation - allocate (xi(3,grid1Red,gridLocal(2),gridLocal(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension !-------------------------------------------------------------------------------------------------- ! tensor MPI fftw plans - planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order + planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock - tensorField_realMPI, tensorField_fourierMPI, &! input data, output data + tensorField_real, tensorField_fourier, & ! input data, output data MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision if (.not. C_ASSOCIATED(planTensorForth)) call IO_error(810, ext_msg='planTensorForth') - planTensorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order + planTensorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock - tensorField_fourierMPI,tensorField_realMPI, &! input data, output data + tensorField_fourier,tensorField_real, & ! input data, output data MPI_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision if (.not. C_ASSOCIATED(planTensorBack)) call IO_error(810, ext_msg='planTensorBack') !-------------------------------------------------------------------------------------------------- ! vector MPI fftw plans - planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order + planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock - vectorField_realMPI, vectorField_fourierMPI, &! input data, output data + vectorField_real, vectorField_fourier, & ! input data, output data MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision if (.not. C_ASSOCIATED(planVectorForth)) call IO_error(810, ext_msg='planVectorForth') - planVectorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order + planVectorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock - vectorField_fourierMPI,vectorField_realMPI, & ! input data, output data + vectorField_fourier,vectorField_real, & ! input data, output data MPI_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision if (.not. C_ASSOCIATED(planVectorBack)) call IO_error(810, ext_msg='planVectorBack') !-------------------------------------------------------------------------------------------------- ! scalar MPI fftw plans - planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order + planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock - scalarField_realMPI, scalarField_fourierMPI, & ! input data, output data + scalarField_real, scalarField_fourier, & ! input data, output data MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision if (.not. C_ASSOCIATED(planScalarForth)) call IO_error(810, ext_msg='planScalarForth') planScalarBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order, no. of transforms scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock - scalarField_fourierMPI,scalarField_realMPI, & ! input data, output data + scalarField_fourier,scalarField_real, & ! input data, output data MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision if (.not. C_ASSOCIATED(planScalarBack)) call IO_error(810, ext_msg='planScalarBack') -!-------------------------------------------------------------------------------------------------- -! depending on debug options, allocate more memory and create additional plans - if (debugDivergence) then - planDiv = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2) ,gridFFTW(1)],vecSize, & - FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & - vectorField_fourierMPI, vectorField_realMPI, & - MPI_COMM_WORLD, fftw_planner_flag) - if (.not. C_ASSOCIATED(planDiv)) call IO_error(810, ext_msg='planDiv') - endif - if (debugFFTW) then - planDebugForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & - scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & - scalarField_realMPI, scalarField_fourierMPI, & - MPI_COMM_WORLD, fftw_planner_flag) - if (.not. C_ASSOCIATED(planDebugForth)) call IO_error(810, ext_msg='planDebugForth') - planDebugBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & - scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & - scalarField_fourierMPI,scalarField_realMPI, & - MPI_COMM_WORLD, fftw_planner_flag) - if (.not. C_ASSOCIATED(planDebugBack)) call IO_error(810, ext_msg='planDebugBack') - endif !-------------------------------------------------------------------------------------------------- ! general initialization of FFTW (see manual on fftw.org for more details) if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(0_pInt,ext_msg='Fortran to C') ! check for correct precision in C @@ -361,21 +333,21 @@ subroutine utilities_init() !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - do k = gridOffset+1_pInt, gridOffset+gridLocal(3) + do k = grid3Offset+1_pInt, grid3Offset+grid3 k_s(3) = k - 1_pInt - if(k > gridGlobal(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - gridGlobal(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 - do j = 1_pInt, gridLocal(2) + if(k > grid(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + do j = 1_pInt, grid(2) k_s(2) = j - 1_pInt - if(j > gridGlobal(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - gridGlobal(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + if(j > grid(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 do i = 1_pInt, grid1Red k_s(1) = i - 1_pInt ! symmetry, junst running from 0,1,...,N/2,N/2+1 - xi(1:3,i,j,k-gridOffset) = real(k_s, pReal)/scaledGeomSize ! if divergence_correction is set, frequencies are calculated on unit length + xi(1:3,i,j,k-grid3Offset) = real(k_s, pReal)/scaledGeomSize ! if divergence_correction is set, frequencies are calculated on unit length enddo; enddo; enddo if(memory_efficient) then ! allocate just single fourth order tensor allocate (gamma_hat(3,3,3,3,1,1,1), source = 0.0_pReal) else ! precalculation of gamma_hat field - allocate (gamma_hat(3,3,3,3,grid1Red,gridLocal(2),gridLocal(3)), source = 0.0_pReal) + allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = 0.0_pReal) endif end subroutine utilities_init @@ -396,8 +368,9 @@ subroutine utilities_updateGamma(C,saveReference) memory_efficient, & worldrank use mesh, only: & - gridOffset, & - gridLocal + grid3Offset, & + grid3,& + grid use math, only: & math_inv33 @@ -421,15 +394,15 @@ subroutine utilities_updateGamma(C,saveReference) endif if(.not. memory_efficient) then - do k = gridOffset+1_pInt, gridOffset+gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, grid1Red + do k = grid3Offset+1_pInt, grid3Offset+grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red if (any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - xiDyad(l,m) = xi(l, i,j,k-gridOffset)*xi(m, i,j,k-gridOffset) + xiDyad(l,m) = xi(l, i,j,k-grid3Offset)*xi(m, i,j,k-grid3Offset) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & temp33_Real(l,m) = sum(C_ref(l,1:3,m,1:3)*xiDyad) temp33_Real = math_inv33(temp33_Real) forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)& - gamma_hat(l,m,n,o,i,j,k-gridOffset) = temp33_Real(l,n)*xiDyad(m,o) + gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_Real(l,n)*xiDyad(m,o) endif enddo; enddo; enddo endif @@ -438,72 +411,25 @@ end subroutine utilities_updateGamma !-------------------------------------------------------------------------------------------------- !> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed -!> @details Does an unweighted FFT transform from real to complex. -!> In case of debugging the FFT, also one component of the tensor (specified by row and column) -!> is independetly transformed complex to complex and compared to the whole tensor transform +!> @details Does an unweighted filtered FFT transform from real to complex. !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTtensorForward() - use math - use numerics, only: & - worldrank use mesh, only: & - gridLocal + grid3, & + grid implicit none - integer(pInt) :: row, column ! if debug FFTW, compare 3D array field of row and column - real(pReal), dimension(2) :: myRand, maxScalarField ! random numbers integer(pInt) :: i, j, k - PetscErrorCode :: ierr - -!-------------------------------------------------------------------------------------------------- -! copy one component of the stress field to to a single FT and check for mismatch - if (debugFFTW) then - if (worldrank == 0_pInt) then - call random_number(myRand) ! two numbers: 0 <= x < 1 - row = nint(myRand(1)*2_pReal + 1_pReal,pInt) - column = nint(myRand(2)*2_pReal + 1_pReal,pInt) - endif - call MPI_Bcast(row ,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) - call MPI_Bcast(column,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) - scalarField_realMPI = 0.0_pReal - scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = & - tensorField_realMPI(row,column,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) ! store the selected component - endif !-------------------------------------------------------------------------------------------------- ! doing the FFT - call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_realMPI,tensorField_fourierMPI) + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) -!-------------------------------------------------------------------------------------------------- -! comparing 1 and 3x3 FT results - if (debugFFTW) then - call fftw_mpi_execute_dft_r2c(planDebugForth,scalarField_realMPI,scalarField_fourierMPI) - where(abs(scalarField_fourierMPI(1:grid1Red,1:gridLocal(2),1:gridLocal(3))) > tiny(1.0_pReal)) ! avoid division by zero - scalarField_fourierMPI(1:grid1Red,1:gridLocal(2),1:gridLocal(3)) = & - (scalarField_fourierMPI(1:grid1Red,1:gridLocal(2),1:gridLocal(3))-& - tensorField_fourierMPI(row,column,1:grid1Red,1:gridLocal(2),1:gridLocal(3)))/& - scalarField_fourierMPI(1:grid1Red,1:gridLocal(2),1:gridLocal(3)) - else where - scalarField_realMPI = cmplx(0.0,0.0,pReal) - end where - maxScalarField(1) = maxval(real (scalarField_fourierMPI(1:grid1Red,1:gridLocal(2), & - 1:gridLocal(3)))) - maxScalarField(2) = maxval(aimag(scalarField_fourierMPI(1:grid1Red,1:gridLocal(2), & - 1:gridLocal(3)))) - call MPI_reduce(MPI_IN_PLACE,maxScalarField,2,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr) - if (worldrank == 0_pInt) then - write(6,'(/,a,i1,1x,i1,a)') ' .. checking FT results of compontent ', row, column, ' ..' - write(6,'(/,a,2(es11.4,1x))') ' max FT relative error = ',& ! print real and imaginary part seperately - maxScalarField(1),maxScalarField(2) - flush(6) - endif - endif - !-------------------------------------------------------------------------------------------------- ! applying filter - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,grid1Red - tensorField_fourierMPI(1:3,1:3,i,j,k) = utilities_getFilter(xi(1:3,i,j,k))* & - tensorField_fourierMPI(1:3,1:3,i,j,k) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red + tensorField_fourier(1:3,1:3,i,j,k) = utilities_getFilter(xi(1:3,i,j,k))* & + tensorField_fourier(1:3,1:3,i,j,k) enddo; enddo; enddo end subroutine utilities_FFTtensorForward @@ -511,146 +437,87 @@ end subroutine utilities_FFTtensorForward !-------------------------------------------------------------------------------------------------- !> @brief backward FFT of data in field_fourier to field_real -!> @details Does an inverse FFT transform from complex to real -!> In case of debugging the FFT, also one component of the tensor (specified by row and column) -!> is independetly transformed complex to complex and compared to the whole tensor transform -!> results is weighted by number of points stored in wgt +!> @details Does an weighted inverse FFT transform from complex to real !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTtensorBackward() - use math - use numerics, only: & - worldrank - use mesh, only: & - gridLocal - implicit none - integer(pInt) :: row, column !< if debug FFTW, compare 3D array field of row and column - real(pReal), dimension(2) :: myRand - real(pReal) :: maxScalarField - PetscErrorCode :: ierr -!-------------------------------------------------------------------------------------------------- -! unpack FFT data for conj complex symmetric part. This data is not transformed when using c2r - if (debugFFTW) then - if (worldrank == 0_pInt) then - call random_number(myRand) ! two numbers: 0 <= x < 1 - row = nint(myRand(1)*2_pReal + 1_pReal,pInt) - column = nint(myRand(2)*2_pReal + 1_pReal,pInt) - endif - call MPI_Bcast(row ,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) - call MPI_Bcast(column,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) - scalarField_fourierMPI(1:grid1Red,1:gridLocal(2),1:gridLocal(3)) = & - tensorField_fourierMPI(row,column,1:grid1Red,1:gridLocal(2),1:gridLocal(3)) - endif - -!-------------------------------------------------------------------------------------------------- -! doing the iFFT - call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourierMPI,tensorField_realMPI) ! back transform of fluct deformation gradient - -!-------------------------------------------------------------------------------------------------- -! comparing 1 and 3x3 inverse FT results - if (debugFFTW) then - call fftw_mpi_execute_dft_c2r(planDebugBack,scalarField_fourierMPI,scalarField_realMPI) - where(abs(real(scalarField_realMPI,pReal)) > tiny(1.0_pReal)) ! avoid division by zero - scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = & - (scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) & - - tensorField_realMPI (row,column,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)))/ & - scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) - else where - scalarField_realMPI = cmplx(0.0,0.0,pReal) - end where - maxScalarField = maxval(real (scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)))) - call MPI_reduce(MPI_IN_PLACE,maxScalarField,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr) - if (worldrank == 0_pInt) then - write(6,'(/,a,i1,1x,i1,a)') ' ... checking iFT results of compontent ', row, column, ' ..' - write(6,'(/,a,es11.4)') ' max iFT relative error = ', maxScalarField - flush(6) - endif - endif - - tensorField_realMPI = tensorField_realMPI * wgt ! normalize the result by number of elements + call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) ! back transform of fluct deformation gradient + tensorField_real = tensorField_real * wgt ! normalize the result by number of elements end subroutine utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- -!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed -!> @details Does an unweighted FFT transform from real to complex. -!> In case of debugging the FFT, also one component of the scalar -!> is independetly transformed complex to complex and compared to the whole scalar transform +!> @brief forward FFT of data in scalarField_real to scalarField_fourier +!> @details Does an unweighted filtered FFT transform from real to complex. !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTscalarForward() - use math use mesh, only: & - gridLocal - + grid3, & + grid + + implicit none integer(pInt) :: i, j, k +!-------------------------------------------------------------------------------------------------- ! doing the scalar FFT - call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_realMPI,scalarField_fourierMPI) + call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) +!-------------------------------------------------------------------------------------------------- ! applying filter - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,grid1Red - scalarField_fourierMPI(i,j,k) = utilities_getFilter(xi(1:3,i,j,k))* & - scalarField_fourierMPI(i,j,k) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red + scalarField_fourier(i,j,k) = utilities_getFilter(xi(1:3,i,j,k))* & + scalarField_fourier(i,j,k) enddo; enddo; enddo end subroutine utilities_FFTscalarForward !-------------------------------------------------------------------------------------------------- -!> @brief backward FFT of data in field_fourier to field_real -!> @details Does an inverse FFT transform from complex to real -!> In case of debugging the FFT, also one component of the scalar -!> is independetly transformed complex to complex and compared to the whole scalar transform -!> results is weighted by number of points stored in wgt +!> @brief backward FFT of data in scalarField_fourier to scalarField_real +!> @details Does an weighted inverse FFT transform from complex to real !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTscalarBackward() - use math - -! doing the scalar iFFT - call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourierMPI,scalarField_realMPI) + implicit none - scalarField_realMPI = scalarField_realMPI * wgt ! normalize the result by number of elements + call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) + + scalarField_real = scalarField_real * wgt ! normalize the result by number of elements end subroutine utilities_FFTscalarBackward !-------------------------------------------------------------------------------------------------- !> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed -!> @details Does an unweighted FFT transform from real to complex. -!> In case of debugging the FFT, also one component of the vector -!> is independetly transformed complex to complex and compared to the whole vector transform +!> @details Does an unweighted filtered FFT transform from real to complex. !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTvectorForward() - use math use mesh, only: & - gridLocal - + grid3, & + grid + + implicit none integer(pInt) :: i, j, k -! doing the vecotr FFT - call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_realMPI,vectorField_fourierMPI) +! doing the vector FFT + call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) ! applying filter - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,grid1Red - vectorField_fourierMPI(1:3,i,j,k) = utilities_getFilter(xi(1:3,i,j,k))* & - vectorField_fourierMPI(1:3,i,j,k) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red + vectorField_fourier(1:3,i,j,k) = utilities_getFilter(xi(1:3,i,j,k))* & + vectorField_fourier(1:3,i,j,k) enddo; enddo; enddo end subroutine utilities_FFTvectorForward !-------------------------------------------------------------------------------------------------- !> @brief backward FFT of data in field_fourier to field_real -!> @details Does an inverse FFT transform from complex to real -!> In case of debugging the FFT, also one component of the vector -!> is independetly transformed complex to complex and compared to the whole vector transform -!> results is weighted by number of points stored in wgt +!> @details Does an weighted inverse FFT transform from complex to real !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTvectorBackward() - use math - -! doing the vector iFFT - call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourierMPI,vectorField_realMPI) + implicit none - vectorField_realMPI = vectorField_realMPI * wgt ! normalize the result by number of elements + call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) + + vectorField_real = vectorField_real * wgt ! normalize the result by number of elements end subroutine utilities_FFTvectorBackward @@ -664,9 +531,10 @@ subroutine utilities_inverseLaplace() use numerics, only: & worldrank use mesh, only: & - gridLocal, & - gridOffset, & - geomSizeGlobal + grid, & + grid3, & + grid3Offset, & + geomSize implicit none integer(pInt) :: i, j, k @@ -677,16 +545,16 @@ subroutine utilities_inverseLaplace() flush(6) endif - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, grid1Red + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red k_s = xi(1:3,i,j,k)*scaledGeomSize - if (any(k_s /= 0_pInt)) tensorField_fourierMPI(1:3,1:3,i,j,k-gridOffset) = & - tensorField_fourierMPI(1:3,1:3,i,j,k-gridOffset)/ & - cmplx(-sum((2.0_pReal*PI*k_s/geomSizeGlobal)* & - (2.0_pReal*PI*k_s/geomSizeGlobal)),0.0_pReal,pReal) + if (any(k_s /= 0_pInt)) tensorField_fourier(1:3,1:3,i,j,k-grid3Offset) = & + tensorField_fourier(1:3,1:3,i,j,k-grid3Offset)/ & + cmplx(-sum((2.0_pReal*PI*k_s/geomSize)* & + (2.0_pReal*PI*k_s/geomSize)),0.0_pReal,pReal) enddo; enddo; enddo - if (gridOffset == 0_pInt) & - tensorField_fourierMPI(1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal) + if (grid3Offset == 0_pInt) & + tensorField_fourier(1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal) end subroutine utilities_inverseLaplace @@ -702,8 +570,9 @@ subroutine utilities_fourierGammaConvolution(fieldAim) use numerics, only: & worldrank use mesh, only: & - gridLocal, & - gridOffset + grid3, & + grid, & + grid3Offset implicit none real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution @@ -721,9 +590,9 @@ subroutine utilities_fourierGammaConvolution(fieldAim) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation (mechanical equilibrium) - if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2) ;do i = 1_pInt, grid1Red - if(any([i,j,k+gridOffset] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + memoryEfficient: if(memory_efficient) then + do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red + if(any([i,j,k+grid3Offset] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & @@ -733,21 +602,21 @@ subroutine utilities_fourierGammaConvolution(fieldAim) gamma_hat(l,m,n,o, 1,1,1) = temp33_Real(l,n)*xiDyad(m,o) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3, 1,1,1) * & - tensorField_fourierMPI(1:3,1:3,i,j,k)) - tensorField_fourierMPI(1:3,1:3,i,j,k) = temp33_Complex + tensorField_fourier(1:3,1:3,i,j,k)) + tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex endif enddo; enddo; enddo - else ! use precalculated gamma-operator - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,grid1Red + else memoryEfficient + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * & - tensorField_fourierMPI(1:3,1:3,i,j,k)) - tensorField_fourierMPI(1:3,1:3,i,j,k) = temp33_Complex + tensorField_fourier(1:3,1:3,i,j,k)) + tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex enddo; enddo; enddo - endif + endif memoryEfficient - if (gridOffset == 0_pInt) & - tensorField_fourierMPI(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + if (grid3Offset == 0_pInt) & + tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 end subroutine utilities_fourierGammaConvolution @@ -756,16 +625,13 @@ end subroutine utilities_fourierGammaConvolution !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT) - use numerics, only: & - memory_efficient use math, only: & math_mul33x3, & PI - use numerics, only: & - worldrank use mesh, only: & - gridLocal, & - geomSizeGlobal + grid, & + grid3, & + geomSize implicit none real(pReal), dimension(3,3), intent(in) :: D_ref !< desired average value of the field after convolution @@ -776,12 +642,12 @@ subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2) ;do i = 1_pInt, grid1Red + do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red k_s = xi(1:3,i,j,k)*scaledGeomSize GreenOp_hat = 1.0_pReal/ & - (mobility_ref + deltaT*sum((2.0_pReal*PI*k_s/geomSizeGlobal)* & - math_mul33x3(D_ref,(2.0_pReal*PI*k_s/geomSizeGlobal))))!< GreenOp_hat = iK^{T} * D_ref * iK, K is frequency - scalarField_fourierMPI(i,j,k) = scalarField_fourierMPI(i,j,k)*GreenOp_hat + (mobility_ref + deltaT*sum((2.0_pReal*PI*k_s/geomSize)* & + math_mul33x3(D_ref,(2.0_pReal*PI*k_s/geomSize)))) !< GreenOp_hat = iK^{T} * D_ref * iK, K is frequency + scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat enddo; enddo; enddo end subroutine utilities_fourierGreenConvolution @@ -789,21 +655,18 @@ end subroutine utilities_fourierGreenConvolution !-------------------------------------------------------------------------------------------------- !> @brief calculate root mean square of divergence of field_fourier !-------------------------------------------------------------------------------------------------- -real(pReal) function utilities_divergenceRMS() - use math +real(pReal) function utilities_divergenceRMS() + use math, only: & + TWOPIIMG, & + math_mul33x3_complex use numerics, only: & worldrank use mesh, only: & - gridLocal, & - gridGlobal + grid, & + grid3 implicit none integer(pInt) :: i, j, k - real(pReal) :: & - err_real_div_RMS, & !< RMS of divergence in real space - err_div_max, & !< maximum value of divergence in Fourier space - err_real_div_max !< maximum value of divergence in real space - complex(pReal), dimension(3) :: temp3_complex PetscErrorCode :: ierr @@ -815,60 +678,28 @@ real(pReal) function utilities_divergenceRMS() !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space utilities_divergenceRMS = 0.0_pReal - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2) do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. utilities_divergenceRMS = utilities_divergenceRMS & - + 2.0_pReal*(sum (real(math_mul33x3_complex(tensorField_fourierMPI(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again + + 2.0_pReal*(sum (real(math_mul33x3_complex(tensorField_fourier(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector - +sum(aimag(math_mul33x3_complex(tensorField_fourierMPI(1:3,1:3,i,j,k),& + +sum(aimag(math_mul33x3_complex(tensorField_fourier(1:3,1:3,i,j,k),& xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) enddo utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1) - + sum( real(math_mul33x3_complex(tensorField_fourierMPI(1:3,1:3,1 ,j,k), & + + sum( real(math_mul33x3_complex(tensorField_fourier(1:3,1:3,1 ,j,k), & xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal) & - + sum(aimag(math_mul33x3_complex(tensorField_fourierMPI(1:3,1:3,1 ,j,k), & + + sum(aimag(math_mul33x3_complex(tensorField_fourier(1:3,1:3,1 ,j,k), & xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal) & - + sum( real(math_mul33x3_complex(tensorField_fourierMPI(1:3,1:3,grid1Red,j,k), & + + sum( real(math_mul33x3_complex(tensorField_fourier(1:3,1:3,grid1Red,j,k), & xi(1:3,grid1Red,j,k))*TWOPIIMG)**2.0_pReal) & - + sum(aimag(math_mul33x3_complex(tensorField_fourierMPI(1:3,1:3,grid1Red,j,k), & + + sum(aimag(math_mul33x3_complex(tensorField_fourier(1:3,1:3,grid1Red,j,k), & xi(1:3,grid1Red,j,k))*TWOPIIMG)**2.0_pReal) enddo; enddo - if(gridGlobal(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space -!-------------------------------------------------------------------------------------------------- -! calculate additional divergence criteria and report - if (debugDivergence) then ! calculate divergence again - err_div_max = 0.0_pReal - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, grid1Red - temp3_Complex = math_mul33x3_complex(tensorField_fourierMPI(1:3,1:3,i,j,k)*wgt,& ! weighting P_fourier - xi(1:3,i,j,k))*TWOPIIMG - err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal)) - vectorField_fourierMPI(1:3,i,j,k) = temp3_Complex ! need divergence NOT squared - enddo; enddo; enddo - - call fftw_mpi_execute_dft_c2r(planDiv,vectorField_fourierMPI,vectorField_realMPI) ! already weighted - - err_real_div_RMS = sum(vectorField_realMPI**2.0_pReal) - call MPI_reduce(MPI_IN_PLACE,err_real_div_RMS,1,MPI_DOUBLE,MPI_SUM,0,PETSC_COMM_WORLD,ierr) - err_real_div_RMS = sqrt(wgt*err_real_div_RMS) ! RMS in real space - - err_real_div_max = maxval(sum(vectorField_realMPI**2.0_pReal,dim=4)) ! max in real space - call MPI_reduce(MPI_IN_PLACE,err_real_div_max,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr) - err_real_div_max = sqrt(err_real_div_max) - - call MPI_reduce(MPI_IN_PLACE,err_div_max,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr) - err_div_max = sqrt( err_div_max) ! max in Fourier space - - if (worldrank == 0_pInt) then - write(6,'(/,1x,a,es11.4)') 'error divergence FT RMS = ',utilities_divergenceRMS - write(6,'(1x,a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS - write(6,'(1x,a,es11.4)') 'error divergence FT max = ',err_div_max - write(6,'(1x,a,es11.4)') 'error divergence Real max = ',err_real_div_max - flush(6) - endif - endif end function utilities_divergenceRMS @@ -881,8 +712,8 @@ real(pReal) function utilities_curlRMS() use numerics, only: & worldrank use mesh, only: & - gridLocal, & - gridGlobal + grid, & + grid3 implicit none integer(pInt) :: i, j, k, l @@ -898,43 +729,43 @@ real(pReal) function utilities_curlRMS() ! calculating max curl criterion in Fourier space utilities_curlRMS = 0.0_pReal - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 2_pInt, grid1Red - 1_pInt do l = 1_pInt, 3_pInt - curl_fourier(l,1) = (+tensorField_fourierMPI(l,3,i,j,k)*xi(2,i,j,k)& - -tensorField_fourierMPI(l,2,i,j,k)*xi(3,i,j,k))*TWOPIIMG - curl_fourier(l,2) = (+tensorField_fourierMPI(l,1,i,j,k)*xi(3,i,j,k)& - -tensorField_fourierMPI(l,3,i,j,k)*xi(1,i,j,k))*TWOPIIMG - curl_fourier(l,3) = (+tensorField_fourierMPI(l,2,i,j,k)*xi(1,i,j,k)& - -tensorField_fourierMPI(l,1,i,j,k)*xi(2,i,j,k))*TWOPIIMG + curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi(2,i,j,k)& + -tensorField_fourier(l,2,i,j,k)*xi(3,i,j,k))*TWOPIIMG + curl_fourier(l,2) = (+tensorField_fourier(l,1,i,j,k)*xi(3,i,j,k)& + -tensorField_fourier(l,3,i,j,k)*xi(1,i,j,k))*TWOPIIMG + curl_fourier(l,3) = (+tensorField_fourier(l,2,i,j,k)*xi(1,i,j,k)& + -tensorField_fourier(l,1,i,j,k)*xi(2,i,j,k))*TWOPIIMG enddo utilities_curlRMS = utilities_curlRMS + & 2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) enddo do l = 1_pInt, 3_pInt - curl_fourier = (+tensorField_fourierMPI(l,3,1,j,k)*xi(2,1,j,k)& - -tensorField_fourierMPI(l,2,1,j,k)*xi(3,1,j,k))*TWOPIIMG - curl_fourier = (+tensorField_fourierMPI(l,1,1,j,k)*xi(3,1,j,k)& - -tensorField_fourierMPI(l,3,1,j,k)*xi(1,1,j,k))*TWOPIIMG - curl_fourier = (+tensorField_fourierMPI(l,2,1,j,k)*xi(1,1,j,k)& - -tensorField_fourierMPI(l,1,1,j,k)*xi(2,1,j,k))*TWOPIIMG + curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi(2,1,j,k)& + -tensorField_fourier(l,2,1,j,k)*xi(3,1,j,k))*TWOPIIMG + curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi(3,1,j,k)& + -tensorField_fourier(l,3,1,j,k)*xi(1,1,j,k))*TWOPIIMG + curl_fourier = (+tensorField_fourier(l,2,1,j,k)*xi(1,1,j,k)& + -tensorField_fourier(l,1,1,j,k)*xi(2,1,j,k))*TWOPIIMG enddo utilities_curlRMS = utilities_curlRMS + & 2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) do l = 1_pInt, 3_pInt - curl_fourier = (+tensorField_fourierMPI(l,3,grid1Red,j,k)*xi(2,grid1Red,j,k)& - -tensorField_fourierMPI(l,2,grid1Red,j,k)*xi(3,grid1Red,j,k))*TWOPIIMG - curl_fourier = (+tensorField_fourierMPI(l,1,grid1Red,j,k)*xi(3,grid1Red,j,k)& - -tensorField_fourierMPI(l,3,grid1Red,j,k)*xi(1,grid1Red,j,k))*TWOPIIMG - curl_fourier = (+tensorField_fourierMPI(l,2,grid1Red,j,k)*xi(1,grid1Red,j,k)& - -tensorField_fourierMPI(l,1,grid1Red,j,k)*xi(2,grid1Red,j,k))*TWOPIIMG + curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi(2,grid1Red,j,k)& + -tensorField_fourier(l,2,grid1Red,j,k)*xi(3,grid1Red,j,k))*TWOPIIMG + curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi(3,grid1Red,j,k)& + -tensorField_fourier(l,3,grid1Red,j,k)*xi(1,grid1Red,j,k))*TWOPIIMG + curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi(1,grid1Red,j,k)& + -tensorField_fourier(l,1,grid1Red,j,k)*xi(2,grid1Red,j,k))*TWOPIIMG enddo utilities_curlRMS = utilities_curlRMS + & 2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) enddo; enddo call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) utilities_curlRMS = sqrt(utilities_curlRMS) * wgt - if(gridGlobal(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 end function utilities_curlRMS @@ -1048,15 +879,18 @@ subroutine utilities_fourierScalarGradient() use math, only: & PI use mesh, only: & - gridLocal, & - geomSizeGlobal + grid3, & + grid, & + geomSize + + implicit none integer(pInt) :: i, j, k - vectorField_fourierMPI = cmplx(0.0_pReal,0.0_pReal,pReal) - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,grid1Red - vectorField_fourierMPI(1:3,i,j,k) = scalarField_fourierMPI(i,j,k)* & + vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red + vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)* & cmplx(0.0_pReal,2.0_pReal*PI*xi(1:3,i,j,k)* & - scaledGeomSize/geomSizeGlobal,pReal) + scaledGeomSize/geomSize,pReal) enddo; enddo; enddo end subroutine utilities_fourierScalarGradient @@ -1068,17 +902,20 @@ subroutine utilities_fourierVectorDivergence() use math, only: & PI use mesh, only: & - gridLocal, & - geomSizeGlobal + grid3, & + grid, & + geomSize + + implicit none integer(pInt) :: i, j, k, m - scalarField_fourierMPI = cmplx(0.0_pReal,0.0_pReal,pReal) - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,grid1Red + scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red do m = 1_pInt, 3_pInt - scalarField_fourierMPI(i,j,k) = & - scalarField_fourierMPI(i,j,k) + & - vectorField_fourierMPI(m,i,j,k)* & - cmplx(0.0_pReal,2.0_pReal*PI*xi(m,i,j,k)*scaledGeomSize(m)/geomSizeGlobal(m),pReal) + scalarField_fourier(i,j,k) = & + scalarField_fourier(i,j,k) + & + vectorField_fourier(m,i,j,k)* & + cmplx(0.0_pReal,2.0_pReal*PI*xi(m,i,j,k)*scaledGeomSize(m)/geomSize(m),pReal) enddo enddo; enddo; enddo end subroutine utilities_fourierVectorDivergence @@ -1098,7 +935,8 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,& math_rotate_forward33, & math_det33 use mesh, only: & - gridLocal + grid,& + grid3 use FEsolving, only: & restartWrite use CPFEM, only: & @@ -1113,7 +951,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,& materialpoint_dPdF implicit none - real(pReal), intent(in), dimension(3,3,gridLocal(1),gridLocal(2),gridLocal(3)) :: & + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & F_lastInc, & !< target deformation gradient F !< previous deformation gradient real(pReal), intent(in) :: timeinc !< loading time @@ -1122,7 +960,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,& 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,gridLocal(1),gridLocal(2),gridLocal(3)) :: P !< PK stress + real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress integer(pInt) :: & calcMode, & !< CPFEM mode for calculation @@ -1130,8 +968,11 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,& real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet PetscErrorCode :: ierr + + external :: & + MPI_Allreduce - if (worldrank == 0_pInt) then + if (worldrank == 0_pInt) then write(6,'(/,a)') ' ... evaluating constitutive response ......................................' flush(6) endif @@ -1139,7 +980,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,& if (forwardData) then ! aging results calcMode = ior(calcMode, CPFEM_AGERESULTS) - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(gridLocal)]) + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) endif if (cutBack) then ! restore saved variables calcMode = iand(calcMode, not(CPFEM_AGERESULTS)) @@ -1148,7 +989,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,& call CPFEM_general(CPFEM_COLLECT,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), & timeinc,1_pInt,1_pInt) - materialpoint_F = reshape(F,[3,3,1,product(gridLocal)]) + materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) call debug_reset() !-------------------------------------------------------------------------------------------------- @@ -1156,7 +997,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,& if(debugGeneral) then defgradDetMax = -huge(1.0_pReal) defgradDetMin = +huge(1.0_pReal) - do j = 1_pInt, product(gridLocal) + do j = 1_pInt, product(grid(1:2))*grid3 defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j)) defgradDetMax = max(defgradDetMax,defgradDet) defgradDetMin = min(defgradDetMin,defgradDet) @@ -1177,7 +1018,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,& max_dPdF_norm = 0.0_pReal min_dPdF = huge(1.0_pReal) min_dPdF_norm = huge(1.0_pReal) - do k = 1_pInt, product(gridLocal) + do k = 1_pInt, product(grid(1:2))*grid3 if (max_dPdF_norm < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then max_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k) max_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal) @@ -1199,7 +1040,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,& restartWrite = .false. ! reset restartWrite status cutBack = .false. ! reset cutBack status - P = reshape(materialpoint_P, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)]) + P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) if (debugRotation .and. worldrank == 0_pInt) & @@ -1220,7 +1061,8 @@ end subroutine utilities_constitutiveResponse !-------------------------------------------------------------------------------------------------- pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,field) use mesh, only: & - gridLocal + grid3, & + grid implicit none real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon @@ -1228,17 +1070,16 @@ pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,fie timeinc_old !< timeinc of last step logical, intent(in) :: & guess !< guess along former trajectory - real(pReal), intent(in), dimension(3,3,gridLocal(1),gridLocal(2),gridLocal(3)) :: & + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & field_lastInc, & !< data of previous step field !< data of current step - real(pReal), dimension(3,3,gridLocal(1),gridLocal(2),gridLocal(3)) :: & + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & utilities_calculateRate if (guess) then utilities_calculateRate = (field-field_lastInc) / timeinc_old else - utilities_calculateRate = spread(spread(spread(avRate,3,gridLocal(1)),4,gridLocal(2)), & - 5,gridLocal(3)) + utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3) endif end function utilities_calculateRate @@ -1250,28 +1091,32 @@ end function utilities_calculateRate !-------------------------------------------------------------------------------------------------- function utilities_forwardField(timeinc,field_lastInc,rate,aim) use mesh, only: & - gridLocal - + grid3, & + grid + implicit none real(pReal), intent(in) :: & timeinc !< timeinc of current step - real(pReal), intent(in), dimension(3,3,gridLocal(1),gridLocal(2),gridLocal(3)) :: & + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & field_lastInc, & !< initial field rate !< rate by which to forward real(pReal), intent(in), optional, dimension(3,3) :: & aim !< average field value aim - real(pReal), dimension(3,3,gridLocal(1),gridLocal(2),gridLocal(3)) :: & + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & utilities_forwardField real(pReal), dimension(3,3) :: fieldDiff !< - aim PetscErrorCode :: ierr + external :: & + MPI_Allreduce + utilities_forwardField = field_lastInc + rate*timeinc 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,PETSC_COMM_WORLD,ierr) fieldDiff = fieldDiff - aim utilities_forwardField = utilities_forwardField - & - spread(spread(spread(fieldDiff,3,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) + spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid3) endif end function utilities_forwardField @@ -1280,44 +1125,41 @@ end function utilities_forwardField !-------------------------------------------------------------------------------------------------- !> @brief calculates filter for fourier convolution depending on type given in numerics.config !-------------------------------------------------------------------------------------------------- -real(pReal) function utilities_getFilter(k) - use IO, only: & - IO_error +complex(pReal) pure function utilities_getFilter(k) use numerics, only: & spectral_filter use math, only: & PI use mesh, only: & - gridGlobal - + grid + implicit none real(pReal),intent(in), dimension(3) :: k !< indices of frequency - utilities_getFilter = 1.0_pReal + utilities_getFilter = (1.0_pReal,0.0_pReal) select case (spectral_filter) case ('none') ! default, no weighting case ('cosine') ! cosine curve with 1 for avg and zero for highest freq - utilities_getFilter = product(1.0_pReal + cos(PI*k*scaledGeomSize/gridGlobal))/8.0_pReal + utilities_getFilter = cmplx(product(1.0_pReal + cos(PI*k*scaledGeomSize/grid))/8.0_pReal,& + 0.0_pReal) case ('gradient') ! gradient, might need grid scaling as for cosine filter - utilities_getFilter = 1.0_pReal/(1.0_pReal + & - (k(1)*k(1) + k(2)*k(2) + k(3)*k(3))) - case default - call IO_error(error_ID = 892_pInt, ext_msg = trim(spectral_filter)) + utilities_getFilter = cmplx(1.0_pReal/(1.0_pReal + (k(1)*k(1) + k(2)*k(2) + k(3)*k(3))),& + 0.0_pReal) end select - if (gridGlobal(1) /= 1_pInt .and. k(1) == real(grid1Red - 1_pInt, pReal)/scaledGeomSize(1)) & - utilities_getFilter = 0.0_pReal - if (gridGlobal(2) /= 1_pInt .and. k(2) == real(gridGlobal(2)/2_pInt, pReal)/scaledGeomSize(2)) & - utilities_getFilter = 0.0_pReal ! do not delete the whole slice in case of 2D calculation - if (gridGlobal(2) /= 1_pInt .and. & - k(2) == real(gridGlobal(2)/2_pInt + mod(gridGlobal(2),2_pInt), pReal)/scaledGeomSize(2)) & - utilities_getFilter = 0.0_pReal ! do not delete the whole slice in case of 2D calculation - if (gridGlobal(3) /= 1_pInt .and. k(3) == real(gridGlobal(3)/2_pInt, pReal)/scaledGeomSize(3)) & - utilities_getFilter = 0.0_pReal ! do not delete the whole slice in case of 2D calculation - if (gridGlobal(3) /= 1_pInt .and. & - k(3) == real(gridGlobal(3)/2_pInt + mod(gridGlobal(3),2_pInt), pReal)/scaledGeomSize(3)) & - utilities_getFilter = 0.0_pReal ! do not delete the whole slice in case of 2D calculation + if (grid(1) /= 1_pInt .and. k(1) == real(grid1Red - 1_pInt, pReal)/scaledGeomSize(1)) & + utilities_getFilter = (0.0_pReal,0.0_pReal) ! do not delete the whole slice in case of 2D calculation + if (grid(2) /= 1_pInt .and. k(2) == real(grid(2)/2_pInt, pReal)/scaledGeomSize(2)) & + utilities_getFilter = (0.0_pReal,0.0_pReal) ! do not delete the whole slice in case of 2D calculation + if (grid(2) /= 1_pInt .and. & + k(2) == real(grid(2)/2_pInt + mod(grid(2),2_pInt), pReal)/scaledGeomSize(2)) & + utilities_getFilter = (0.0_pReal,0.0_pReal) ! do not delete the whole slice in case of 2D calculation + if (grid(3) /= 1_pInt .and. k(3) == real(grid(3)/2_pInt, pReal)/scaledGeomSize(3)) & + utilities_getFilter = (0.0_pReal,0.0_pReal) ! do not delete the whole slice in case of 2D calculation + if (grid(3) /= 1_pInt .and. & + k(3) == real(grid(3)/2_pInt + mod(grid(3),2_pInt), pReal)/scaledGeomSize(3)) & + utilities_getFilter = (0.0_pReal,0.0_pReal) ! do not delete the whole slice in case of 2D calculation end function utilities_getFilter @@ -1330,9 +1172,6 @@ subroutine utilities_destroy() implicit none - if (debugDivergence) call fftw_destroy_plan(planDiv) - if (debugFFTW) call fftw_destroy_plan(planDebugForth) - if (debugFFTW) call fftw_destroy_plan(planDebugBack) call fftw_destroy_plan(planTensorForth) call fftw_destroy_plan(planTensorBack) call fftw_destroy_plan(planVectorForth) @@ -1349,58 +1188,62 @@ end subroutine utilities_destroy ! convolution !-------------------------------------------------------------------------------------------------- subroutine utilities_updateIPcoords(F) - use math + use math, only: & + PI, & + math_mul33x3 use mesh, only: & - gridGlobal, & - gridLocal, & - gridOffset, & - geomSizeGlobal, & + grid, & + grid3, & + grid3Offset, & + geomSize, & mesh_ipCoordinates implicit none - real(pReal), dimension(3,3,gridLocal(1),gridLocal(2),gridLocal(3)), intent(in) :: F + real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F integer(pInt) :: i, j, k, m real(pReal), dimension(3) :: step, offset_coords, integrator real(pReal), dimension(3,3) :: Favg PetscErrorCode :: ierr + external & + MPI_Bcast - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F call utilities_FFTtensorForward() - integrator = geomSizeGlobal * 0.5_pReal / PI - step = geomSizeGlobal/real(gridGlobal, pReal) + integrator = geomSize * 0.5_pReal / PI + step = geomSize/real(grid, pReal) !-------------------------------------------------------------------------------------------------- ! average F - if (gridOffset == 0_pInt) Favg = real(tensorField_fourierMPI(1:3,1:3,1,1,1),pReal)*wgt + if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! integration in Fourier space - vectorField_fourierMPI = cmplx(0.0_pReal, 0.0_pReal, pReal) - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,grid1Red + vectorField_fourier = cmplx(0.0_pReal, 0.0_pReal, pReal) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red do m = 1_pInt,3_pInt - vectorField_fourierMPI(m,i,j,k) = sum(tensorField_fourierMPI(m,1:3,i,j,k)*& + vectorField_fourier(m,i,j,k) = sum(tensorField_fourier(m,1:3,i,j,k)*& cmplx(0.0_pReal,xi(1:3,i,j,k)*scaledGeomSize*integrator,pReal)) enddo - if (any(xi(1:3,i,j,k) /= 0.0_pReal)) & - vectorField_fourierMPI(1:3,i,j,k) = & - vectorField_fourierMPI(1:3,i,j,k)/cmplx(-sum(xi(1:3,i,j,k)*scaledGeomSize*xi(1:3,i,j,k)* & + if (any(abs(xi(1:3,i,j,k)) > tiny(0.0_pReal))) & + vectorField_fourier(1:3,i,j,k) = & + vectorField_fourier(1:3,i,j,k)/cmplx(-sum(xi(1:3,i,j,k)*scaledGeomSize*xi(1:3,i,j,k)* & scaledGeomSize),0.0_pReal,pReal) enddo; enddo; enddo - call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourierMPI,vectorField_realMPI) + call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) !-------------------------------------------------------------------------------------------------- ! add average to fluctuation and put (0,0,0) on (0,0,0) - if (gridOffset == 0_pInt) offset_coords = vectorField_realMPI(1:3,1,1,1) + if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1) call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords m = 1_pInt - do k = 1_pInt,gridLocal(3); do j = 1_pInt,gridLocal(2); do i = 1_pInt,gridLocal(1) - mesh_ipCoordinates(1:3,1,m) = vectorField_realMPI(1:3,i,j,k) & + do k = 1_pInt,grid3; do j = 1_pInt,grid(2); do i = 1_pInt,grid(1) + mesh_ipCoordinates(1:3,1,m) = vectorField_real(1:3,i,j,k) & + offset_coords & - + math_mul33x3(Favg,step*real([i,j,k+gridOffset]-1_pInt,pReal)) + + math_mul33x3(Favg,step*real([i,j,k+grid3Offset]-1_pInt,pReal)) m = m+1_pInt enddo; enddo; enddo diff --git a/code/mesh.f90 b/code/mesh.f90 index b24501c5f..39b2cc57b 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -30,15 +30,17 @@ module mesh mesh_Nelems !< total number of elements in mesh #ifdef Spectral + integer(pInt), dimension(3), public, protected :: & + grid !< (global) grid integer(pInt), public, protected :: & mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh - gridGlobal(3), & - gridLocal (3), & - gridOffset + grid3, & !< (local) grid in 3rd direction + grid3Offset !< (local) grid offset in 3rd direction + real(pReal), dimension(3), public, protected :: & + geomSize real(pReal), public, protected :: & - geomSizeGlobal(3), & - geomSizeLocal (3), & - geomSizeOffset + size3, & !< (local) size in 3rd direction + size3offset !< (local) size offset in 3rd direction #endif integer(pInt), dimension(:,:), allocatable, public, protected :: & @@ -565,31 +567,21 @@ subroutine mesh_init(ip,el) #ifdef Spectral #ifdef PETSc call fftw_mpi_init() +#endif call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file... if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6) - - gridGlobal = mesh_spectral_getGrid(fileUnit) - gridMPI = int(gridGlobal,C_INTPTR_T) + grid = mesh_spectral_getGrid(fileUnit) + geomSize = mesh_spectral_getSize(fileUnit) + +#ifdef PETSc + gridMPI = int(grid,C_INTPTR_T) alloc_local = fftw_mpi_local_size_3d(gridMPI(3), gridMPI(2), gridMPI(1)/2 +1, & MPI_COMM_WORLD, local_K, local_K_offset) - gridLocal = [gridGlobal(1:2),int(local_K,pInt)] - gridOffset = int(local_K_offset,pInt) + grid3 = int(local_K,pInt) + grid3Offset = int(local_K_offset,pInt) - geomSizeGlobal = mesh_spectral_getSize(fileUnit) - geomSizeLocal = [geomSizeGlobal(1:2), & - geomSizeGlobal(3)*real(gridLocal(3),pReal)/real(gridGlobal(3),pReal)] - geomSizeOffset = geomSizeGlobal(3)*real(gridOffset,pReal) /real(gridGlobal(3),pReal) -#else - call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file... - if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6) - - gridGlobal = mesh_spectral_getGrid(fileUnit) - gridLocal = gridGlobal - gridOffset = 0_pInt - - geomSizeGlobal = mesh_spectral_getSize(fileUnit) - geomSizeLocal = geomSizeGlobal - geomSizeOffset = 0.0_pReal + size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) + size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal) #endif if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6) call mesh_spectral_count() @@ -1256,11 +1248,11 @@ subroutine mesh_spectral_count() implicit none - mesh_Nelems = product(gridLocal) + mesh_Nelems = product(grid(1:2))*grid3 mesh_NcpElems= mesh_Nelems - mesh_Nnodes = product(gridLocal + 1_pInt) + mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt) - mesh_NcpElemsGlobal = product(gridGlobal) + mesh_NcpElemsGlobal = product(grid) end subroutine mesh_spectral_count @@ -1319,15 +1311,15 @@ subroutine mesh_spectral_build_nodes() forall (n = 0_pInt:mesh_Nnodes-1_pInt) mesh_node0(1,n+1_pInt) = mesh_unitlength * & - geomSizeLocal(1)*real(mod(n,(gridLocal(1)+1_pInt) ),pReal) & - / real(gridLocal(1),pReal) + geomSize(1)*real(mod(n,(grid(1)+1_pInt) ),pReal) & + / real(grid(1),pReal) mesh_node0(2,n+1_pInt) = mesh_unitlength * & - geomSizeLocal(2)*real(mod(n/(gridLocal(1)+1_pInt),(gridLocal(2)+1_pInt)),pReal) & - / real(gridLocal(2),pReal) + geomSize(2)*real(mod(n/(grid(1)+1_pInt),(grid(2)+1_pInt)),pReal) & + / real(grid(2),pReal) mesh_node0(3,n+1_pInt) = mesh_unitlength * & - geomSizeLocal(3)*real(mod(n/(gridLocal(1)+1_pInt)/(gridLocal(2)+1_pInt),(gridLocal(3)+1_pInt)),pReal) & - / real(gridLocal(3),pReal) + & - geomSizeOffset + size3*real(mod(n/(grid(1)+1_pInt)/(grid(2)+1_pInt),(grid3+1_pInt)),pReal) & + / real(grid3,pReal) + & + size3offset end forall mesh_node = mesh_node0 @@ -1422,7 +1414,7 @@ subroutine mesh_spectral_build_elements(fileUnit) enddo elemType = FE_mapElemtype('C3D8R') - elemOffset = gridLocal(1)*gridLocal(2)*gridOffset + elemOffset = product(grid(1:2))*grid3Offset e = 0_pInt do while (e < mesh_NcpElems) ! fill expected number of elements, stop at end of data (or blank line!) e = e+1_pInt ! valid element entry @@ -1430,15 +1422,15 @@ subroutine mesh_spectral_build_elements(fileUnit) mesh_element( 2,e) = elemType ! elem type mesh_element( 3,e) = homog ! homogenization mesh_element( 4,e) = mesh_microGlobal(e+elemOffset) ! microstructure - mesh_element( 5,e) = e + (e-1_pInt)/gridLocal(1) + & - ((e-1_pInt)/(gridLocal(1)*gridLocal(2)))*(gridLocal(1)+1_pInt) ! base node + mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + & + ((e-1_pInt)/(grid(1)*grid(2)))*(grid(1)+1_pInt) ! base node mesh_element( 6,e) = mesh_element(5,e) + 1_pInt - mesh_element( 7,e) = mesh_element(5,e) + gridLocal(1) + 2_pInt - mesh_element( 8,e) = mesh_element(5,e) + gridLocal(1) + 1_pInt - mesh_element( 9,e) = mesh_element(5,e) +(gridLocal(1) + 1_pInt) * (gridLocal(2) + 1_pInt) ! second floor base node + mesh_element( 7,e) = mesh_element(5,e) + grid(1) + 2_pInt + mesh_element( 8,e) = mesh_element(5,e) + grid(1) + 1_pInt + mesh_element( 9,e) = mesh_element(5,e) +(grid(1) + 1_pInt) * (grid(2) + 1_pInt) ! second floor base node mesh_element(10,e) = mesh_element(9,e) + 1_pInt - mesh_element(11,e) = mesh_element(9,e) + gridLocal(1) + 2_pInt - mesh_element(12,e) = mesh_element(9,e) + gridLocal(1) + 1_pInt + mesh_element(11,e) = mesh_element(9,e) + grid(1) + 2_pInt + mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) ! needed for statistics mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) enddo @@ -1465,32 +1457,32 @@ subroutine mesh_spectral_build_ipNeighborhood(fileUnit) allocate(mesh_ipNeighborhood(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems),source=0_pInt) e = 0_pInt - do z = 0_pInt,gridLocal(3)-1_pInt - do y = 0_pInt,gridLocal(2)-1_pInt - do x = 0_pInt,gridLocal(1)-1_pInt + do z = 0_pInt,grid3-1_pInt + do y = 0_pInt,grid(2)-1_pInt + do x = 0_pInt,grid(1)-1_pInt e = e + 1_pInt - mesh_ipNeighborhood(1,1,1,e) = z * gridLocal(1) * gridLocal(2) & - + y * gridLocal(1) & - + modulo(x+1_pInt,gridLocal(1)) & + mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x+1_pInt,grid(1)) & + 1_pInt - mesh_ipNeighborhood(1,2,1,e) = z * gridLocal(1) * gridLocal(2) & - + y * gridLocal(1) & - + modulo(x-1_pInt,gridLocal(1)) & + mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x-1_pInt,grid(1)) & + 1_pInt - mesh_ipNeighborhood(1,3,1,e) = z * gridLocal(1) * gridLocal(2) & - + modulo(y+1_pInt,gridLocal(2)) * gridLocal(1) & + mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) & + + modulo(y+1_pInt,grid(2)) * grid(1) & + x & + 1_pInt - mesh_ipNeighborhood(1,4,1,e) = z * gridLocal(1) * gridLocal(2) & - + modulo(y-1_pInt,gridLocal(2)) * gridLocal(1) & + mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) & + + modulo(y-1_pInt,grid(2)) * grid(1) & + x & + 1_pInt - mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,gridLocal(3)) * gridLocal(1) * gridLocal(2) & - + y * gridLocal(1) & + mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,grid3) * grid(1) * grid(2) & + + y * grid(1) & + x & + 1_pInt - mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,gridLocal(3)) * gridLocal(1) * gridLocal(2) & - + y * gridLocal(1) & + mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,grid3) * grid(1) * grid(2) & + + y * grid(1) & + x & + 1_pInt mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt diff --git a/code/spectral_damage.f90 b/code/spectral_damage.f90 index 6e90e6aaf..a358250ef 100644 --- a/code/spectral_damage.f90 +++ b/code/spectral_damage.f90 @@ -85,8 +85,8 @@ subroutine spectral_damage_init() use DAMASK_spectral_Utilities, only: & wgt use mesh, only: & - gridLocal, & - gridGlobal + grid, & + grid3 use damage_nonlocal, only: & damage_nonlocal_getDiffusion33, & damage_nonlocal_getMobility @@ -112,17 +112,17 @@ subroutine spectral_damage_init() ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 do proc = 1, worldsize call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) enddo call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & !< cut off stencil at boundary DMDA_STENCIL_BOX, & !< Moore (26) neighborhood around central point - gridGlobal(1),gridGlobal(2),gridGlobal(3), & !< global grid + grid(1),grid(2),grid(3), & !< global grid 1, 1, worldsize, & 1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap) - gridLocal (1),gridLocal (2),localK, & !< local grid + grid(1),grid(2),localK, & !< local grid damage_grid,ierr) !< handle, error CHKERRQ(ierr) call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da @@ -150,16 +150,16 @@ subroutine spectral_damage_init() yend = ystart + yend - 1 zend = zstart + zend - 1 call VecSet(solution,1.0,ierr); CHKERRQ(ierr) - allocate(damage_current(gridLocal(1),gridLocal(2),gridLocal(3)), source=1.0_pReal) - allocate(damage_lastInc(gridLocal(1),gridLocal(2),gridLocal(3)), source=1.0_pReal) - allocate(damage_stagInc(gridLocal(1),gridLocal(2),gridLocal(3)), source=1.0_pReal) + 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) !-------------------------------------------------------------------------------------------------- ! damage reference diffusion update cell = 0_pInt D_ref = 0.0_pReal mobility_ref = 0.0_pReal - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) @@ -184,7 +184,8 @@ type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old Utilities_maskedCompliance, & Utilities_updateGamma use mesh, only: & - gridLocal + grid, & + grid3 use damage_nonlocal, only: & damage_nonlocal_putNonLocalDamage @@ -234,7 +235,7 @@ type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old !-------------------------------------------------------------------------------------------------- ! updating damage state cell = 0_pInt !< material point = 0 - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt !< material point increase call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) enddo; enddo; enddo @@ -260,12 +261,13 @@ subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr) use numerics, only: & residualStiffness use mesh, only: & - gridLocal + grid, & + grid3 use math, only: & math_mul33x3 use DAMASK_spectral_Utilities, only: & - scalarField_realMPI, & - vectorField_realMPI, & + scalarField_real, & + vectorField_real, & utilities_FFTvectorForward, & utilities_FFTvectorBackward, & utilities_FFTscalarForward, & @@ -295,30 +297,30 @@ subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr) damage_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_realMPI = 0.0_pReal - scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = damage_current + scalarField_real = 0.0_pReal + scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current call utilities_FFTscalarForward() call utilities_fourierScalarGradient() !< calculate gradient of damage field call utilities_FFTvectorBackward() cell = 0_pInt - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt - vectorField_realMPI(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, & - vectorField_realMPI(1:3,i,j,k)) + vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_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_FFTscalarBackward() cell = 0_pInt - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell) mobility = damage_nonlocal_getMobility(1,cell) - scalarField_realMPI(i,j,k) = params%timeinc*scalarField_realMPI(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) + & + params%timeinc*phiDot + & + mobility*damage_lastInc(i,j,k) - & + mobility*damage_current(i,j,k) + & + mobility_ref*damage_current(i,j,k) enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- @@ -326,13 +328,12 @@ subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr) call utilities_FFTscalarForward() call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) call utilities_FFTscalarBackward() - where(scalarField_realMPI > damage_lastInc) scalarField_realMPI = damage_lastInc - where(scalarField_realMPI < residualStiffness) scalarField_realMPI = residualStiffness + where(scalarField_real > damage_lastInc) scalarField_real = damage_lastInc + where(scalarField_real < residualStiffness) scalarField_real = residualStiffness !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) - & - damage_current + f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current end subroutine spectral_damage_formResidual @@ -341,7 +342,8 @@ end subroutine spectral_damage_formResidual !-------------------------------------------------------------------------------------------------- subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime) use mesh, only: & - gridLocal + grid, & + grid3 use DAMASK_spectral_Utilities, only: & cutBack, & wgt @@ -371,7 +373,7 @@ subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime) call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) enddo; enddo; enddo @@ -382,7 +384,7 @@ subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime) cell = 0_pInt D_ref = 0.0_pReal mobility_ref = 0.0_pReal - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) diff --git a/code/spectral_thermal.f90 b/code/spectral_thermal.f90 index 59df82baa..831df5587 100644 --- a/code/spectral_thermal.f90 +++ b/code/spectral_thermal.f90 @@ -85,8 +85,8 @@ subroutine spectral_thermal_init use DAMASK_spectral_Utilities, only: & wgt use mesh, only: & - gridLocal, & - gridGlobal + grid, & + grid3 use thermal_conduction, only: & thermal_conduction_getConductivity33, & thermal_conduction_getMassDensity, & @@ -116,17 +116,17 @@ subroutine spectral_thermal_init ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 do proc = 1, worldsize call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) enddo call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid + grid(1),grid(2),grid(3), & ! global grid 1, 1, worldsize, & 1, 0, & ! #dof (temperature field), ghost boundary width (domain overlap) - gridLocal (1),gridLocal (2),localK, & ! local grid + grid (1),grid(2),localK, & ! local grid thermal_grid,ierr) ! handle, error CHKERRQ(ierr) call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da @@ -142,11 +142,11 @@ subroutine spectral_thermal_init xend = xstart + xend - 1 yend = ystart + yend - 1 zend = zstart + zend - 1 - allocate(temperature_current(gridLocal(1),gridLocal(2),gridLocal(3)), source=0.0_pReal) - allocate(temperature_lastInc(gridLocal(1),gridLocal(2),gridLocal(3)), source=0.0_pReal) - allocate(temperature_stagInc(gridLocal(1),gridLocal(2),gridLocal(3)), source=0.0_pReal) + 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) cell = 0_pInt - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt temperature_current(i,j,k) = temperature(mappingHomogenization(2,1,cell))% & p(thermalMapping(mappingHomogenization(2,1,cell))%p(1,cell)) @@ -160,7 +160,7 @@ subroutine spectral_thermal_init cell = 0_pInt D_ref = 0.0_pReal mobility_ref = 0.0_pReal - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* & @@ -186,7 +186,8 @@ type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_ol Utilities_maskedCompliance, & Utilities_updateGamma use mesh, only: & - gridLocal + grid, & + grid3 use thermal_conduction, only: & thermal_conduction_putTemperatureAndItsRate @@ -236,7 +237,7 @@ type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_ol !-------------------------------------------------------------------------------------------------- ! updating thermal state cell = 0_pInt !< material point = 0 - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt !< material point increase call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & (temperature_current(i,j,k)-temperature_lastInc(i,j,k))/params%timeinc, & @@ -262,12 +263,13 @@ end function spectral_thermal_solution !-------------------------------------------------------------------------------------------------- subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr) use mesh, only: & - gridLocal + grid, & + grid3 use math, only: & math_mul33x3 use DAMASK_spectral_Utilities, only: & - scalarField_realMPI, & - vectorField_realMPI, & + scalarField_real, & + vectorField_real, & utilities_FFTvectorForward, & utilities_FFTvectorBackward, & utilities_FFTscalarForward, & @@ -298,30 +300,30 @@ subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr) temperature_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_realMPI = 0.0_pReal - scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = temperature_current + scalarField_real = 0.0_pReal + scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current call utilities_FFTscalarForward() call utilities_fourierScalarGradient() !< calculate gradient of damage field call utilities_FFTvectorBackward() cell = 0_pInt - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt - vectorField_realMPI(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & - vectorField_realMPI(1:3,i,j,k)) + vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_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_FFTscalarBackward() cell = 0_pInt - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell) - scalarField_realMPI(i,j,k) = params%timeinc*scalarField_realMPI(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) + 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) enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- @@ -332,7 +334,7 @@ subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = temperature_current - scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) + f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) end subroutine spectral_thermal_formResidual @@ -341,7 +343,8 @@ end subroutine spectral_thermal_formResidual !-------------------------------------------------------------------------------------------------- subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime) use mesh, only: & - gridLocal + grid, & + grid3 use DAMASK_spectral_Utilities, only: & cutBack, & wgt @@ -373,7 +376,7 @@ subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime) call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt !< material point increase call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & (temperature_current(i,j,k) - & @@ -387,7 +390,7 @@ subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime) cell = 0_pInt D_ref = 0.0_pReal mobility_ref = 0.0_pReal - do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) cell = cell + 1_pInt D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* &