diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index 7d739b35d..c325431d1 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -51,12 +51,13 @@ program DAMASK_spectral_Driver IO_read_jobBinaryFile, & IO_write_jobBinaryFile, & IO_intOut, & - IO_warning + IO_warning, & + IO_timeStamp + use debug, only: & + debug_level, & + debug_spectral, & + debug_levelBasic use math ! need to include the whole module for FFTW - use mesh, only : & - res, & - geomdim, & - mesh_NcpElems use CPFEM, only: & CPFEM_initAll use FEsolving, only: & @@ -71,10 +72,10 @@ program DAMASK_spectral_Driver materialpoint_sizeResults, & materialpoint_results use DAMASK_spectral_Utilities, only: & + grid, & + geomSize, & tBoundaryCondition, & tSolutionState, & - debugGeneral, & - debugDivergence, & cutBack use DAMASK_spectral_SolverBasic #ifdef PETSc @@ -150,8 +151,9 @@ program DAMASK_spectral_Driver !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) call CPFEM_initAll(temperature = 300.0_pReal, element = 1_pInt, IP= 1_pInt) - write(6,'(/,a)') ' <<<+- DAMASK_spectral_driver init -+>>>' - write(6,'(a)') ' $Id$' + write(6,'(/,a)') ' <<<+- DAMASK_spectral_driver init -+>>>' + write(6,'(a)') ' $Id$' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" !-------------------------------------------------------------------------------------------------- @@ -327,7 +329,8 @@ program DAMASK_spectral_Driver case (DAMASK_spectral_SolverBasicPETSc_label) call basicPETSc_init(loadCases(1)%temperature) case (DAMASK_spectral_SolverAL_label) - if(debugDivergence) call IO_warning(42_pInt, ext_msg='debug Divergence') + if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & + call IO_warning(42_pInt, ext_msg='debug Divergence') call AL_init(loadCases(1)%temperature) #endif case default @@ -349,10 +352,10 @@ 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', res - !write(resUnit) 'size', geomdim - write(resUnit) 'resolution', res - write(resUnit) 'dimension', geomdim + !write(resUnit) 'grid', grid + !write(resUnit) 'size', geomSize + write(resUnit) 'resolution', grid + write(resUnit) 'dimension', geomSize write(resUnit) 'materialpoint_sizeResults', materialpoint_sizeResults write(resUnit) 'loadcases', size(loadCases) write(resUnit) 'frequencies', loadCases%outputfrequency ! one entry per currentLoadCase @@ -361,11 +364,12 @@ program DAMASK_spectral_Driver write(resUnit) 'increments', loadCases%incs ! one entry per currentLoadCase write(resUnit) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc write(resUnit) 'eoh' ! end of header - write(resUnit) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! initial (non-deformed or read-in) results + write(resUnit) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:product(grid)) ! initial (non-deformed or read-in) results open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& '.sta',form='FORMATTED',status='REPLACE') write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file - if (debugGeneral) write(6,'(/,a)') ' header of result file written out' + if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) & + write(6,'(/,a)') ' header of result file written out' flush(6) endif diff --git a/code/DAMASK_spectral_interface.f90 b/code/DAMASK_spectral_interface.f90 index b6e676822..df43ad471 100644 --- a/code/DAMASK_spectral_interface.f90 +++ b/code/DAMASK_spectral_interface.f90 @@ -250,7 +250,7 @@ character(len=1024) function storeWorkingDirectory(workingDirectoryArg,geometryA endif if (storeWorkingDirectory(len(trim(storeWorkingDirectory)):len(trim(storeWorkingDirectory))) & ! if path seperator is not given, append it /= pathSep) storeWorkingDirectory = trim(storeWorkingDirectory)//pathSep - !here check if exists and use chdir! + !> @ToDO here check if exists and use chdir! else ! using path to geometry file as working dir if (geometryArg(1:1) == pathSep) then ! absolute path given as command line argument storeWorkingDirectory = geometryArg(1:scan(geometryArg,pathSep,back=.true.)) @@ -383,8 +383,8 @@ function rectifyPath(path) l = len_trim(path) rectifyPath = path do i = l,3,-1 - if ( rectifyPath(i-1:i) == '.'//pathSep .and. rectifyPath(i-2:i-2) /= '.' ) & - rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' + if ( rectifyPath(i-1:i) == '.'//pathSep .and. rectifyPath(i-2:i-2) /= '.' ) & + rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' enddo !remove ../ and corresponding directory from rectifyPath diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index aad6b9632..9b27de989 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -94,22 +94,22 @@ module DAMASK_spectral_solverAL AL_solution, & AL_destroy external :: & - VecDestroy, & - DMDestroy, & - DMDACreate3D, & - DMCreateGlobalVector, & - DMDASetLocalFunction, & - PETScFinalize, & - SNESDestroy, & - SNESGetNumberFunctionEvals, & - SNESGetIterationNumber, & - SNESSolve, & - SNESSetDM, & - SNESGetConvergedReason, & - SNESSetConvergenceTest, & - SNESSetFromOptions, & - SNESCreate, & - MPI_Abort + VecDestroy, & + DMDestroy, & + DMDACreate3D, & + DMCreateGlobalVector, & + DMDASetLocalFunction, & + PETScFinalize, & + SNESDestroy, & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber, & + SNESSolve, & + SNESSetDM, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions, & + SNESCreate, & + MPI_Abort contains @@ -119,9 +119,14 @@ contains subroutine AL_init(temperature) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use IO, only: & + IO_intOut, & IO_read_JobBinaryFile, & IO_write_JobBinaryFile, & IO_timeStamp + use debug, only : & + debug_level, & + debug_spectral, & + debug_spectralRestart use FEsolving, only: & restartInc use DAMASK_interface, only: & @@ -130,14 +135,12 @@ subroutine AL_init(temperature) Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & - debugRestart + grid, & + geomSize, & + wgt use numerics, only: & petsc_options use mesh, only: & - res, & - geomdim, & - wgt, & - mesh_NcpElems, & mesh_ipCoordinates, & mesh_deformedCoordsFFT use math, only: & @@ -148,7 +151,7 @@ subroutine AL_init(temperature) temperature #include #include - real(pReal), dimension(3,3, res(1), res(2),res(3)) :: P + real(pReal), dimension(:,:,:,:,:), allocatable :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal real(pReal), dimension(3,3,3,3) :: & @@ -164,18 +167,21 @@ subroutine AL_init(temperature) write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' write(6,'(a16,a)') ' Current time : ',IO_timeStamp() #include "compilation_info.f90" - - allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) - allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) !< @Todo sourced allocation allocate(Fdot,source = F_lastInc) - allocate (F_tau_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal) - allocate (F_tauDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal) + + allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) +!-------------------------------------------------------------------------------------------------- +! allocate global fields + allocate (F_lastInc (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) !< @Todo sourced allocation allocate(Fdot,source = F_lastInc) + allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (F_tau_lastInc(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (F_tauDot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! PETSc Init call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call DMDACreate3d(PETSC_COMM_WORLD, & DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & - DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & + DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & 18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) @@ -191,13 +197,14 @@ subroutine AL_init(temperature) 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,res(1)),4,res(2)),5,res(3)) ! initialize to identity + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity F_tau_lastInc = F_lastInc - F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) F_tau = F - elseif (restartInc > 1_pInt) then ! using old values from file - if (debugRestart) write(6,'(/,a,i6,a)') ' reading values of increment ',& - restartInc - 1_pInt,' from file' + elseif (restartInc > 1_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading values of increment', restartInc - 1_pInt, 'from file' flush(6) call IO_read_jobBinaryFile(777,'F',& trim(getSolverJobName()),size(F)) @@ -230,8 +237,8 @@ subroutine AL_init(temperature) read (777,rec=1) C_minmaxAvg close (777) endif - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(& - F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real,temp3333_Real2,& temp33_Real,.false.,math_I3) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) @@ -263,14 +270,13 @@ type(tSolutionState) function & math_rotate_backward33, & math_invSym3333 use mesh, only: & - res, & - geomdim, & - mesh_NcpElems, & mesh_ipCoordinates, & mesh_deformedCoordsFFT use IO, only: & IO_write_JobBinaryFile use DAMASK_spectral_Utilities, only: & + grid, & + geomSize, & tBoundaryCondition, & Utilities_forwardField, & Utilities_calculateRate, & @@ -342,8 +348,8 @@ use mesh, only: & if ( cutBack) then F_aim = F_aim_lastInc - F_tau= reshape(F_tau_lastInc,[9,res(1),res(2),res(3)]) - F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid(3)]) + F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)]) C = C_lastInc else C_lastInc = C @@ -359,23 +365,23 @@ use mesh, only: & !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(& - F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) + timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) F_tauDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,res(1),res(2),res(3)])) + timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid(3)])) - F_lastInc = reshape(F, [3,3,res(1),res(2),res(3)]) - F_tau_lastInc = reshape(F_tau,[3,3,res(1),res(2),res(3)]) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) + F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid(3)]) 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,res(1),res(2),res(3)]) - F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), [9,res(1),res(2),res(3)]) ! does not have any average value as boundary condition + math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)]) + F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) CHKERRQ(ierr) @@ -422,22 +428,25 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) itmin, & polarAlpha, & polarBeta + use IO, only: & + IO_intOut use math, only: & math_rotate_backward33, & math_transpose33, & math_mul3333xx33, & math_invSym3333 - use mesh, only: & - res, & - wgt use DAMASK_spectral_Utilities, only: & + grid, & + wgt, & field_real, & Utilities_FFTforward, & Utilities_fourierConvolution, & Utilities_FFTbackward, & - Utilities_constitutiveResponse, & - debugRotation - use IO, only: IO_intOut + Utilities_constitutiveResponse + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation use homogenization, only: & materialpoint_P, & materialpoint_dPdF @@ -490,7 +499,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) if (callNo == 0 .or. mod(callNo,2) == 1_pInt) then write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & ' @ Iteration ', itmin, '≤',reportIter, '≤', itmax - if (debugRotation) & + 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 =', & @@ -503,7 +512,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! field_real = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_scale,(polarAlpha + polarBeta)*F(1:3,1:3,i,j,k) - & (polarAlpha)*F_tau(1:3,1:3,i,j,k)) enddo; enddo; enddo @@ -516,8 +525,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_tau = polarBeta*F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),& - [3,3,res(1),res(2),res(3)],order=[3,4,5,1,2]) + residual_F_tau = polarBeta*F - reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),& + [3,3,grid(1),grid(2),grid(3)],order=[3,4,5,1,2]) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -534,7 +543,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) ! constructing residual n_ele = 0_pInt err_p = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) n_ele = n_ele + 1_pInt err_p = err_p + sum((math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) - & (F_tau(1:3,1:3,i,j,k) - & diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index 44ec7f164..efa98fb0e 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -73,6 +73,10 @@ subroutine basic_init(temperature) IO_write_JobBinaryFile, & IO_intOut, & IO_timeStamp + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRestart use FEsolving, only: & restartInc use DAMASK_interface, only: & @@ -81,20 +85,17 @@ subroutine basic_init(temperature) Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & - debugRestart - use mesh, only: & - res, & + grid, & wgt, & - geomdim, & - scaledDim, & + geomSize + use mesh, only: & mesh_ipCoordinates, & - mesh_NcpElems, & mesh_deformedCoordsFFT implicit none real(pReal), intent(inout) :: & temperature - real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P + real(pReal), dimension(:,:,:,:,:), allocatable :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal real(pReal), dimension(3,3,3,3) :: & @@ -105,22 +106,23 @@ subroutine basic_init(temperature) write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - write(6,'(a,3(f12.5)/)') ' scaledDim x y z:', scaledDim + allocate (P (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F (3,3,res(1), res(2),res(3)), source = 0.0_pReal) - allocate (F_lastInc (3,3,res(1), res(2),res(3)), source = 0.0_pReal) - allocate (Fdot (3,3,res(1), res(2),res(3)), source = 0.0_pReal) + allocate (F (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal) + allocate (F_lastInc (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal) + allocate (Fdot (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! init fields and average quantities if (restartInc == 1_pInt) then ! no deformation (no restart) - F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity + F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity F_lastInc = F elseif (restartInc > 1_pInt) then ! using old values from file - if (debugRestart) write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading values of increment', restartInc - 1_pInt, 'from file' + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading values of increment', restartInc - 1_pInt, 'from file' flush(6) call IO_read_jobBinaryFile(777,'F',& trim(getSolverJobName()),size(F)) @@ -147,8 +149,9 @@ subroutine basic_init(temperature) read (777,rec=1) temp3333_Real close (777) endif - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,F),[3,1,mesh_NcpElems]) - call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,F),[3,1,product(grid)]) + call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,C_minmaxAvg,& + temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness temp3333_Real = C_minmaxAvg endif @@ -173,15 +176,15 @@ type(tSolutionState) function & math_transpose33, & math_mul3333xx33 use mesh, only: & - res,& - geomdim, & - wgt, & mesh_ipCoordinates,& - mesh_NcpElems, & mesh_deformedCoordsFFT use IO, only: & IO_write_JobBinaryFile, & IO_intOut + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation use DAMASK_spectral_Utilities, only: & tBoundaryCondition, & field_real, & @@ -194,7 +197,9 @@ type(tSolutionState) function & Utilities_updateGamma, & Utilities_constitutiveResponse, & Utilities_calculateRate, & - debugRotation + grid,& + geomSize, & + wgt use FEsolving, only: & restartWrite, & restartRead, & @@ -224,7 +229,7 @@ type(tSolutionState) function & real(pReal), dimension(3,3) :: & F_aim_lastIter, & !< aim of last iteration P_av - real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P + real(pReal), dimension(3,3,grid(1),grid(2),grid(3)) :: P !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. real(pReal) :: err_div, err_stress @@ -261,7 +266,7 @@ type(tSolutionState) function & C = C_lastInc else C_lastInc = C - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,F),[3,1,mesh_NcpElems]) + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,F),[3,1,product(grid)]) !-------------------------------------------------------------------------------------------------- ! calculate rate for aim @@ -298,7 +303,7 @@ type(tSolutionState) function & ! report begin of new iteration write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & ' @ Iteration ', itmin, '≤',iter, '≤', itmax - if (debugRotation) & + 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,rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & @@ -322,13 +327,13 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme field_real = 0.0_pReal - field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(P,[res(1),res(2),res(3),3,3],& - order=[4,5,1,2,3]) ! field real has a different order + field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(P,[grid(1),grid(2),grid(3),3,3],& + order=[4,5,1,2,3]) ! field real has a different order call Utilities_FFTforward() err_div = Utilities_divergenceRMS() call Utilities_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,rotation_BC)) call Utilities_FFTbackward() - F = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization + F = F - reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av) write(6,'(/,a)') ' ==========================================================================' flush(6) @@ -368,7 +373,7 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress) err_stress_tol, & pAvgDivL2 - pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(pAvgDiv,math_transpose33(pAvgDiv))))) ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html) + pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(pAvgDiv,math_transpose33(pAvgDiv))))) ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html) err_stress_tol = max(maxval(abs(pAvgStress))*err_stress_tolrel,err_stress_tolabs) basic_Converged = all([ err_div/pAvgDivL2/err_div_tol,& diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 3ba457ce2..70d9000f8 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -84,22 +84,22 @@ module DAMASK_spectral_SolverBasicPETSc basicPETSc_solution ,& basicPETSc_destroy external :: & - VecDestroy, & - DMDestroy, & - DMDACreate3D, & - DMCreateGlobalVector, & - DMDASetLocalFunction, & - PETScFinalize, & - SNESDestroy, & - SNESGetNumberFunctionEvals, & - SNESGetIterationNumber, & - SNESSolve, & - SNESSetDM, & - SNESGetConvergedReason, & - SNESSetConvergenceTest, & - SNESSetFromOptions, & - SNESCreate, & - MPI_Abort + VecDestroy, & + DMDestroy, & + DMDACreate3D, & + DMCreateGlobalVector, & + DMDASetLocalFunction, & + PETScFinalize, & + SNESDestroy, & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber, & + SNESSolve, & + SNESSetDM, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions, & + SNESCreate, & + MPI_Abort contains @@ -111,7 +111,12 @@ subroutine basicPETSc_init(temperature) use IO, only: & IO_read_JobBinaryFile, & IO_write_JobBinaryFile, & + IO_intOut, & IO_timeStamp + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRestart use FEsolving, only: & restartInc use DAMASK_interface, only: & @@ -120,15 +125,12 @@ subroutine basicPETSc_init(temperature) Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & - debugRestart + grid, & + wgt, & + geomSize use numerics, only: & petsc_options use mesh, only: & - res, & - wgt, & - geomdim, & - scaledDim, & - mesh_NcpElems, & mesh_ipCoordinates, & mesh_deformedCoordsFFT use math, only: & @@ -139,7 +141,7 @@ subroutine basicPETSc_init(temperature) temperature #include #include - real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P + real(pReal), dimension(:,:,:,:,:), allocatable :: P PetscScalar, dimension(:,:,:,:), pointer :: F PetscErrorCode :: ierr PetscObject :: dummy @@ -153,19 +155,19 @@ subroutine basicPETSc_init(temperature) write(6,'(a)') ' $Id: DAMASK_spectral_SolverBasicPETSC.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - write(6,'(a,3(f12.5)/)') ' scaledDim x y z:', scaledDim + allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc(3,3,res(1),res(2),res(3)), source = 0.0_pReal) - allocate (Fdot (3,3,res(1),res(2),res(3)), source = 0.0_pReal) + allocate (F_lastInc(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call DMDACreate3d(PETSC_COMM_WORLD, & DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & - DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & + DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & 9,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) @@ -180,11 +182,12 @@ subroutine basicPETSc_init(temperature) 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,res(1)),4,res(2)),5,res(3)) ! initialize to identity - F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) elseif (restartInc > 1_pInt) then ! using old values from file - if (debugRestart) write(6,'(/,a,i6,a)') ' reading values of increment ',& - restartInc - 1_pInt,' from file' + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading values of increment', restartInc - 1_pInt, 'from file' flush(6) call IO_read_jobBinaryFile(777,'F',& trim(getSolverJobName()),size(F)) @@ -210,12 +213,12 @@ subroutine basicPETSc_init(temperature) read (777,rec=1) temp3333_Real close (777) endif - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(& - F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) call Utilities_constitutiveResponse(& - reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& - reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& - temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3) + reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),& + reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),& + temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back into PETSc if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness temp3333_Real = C_minmaxAvg @@ -237,14 +240,13 @@ type(tSolutionState) function & math_mul33x33 ,& math_rotate_backward33 use mesh, only: & - res,& - geomdim,& mesh_ipCoordinates,& - mesh_NcpElems, & mesh_deformedCoordsFFT use IO, only: & IO_write_JobBinaryFile use DAMASK_spectral_Utilities, only: & + grid, & + geomSize, & tBoundaryCondition, & Utilities_calculateRate, & Utilities_forwardField, & @@ -295,16 +297,16 @@ type(tSolutionState) function & write (777,rec=1) C_lastInc close(777) endif - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(& - F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) if ( cutBack) then F_aim = F_aim_lastInc - F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) C = C_lastInc else C_lastInc = C - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(& - F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) !-------------------------------------------------------------------------------------------------- ! calculate rate for aim @@ -319,13 +321,13 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update rate and forward last inc Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,params%rotation_BC), & - timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) - F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) + timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) + F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)]) endif F_aim = F_aim + f_aimDot * timeinc F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, & - rotation_BC)),[9,res(1),res(2),res(3)]) + rotation_BC)),[9,grid(1),grid(2),grid(3)]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -370,23 +372,26 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) math_rotate_backward33, & math_transpose33, & math_mul3333xx33 - use mesh, only: & - res, & - wgt + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation use DAMASK_spectral_Utilities, only: & + grid, & + wgt, & field_real, & Utilities_FFTforward, & Utilities_FFTbackward, & Utilities_fourierConvolution, & Utilities_constitutiveResponse, & - Utilities_divergenceRMS, & - debugRotation - use IO, only : IO_intOut + Utilities_divergenceRMS + use IO, only: & + IO_intOut implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in - PetscScalar, dimension(3,3,res(1),res(2),res(3)) :: & + PetscScalar, dimension(3,3,grid(1),grid(2),grid(3)) :: & x_scal, & f_scal PetscInt :: & @@ -408,7 +413,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) if (callNo == 0 .or. mod(callNo,2) == 1_pInt) then write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & ' @ Iteration ', itmin, '≤',reportIter, '≤', itmax - if (debugRotation) & + 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 =', & @@ -433,7 +438,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme field_real = 0.0_pReal - field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(f_scal,[res(1),res(2),res(3),3,3],& + field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(f_scal,[grid(1),grid(2),grid(3),3,3],& order=[4,5,1,2,3]) ! field real has a different order call Utilities_FFTforward() err_div = Utilities_divergenceRMS() @@ -442,7 +447,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2]) + f_scal = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2]) end subroutine BasicPETSc_formResidual diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 3f980c51f..ad9cb86c0 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -36,14 +36,21 @@ module DAMASK_spectral_utilities #include #endif logical, public :: cutBack =.false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill - +!-------------------------------------------------------------------------------------------------- +! grid related information information + integer(pInt), public, dimension(3) :: grid !< grid points as specified in geometry file + real(pReal), public :: wgt !< weighting factor 1/Nelems + real(pReal), public, dimension(3) :: geomSize !< size of geometry as specified in geometry file + !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW + integer(pInt), public :: grid1Red !< grid(1)/2 real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress or deformation) of field_fourier complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_fourier !< field on which the Fourier transform operates 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 !< reference stiffness + real(pReal), private, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence (Basic, Basic PETSc) !-------------------------------------------------------------------------------------------------- ! debug fftw @@ -66,10 +73,9 @@ module DAMASK_spectral_utilities !-------------------------------------------------------------------------------------------------- ! variables controlling debugging - logical, public :: & + logical, private :: & debugGeneral, & !< general debugging of spectral solver debugDivergence, & !< debugging of divergence calculation (comparison to function used for post processing) - debugRestart, & !< debbuging of restart features 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 @@ -118,22 +124,25 @@ contains !-------------------------------------------------------------------------------------------------- subroutine utilities_init() use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use DAMASK_interface, only: & + geometryFile use IO, only: & IO_error, & IO_warning, & - IO_timeStamp + IO_timeStamp, & + IO_open_file use numerics, only: & DAMASK_NumThreadsInt, & fftw_planner_flag, & fftw_timelimit, & memory_efficient, & - petsc_options + petsc_options, & + divergence_correction use debug, only: & debug_level, & debug_spectral, & debug_levelBasic, & debug_spectralDivergence, & - debug_spectralRestart, & debug_spectralFFTW, & debug_spectralPETSc, & debug_spectralRotation @@ -142,11 +151,10 @@ subroutine utilities_init() debug_spectralPETSc, & PETScDebug #endif - use mesh, only: & - res, & - res1_red, & - scaledDim use math ! must use the whole module for use of FFTW + use mesh, only: & + mesh_spectral_getSize, & + mesh_spectral_getGrid implicit none #ifdef PETSc @@ -157,6 +165,7 @@ subroutine utilities_init() PetscErrorCode :: ierr #endif integer(pInt) :: i, j, k + integer(pInt), parameter :: fileUnit = 228_pInt integer(pInt), dimension(3) :: k_s type(C_PTR) :: & tensorField, & !< field cotaining data for FFTW in real and fourier space (in place) @@ -172,7 +181,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 - debugRestart = iand(debug_level(debug_spectral),debug_spectralRestart) /= 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 @@ -187,12 +195,41 @@ subroutine utilities_init() call IO_warning(41_pInt, ext_msg='debug PETSc') #endif + call IO_open_file(fileUnit,geometryFile) ! parse info from geometry file... + grid = mesh_spectral_getGrid(fileUnit) + grid1Red = grid(1)/2_pInt + 1_pInt + wgt = 1.0/real(product(grid),pReal) + geomSize = mesh_spectral_getSize(fileUnit) + close(fileUnit) + + write(6,'(a,3(i12 ))') ' grid a b c: ', grid + write(6,'(a,3(f12.5))') ' size x y z: ', geomSize + +!-------------------------------------------------------------------------------------------------- +! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso- +! lution-independent divergence + if (divergence_correction == 1_pInt) then + do j = 1_pInt, 3_pInt + 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(geomSize/grid,1) .and. j /= maxloc(geomSize/grid,1)) & + scaledGeomSize = geomSize/geomSize(j)*grid(j) + enddo + else + scaledGeomSize = geomSize + endif + + + !-------------------------------------------------------------------------------------------------- ! allocation - allocate (xi(3,res1_red,res(2),res(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension - tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8) - call c_f_pointer(tensorField, field_real, [res(1)+2_pInt-mod(res(1),2_pInt),res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField - call c_f_pointer(tensorField, field_fourier,[res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField + allocate (xi(3,grid1Red,grid(2),grid(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension + tensorField = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8) + call c_f_pointer(tensorField, field_real, [grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3),3,3])! place a pointer for a real representation on tensorField + call c_f_pointer(tensorField, field_fourier,[grid1Red, grid(2),grid(3),3,3])! place a pointer for a complex representation on tensorField !-------------------------------------------------------------------------------------------------- ! general initialization of FFTW (see manual on fftw.org for more details) @@ -206,41 +243,41 @@ subroutine utilities_init() !-------------------------------------------------------------------------------------------------- ! creating plans for the convolution - planForth = fftw_plan_many_dft_r2c(3,[res(3),res(2) ,res(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms - field_real,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & ! input data, physical length in each dimension in reversed order - 1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions - field_fourier,[res(3),res(2) ,res1_red], & ! output data, physical length in each dimension in reversed order - 1, res(3)*res(2)* res1_red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision + planForth = fftw_plan_many_dft_r2c(3,[grid(3),grid(2) ,grid(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms + field_real,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], & ! input data, physical length in each dimension in reversed order + 1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions + field_fourier,[grid(3),grid(2) ,grid1Red], & ! output data, physical length in each dimension in reversed order + 1, grid(3)*grid(2)* grid1Red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision - planBack = fftw_plan_many_dft_c2r(3,[res(3),res(2) ,res(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms - field_fourier,[res(3),res(2) ,res1_red], & ! input data, physical length in each dimension in reversed order - 1, res(3)*res(2)* res1_red, & ! striding, product of physical length in the 3 dimensions - field_real,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & ! output data, physical length in each dimension in reversed order - 1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions + planBack = fftw_plan_many_dft_c2r(3,[grid(3),grid(2) ,grid(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms + field_fourier,[grid(3),grid(2) ,grid1Red], & ! input data, physical length in each dimension in reversed order + 1, grid(3)*grid(2)* grid1Red, & ! striding, product of physical length in the 3 dimensions + field_real,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], & ! output data, physical length in each dimension in reversed order + 1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions fftw_planner_flag) ! planner precision !-------------------------------------------------------------------------------------------------- ! depending on debug options, allocate more memory and create additional plans if (debugDivergence) then - div = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T)) - call c_f_pointer(div,divReal, [res(1)+2_pInt-mod(res(1),2_pInt),res(2),res(3),3]) - call c_f_pointer(div,divFourier,[res1_red, res(2),res(3),3]) - planDiv = fftw_plan_many_dft_c2r(3,[res(3),res(2) ,res(1)],3,& - divFourier,[res(3),res(2) ,res1_red],& - 1, res(3)*res(2)* res1_red,& - divReal,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & - 1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & + div = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3)*3_pInt,C_SIZE_T)) + call c_f_pointer(div,divReal, [grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3),3]) + call c_f_pointer(div,divFourier,[grid1Red, grid(2),grid(3),3]) + planDiv = fftw_plan_many_dft_c2r(3,[grid(3),grid(2) ,grid(1)],3,& + divFourier,[grid(3),grid(2) ,grid1Red],& + 1, grid(3)*grid(2)* grid1Red,& + divReal,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], & + 1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & fftw_planner_flag) endif if (debugFFTW) then - scalarField_realC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! allocate data for real representation (no in place transform) - scalarField_fourierC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! allocate data for fourier representation (no in place transform) - call c_f_pointer(scalarField_realC, scalarField_real, [res(1),res(2),res(3)]) ! place a pointer for a real representation - call c_f_pointer(scalarField_fourierC, scalarField_fourier, [res(1),res(2),res(3)]) ! place a pointer for a fourier representation - planDebugForth = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) + scalarField_realC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allocate data for real representation (no in place transform) + scalarField_fourierC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allocate data for fourier representation (no in place transform) + call c_f_pointer(scalarField_realC, scalarField_real, grid) ! place a pointer for a real representation + call c_f_pointer(scalarField_fourierC, scalarField_fourier, grid) ! place a pointer for a fourier representation + planDebugForth = fftw_plan_dft_3d(grid(3),grid(2),grid(1),& ! reversed order (C style) scalarField_real,scalarField_fourier,-1,fftw_planner_flag) ! input, output, forward FFT(-1), planner precision - planDebugBack = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) + planDebugBack = fftw_plan_dft_3d(grid(3),grid(2),grid(1),& ! reversed order (C style) scalarField_fourier,scalarField_real,+1,fftw_planner_flag) ! input, output, backward (1), planner precision endif @@ -249,21 +286,21 @@ subroutine utilities_init() !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - do k = 1_pInt, res(3) + do k = 1_pInt, grid(3) k_s(3) = k - 1_pInt - if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 - do j = 1_pInt, res(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 > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 - do i = 1_pInt, res1_red + 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) = real(k_s, pReal)/scaledDim ! if divergence_correction is set, frequencies are calculated on unit length + xi(1:3,i,j,k) = 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,res1_red ,res(2),res(3)), source =0.0_pReal) + allocate (gamma_hat(3,3,3,3,grid1Red ,grid(2),grid(3)), source = 0.0_pReal) endif end subroutine utilities_init @@ -284,9 +321,6 @@ subroutine utilities_updateGamma(C,saveReference) memory_efficient use math, only: & math_inv33 - use mesh, only: & - res, & - res1_red implicit none real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness @@ -307,7 +341,7 @@ subroutine utilities_updateGamma(C,saveReference) endif if(.not. memory_efficient) then - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red + do k = 1_pInt, grid(3); 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)*xi(m, i,j,k) @@ -333,10 +367,6 @@ end subroutine utilities_updateGamma !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTforward() !< @ToDo make row and column between randomly between 1 and 3 use math - use mesh, only : & - scaledDim, & - res, & - res1_red implicit none integer(pInt) :: row, column ! if debug FFTW, compare 3D array field of row and column @@ -349,7 +379,7 @@ subroutine utilities_FFTforward() 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) - scalarField_real = cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) ! store the selected component + scalarField_real = cmplx(field_real(1:grid(1),1:grid(2),1:grid(3),row,column),0.0_pReal,pReal) ! store the selected component endif !-------------------------------------------------------------------------------------------------- @@ -360,34 +390,34 @@ subroutine utilities_FFTforward() ! comparing 1 and 3x3 FT results if (debugFFTW) then call fftw_execute_dft(planDebugForth,scalarField_real,scalarField_fourier) - where(abs(scalarField_fourier(1:res1_red,1:res(2),1:res(3))) > tiny(1.0_pReal)) ! avoid division by zero - scalarField_fourier(1:res1_red,1:res(2),1:res(3)) = & - (scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& - field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& - scalarField_fourier(1:res1_red,1:res(2),1:res(3)) + where(abs(scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3))) > tiny(1.0_pReal)) ! avoid division by zero + scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3)) = & + (scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3))-& + field_fourier(1:grid1Red,1:grid(2),1:grid(3),row,column))/& + scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3)) else where scalarField_real = cmplx(0.0,0.0,pReal) end where 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 - maxval(real (scalarField_fourier(1:res1_red,1:res(2),1:res(3)))),& - maxval(aimag(scalarField_fourier(1:res1_red,1:res(2),1:res(3)))) + maxval(real (scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3)))),& + maxval(aimag(scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3)))) flush(6) endif !-------------------------------------------------------------------------------------------------- ! removing highest frequencies - Nyquist(2,1:2) = [res(2)/2_pInt + 1_pInt, res(2)/2_pInt + 1_pInt + mod(res(2),2_pInt)] - Nyquist(3,1:2) = [res(3)/2_pInt + 1_pInt, res(3)/2_pInt + 1_pInt + mod(res(3),2_pInt)] + Nyquist(2,1:2) = [grid(2)/2_pInt + 1_pInt, grid(2)/2_pInt + 1_pInt + mod(grid(2),2_pInt)] + Nyquist(3,1:2) = [grid(3)/2_pInt + 1_pInt, grid(3)/2_pInt + 1_pInt + mod(grid(3),2_pInt)] - if(res(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation - field_fourier (res1_red, 1:res(2), 1:res(3), 1:3,1:3) & + if(grid(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + field_fourier (grid1Red, 1:grid(2), 1:grid(3), 1:3,1:3) & = cmplx(0.0_pReal,0.0_pReal,pReal) - if(res(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation - field_fourier (1:res1_red,Nyquist(2,1):Nyquist(2,2),1:res(3), 1:3,1:3) & + if(grid(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + field_fourier (1:grid1Red,Nyquist(2,1):Nyquist(2,2),1:grid(3), 1:3,1:3) & = cmplx(0.0_pReal,0.0_pReal,pReal) - if(res(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation - field_fourier (1:res1_red,1:res(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) & + if(grid(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + field_fourier (1:grid1Red,1:grid(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) & = cmplx(0.0_pReal,0.0_pReal,pReal) end subroutine utilities_FFTforward @@ -401,10 +431,6 @@ end subroutine utilities_FFTforward !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTbackward() use math !< must use the whole module for use of FFTW - use mesh, only: & - wgt, & - res, & - res1_red implicit none integer(pInt) :: row, column !< if debug FFTW, compare 3D array field of row and column @@ -417,18 +443,18 @@ subroutine utilities_FFTbackward() 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) - scalarField_fourier(1:res1_red,1:res(2),1:res(3)) & - = field_fourier(1:res1_red,1:res(2),1:res(3),row,column) - do i = 0_pInt, res(1)/2_pInt-2_pInt + mod(res(1),2_pInt) + scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3)) & + = field_fourier(1:grid1Red,1:grid(2),1:grid(3),row,column) + do i = 0_pInt, grid(1)/2_pInt-2_pInt + mod(grid(1),2_pInt) m = 1_pInt - do k = 1_pInt, res(3) + do k = 1_pInt, grid(3) n = 1_pInt - do j = 1_pInt, res(2) - scalarField_fourier(res(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m)) - if(n == 1_pInt) n = res(2) + 1_pInt + do j = 1_pInt, grid(2) + scalarField_fourier(grid(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m)) + if(n == 1_pInt) n = grid(2) + 1_pInt n = n-1_pInt enddo - if(m == 1_pInt) m = res(3) + 1_pInt + if(m == 1_pInt) m = grid(3) + 1_pInt m = m -1_pInt enddo; enddo endif @@ -443,7 +469,7 @@ subroutine utilities_FFTbackward() call fftw_execute_dft(planDebugBack,scalarField_fourier,scalarField_real) where(abs(real(scalarField_real,pReal)) > tiny(1.0_pReal)) ! avoid division by zero scalarField_real = (scalarField_real & - - cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column), 0.0, pReal))/ & + - cmplx(field_real(1:grid(1),1:grid(2),1:grid(3),row,column), 0.0, pReal))/ & scalarField_real else where scalarField_real = cmplx(0.0,0.0,pReal) @@ -466,10 +492,6 @@ subroutine utilities_fourierConvolution(fieldAim) memory_efficient use math, only: & math_inv33 - use mesh, only: & - mesh_NcpElems, & - res, & - res1_red implicit none real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution @@ -486,7 +508,7 @@ subroutine utilities_fourierConvolution(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, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red + do k = 1_pInt, grid(3); 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)*xi(m, i,j,k) @@ -502,13 +524,13 @@ subroutine utilities_fourierConvolution(fieldAim) endif enddo; enddo; enddo else ! use precalculated gamma-operator - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt,res1_red + do k = 1_pInt, grid(3); 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) * field_fourier(i,j,k,1:3,1:3)) field_fourier(i,j,k, 1:3,1:3) = temp33_Complex enddo; enddo; enddo endif - field_fourier(1,1,1,1:3,1:3) = cmplx(fieldAim*real(mesh_NcpElems,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + field_fourier(1,1,1,1:3,1:3) = cmplx(fieldAim*real(product(grid),pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 end subroutine utilities_fourierConvolution @@ -518,10 +540,6 @@ end subroutine utilities_fourierConvolution !-------------------------------------------------------------------------------------------------- real(pReal) function utilities_divergenceRMS() use math !< must use the whole module for use of FFTW - use mesh, only: & - wgt, & - res, & - res1_red implicit none integer(pInt) :: i, j, k @@ -537,32 +555,32 @@ real(pReal) function utilities_divergenceRMS() !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space utilities_divergenceRMS = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2) - do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. + do k = 1_pInt, grid(3); 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(field_fourier(i,j,k,1:3,1:3),& ! (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(field_fourier(i,j,k,1:3,1:3),& 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 res(1) /= 1) + 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(field_fourier(1 ,j,k,1:3,1:3),& xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& + sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& - + sum( real(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& - xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)& - + sum(aimag(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& - xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) + + sum( real(math_mul33x3_complex(field_fourier(grid1Red,j,k,1:3,1:3),& + xi(1:3,grid1Red,j,k))*TWOPIIMG)**2.0_pReal)& + + sum(aimag(math_mul33x3_complex(field_fourier(grid1Red,j,k,1:3,1:3),& + xi(1:3,grid1Red,j,k))*TWOPIIMG)**2.0_pReal) enddo; enddo - if(res(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of res(1) == 1 + if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 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, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red + do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*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)) @@ -590,10 +608,6 @@ end function utilities_divergenceRMS !-------------------------------------------------------------------------------------------------- real(pReal) function utilities_curlRMS() use math !< must use the whole module for use of FFTW - use mesh, only: & - res, & - res1_red, & - wgt implicit none integer(pInt) :: i, j, k, l @@ -605,8 +619,8 @@ real(pReal) function utilities_curlRMS() ! calculating max curl criterion in Fourier space utilities_curlRMS = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); - do i = 2_pInt, res1_red - 1_pInt + do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); + do i = 2_pInt, grid1Red - 1_pInt do l = 1_pInt, 3_pInt curl_fourier(l,1) = (field_fourier(i,j,k,l,3)*xi(2,i,j,k)& - field_fourier(i,j,k,l,2)*xi(3,i,j,k))*TWOPIIMG @@ -629,12 +643,12 @@ real(pReal) function utilities_curlRMS() 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 = (field_fourier(res1_red,j,k,l,3)*xi(2,res1_red,j,k)& - - field_fourier(res1_red,j,k,l,2)*xi(3,res1_red,j,k))*TWOPIIMG - curl_fourier = (-field_fourier(res1_red,j,k,l,3)*xi(1,res1_red,j,k)& - +field_fourier(res1_red,j,k,l,1)*xi(3,res1_red,j,k) )*TWOPIIMG - curl_fourier = ( field_fourier(res1_red,j,k,l,2)*xi(1,res1_red,j,k)& - -field_fourier(res1_red,j,k,l,1)*xi(2,res1_red,j,k) )*TWOPIIMG + curl_fourier = ( field_fourier(grid1Red,j,k,l,3)*xi(2,grid1Red,j,k)& + -field_fourier(grid1Red,j,k,l,2)*xi(3,grid1Red,j,k))*TWOPIIMG + curl_fourier = (-field_fourier(grid1Red,j,k,l,3)*xi(1,grid1Red,j,k)& + +field_fourier(grid1Red,j,k,l,1)*xi(3,grid1Red,j,k))*TWOPIIMG + curl_fourier = ( field_fourier(grid1Red,j,k,l,2)*xi(1,grid1Red,j,k)& + -field_fourier(grid1Red,j,k,l,1)*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) @@ -758,10 +772,6 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& math_det33 use FEsolving, only: & restartWrite - use mesh, only: & - res, & - wgt, & - mesh_NcpElems use CPFEM, only: & CPFEM_general, & CPFEM_COLLECT, & @@ -778,16 +788,16 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& implicit none real(pReal), intent(inout) :: temperature !< temperature (no field) - real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid(3)) :: & F_lastInc, & !< target deformation gradient F !< previous deformation gradient real(pReal), intent(in) :: timeinc !< loading time logical, intent(in) :: forwardData !< age results - real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame + real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame - 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,res(1),res(2),res(3)) :: P !< PK stress + real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness + real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress + real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid(3)) :: P !< PK stress integer(pInt) :: & calcMode, & !< CPFEM mode for calculation @@ -812,8 +822,8 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& call CPFEM_general(collectMode,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), & ! collect mode handles Jacobian backup / restoration temperature,timeinc,1_pInt,1_pInt) - materialpoint_F0 = reshape(F_lastInc, [3,3,1,mesh_NcpElems]) - materialpoint_F = reshape(F, [3,3,1,mesh_NcpElems]) + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)]) + materialpoint_F = reshape(F, [3,3,1,product(grid)]) materialpoint_Temperature = temperature call debug_reset() @@ -823,7 +833,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& if(debugGeneral) then defgradDetMax = -huge(1.0_pReal) defgradDetMin = +huge(1.0_pReal) - do j = 1_pInt, mesh_NcpElems + do j = 1_pInt, product(grid) defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j)) defgradDetMax = max(defgradDetMax,defgradDet) defgradDetMin = min(defgradDetMin,defgradDet) @@ -840,7 +850,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& max_dPdF_norm = 0.0_pReal min_dPdF = huge(1.0_pReal) min_dPdF_norm = huge(1.0_pReal) - do k = 1_pInt, mesh_NcpElems + do k = 1_pInt, product(grid) 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) @@ -851,7 +861,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& endif end do - P = reshape(materialpoint_P, [3,3,res(1),res(2),res(3)]) + P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid(3)]) C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF) @@ -875,24 +885,22 @@ end subroutine utilities_constitutiveResponse !> @brief calculates forward rate, either guessing or just add delta/timeinc !-------------------------------------------------------------------------------------------------- pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,field) - use mesh, only: & - res - + implicit none real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon real(pReal), intent(in) :: & timeinc_old !< timeinc of last step logical, intent(in) :: & guess !< guess along former trajectory - real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid(3)) :: & field_lastInc, & !< data of previous step field !< data of current step - real(pReal), dimension(3,3,res(1),res(2),res(3)) :: utilities_calculateRate + real(pReal), dimension(3,3,grid(1),grid(2),grid(3)) :: utilities_calculateRate if(guess) then utilities_calculateRate = (field-field_lastInc) / timeinc_old else - utilities_calculateRate = spread(spread(spread(avRate,3,res(1)),4,res(2)),5,res(3)) + utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid(3)) endif end function utilities_calculateRate @@ -903,26 +911,23 @@ end function utilities_calculateRate !> ensures that the average matches the aim !-------------------------------------------------------------------------------------------------- pure function utilities_forwardField(timeinc,field_lastInc,rate,aim) - use mesh, only: & - res, & - wgt implicit none real(pReal), intent(in) :: & timeinc !< timeinc of current step - real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid(3)) :: & 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,res(1),res(2),res(3)) :: utilities_forwardField + real(pReal), dimension(3,3,grid(1),grid(2),grid(3)) :: utilities_forwardField real(pReal), dimension(3,3) :: fieldDiff !< - aim 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 - aim utilities_forwardField = utilities_forwardField - & - spread(spread(spread(fieldDiff,3,res(1)),4,res(2)),5,res(3)) + spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid(3)) endif end function utilities_forwardField @@ -936,8 +941,6 @@ real(pReal) function utilities_getFilter(k) IO_error use numerics, only: & myfilter - use mesh, only: & - res use math, only: & PI @@ -949,9 +952,9 @@ real(pReal) function utilities_getFilter(k) select case (myfilter) case ('none') !< default is already nothing (1.0_pReal) case ('cosine') !< cosine curve with 1 for avg and zero for highest freq - utilities_getFilter = (1.0_pReal + cos(PI*k(3)/res(3))) & - *(1.0_pReal + cos(PI*k(2)/res(2))) & - *(1.0_pReal + cos(PI*k(1)/res(1)))/8.0_pReal + utilities_getFilter = (1.0_pReal + cos(PI*k(3)/grid(3))) & + *(1.0_pReal + cos(PI*k(2)/grid(2))) & + *(1.0_pReal + cos(PI*k(1)/grid(1)))/8.0_pReal case default call IO_error(error_ID = 892_pInt, ext_msg = trim(myfilter)) end select diff --git a/code/lattice.f90 b/code/lattice.f90 index f3956e11b..89238e45f 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -1018,7 +1018,8 @@ pure function lattice_symmetrizeC66(structName,C66) ! TwinTwinInteraction !-------------------------------------------------------------------------------------------------- function lattice_configNchunks(struct) - use prec, only: pReal,pInt + use prec, only: & + pInt implicit none integer(pInt), dimension(6) :: lattice_configNchunks diff --git a/code/math.f90 b/code/math.f90 index f0daa583f..48307503b 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -727,14 +727,14 @@ pure function math_exp33(A,n) order = 5 if (present(n)) order = n - B = math_identity2nd(3) ! init - invfac = 1.0_pReal ! 0! - math_exp33 = B ! A^0 = eye2 + B = math_identity2nd(3) ! init + invfac = 1.0_pReal ! 0! + math_exp33 = B ! A^0 = eye2 do i = 1_pInt,n - invfac = invfac/real(i) ! invfac = 1/i! + invfac = invfac/real(i) ! invfac = 1/i! B = math_mul33x33(B,A) - math_exp33 = math_exp33 + invfac*B ! exp = SUM (A^i)/i! + math_exp33 = math_exp33 + invfac*B ! exp = SUM (A^i)/i! enddo end function math_exp33 diff --git a/code/mesh.f90 b/code/mesh.f90 index 5f9205279..53771fb0a 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -121,16 +121,6 @@ module mesh #ifdef Spectral include 'fftw3.f03' - real(pReal), dimension(3), public, protected :: & - geomdim, & !< physical dimension of volume element per direction - scaledDim !< scaled dimension of volume element, depending on selected divergence calculation - integer(pInt), dimension(3), public, protected :: & - res !< resolution, e.g. number of Fourier points in each direction - real(pReal), public, protected :: & - wgt - integer(pInt), public, protected :: & - res1_red, & - homog #endif ! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) @@ -428,6 +418,8 @@ module mesh mesh_get_nodeAtIP #ifdef Spectral public :: & + mesh_spectral_getGrid, & + mesh_spectral_getSize, & mesh_regrid, & mesh_nodesAroundCentres, & mesh_deformedCoordsFFT, & @@ -438,13 +430,9 @@ module mesh private :: & #ifdef Spectral - mesh_spectral_getGrid, & - mesh_spectral_getSize, & mesh_spectral_getHomogenization, & - mesh_spectral_count_nodesAndElements, & - mesh_spectral_count_cpElements, & - mesh_spectral_map_elements, & - mesh_spectral_map_nodes, & + mesh_spectral_count, & + mesh_spectral_mapNodesAndElems, & mesh_spectral_count_cpSizes, & mesh_spectral_build_nodes, & mesh_spectral_build_elements, & @@ -554,51 +542,23 @@ subroutine mesh_init(ip,el) if (allocated(FE_ipNeighbor)) deallocate(FE_ipNeighbor) if (allocated(FE_cellnodeParentnodeWeights)) deallocate(FE_cellnodeParentnodeWeights) if (allocated(FE_subNodeOnIPFace)) deallocate(FE_subNodeOnIPFace) - call mesh_build_FEdata ! get properties of the different types of elements mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh #ifdef Spectral call IO_open_file(fileUnit,geometryFile) ! parse info from geometry file... - res = mesh_spectral_getGrid(fileUnit) - res1_red = res(1)/2_pInt + 1_pInt - wgt = 1.0/real(product(res),pReal) - geomdim = mesh_spectral_getSize(fileUnit) - homog = mesh_spectral_getHomogenization(fileUnit) - -!-------------------------------------------------------------------------------------------------- -! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso- -! lution-independent divergence - if (divergence_correction == 1_pInt) then - do j = 1_pInt, 3_pInt - if (j /= minloc(geomdim,1) .and. j /= maxloc(geomdim,1)) scaledDim = geomdim/geomdim(j) - enddo - elseif (divergence_correction == 2_pInt) then - do j = 1_pInt, 3_pInt - if (j /= minloc(geomdim/res,1) .and. j /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(j)*res(j) - enddo - else - scaledDim = geomdim - endif - write(6,'(a,3(i12 ))') ' grid a b c: ', res - write(6,'(a,3(f12.5))') ' size x y z: ', geomdim - write(6,'(a,i5,/)') ' homogenization: ', homog - - call mesh_spectral_count_nodesAndElements - call mesh_spectral_count_cpElements - call mesh_spectral_map_elements - call mesh_spectral_map_nodes + call mesh_spectral_count(fileUnit) + call mesh_spectral_mapNodesAndElems call mesh_spectral_count_cpSizes - call mesh_spectral_build_nodes + call mesh_spectral_build_nodes(fileUnit) call mesh_spectral_build_elements(fileUnit) call mesh_get_damaskOptions(fileUnit) - close (fileUnit) call mesh_build_cellconnectivity mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) call mesh_build_ipCoordinates call mesh_build_ipVolumes call mesh_build_ipAreas - call mesh_spectral_build_ipNeighborhood + call mesh_spectral_build_ipNeighborhood(fileUnit) #endif #ifdef Marc4DAMASK call IO_open_inputFile(fileUnit,modelName) ! parse info from input file... @@ -613,7 +573,6 @@ subroutine mesh_init(ip,el) call mesh_marc_count_cpSizes(fileunit) call mesh_marc_build_elements(fileUnit) call mesh_get_damaskOptions(fileUnit) - close (fileUnit) call mesh_build_cellconnectivity mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) call mesh_build_ipCoordinates @@ -638,7 +597,6 @@ subroutine mesh_init(ip,el) call mesh_abaqus_count_cpSizes(fileunit) call mesh_abaqus_build_elements(fileUnit) call mesh_get_damaskOptions(fileUnit) - close (fileUnit) call mesh_build_cellconnectivity mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) call mesh_build_ipCoordinates @@ -649,6 +607,7 @@ subroutine mesh_init(ip,el) call mesh_build_ipNeighborhood #endif + close (fileUnit) call mesh_tell_statistics call mesh_write_meshfile call mesh_write_cellGeom @@ -950,25 +909,24 @@ end subroutine mesh_build_ipCoordinates !-------------------------------------------------------------------------------------------------- pure function mesh_cellCenterCoordinates(ip,el) -implicit none + implicit none + integer(pInt), intent(in) :: el, & !< element number + ip !< integration point number + real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell -integer(pInt), intent(in) :: el, & !< element number - ip !< integration point number -real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell - -integer(pInt) :: t,g,c,n + integer(pInt) :: t,g,c,n -t = mesh_element(2_pInt,el) ! get element type -g = FE_geomtype(t) ! get geometry type -c = FE_celltype(g) ! get cell type -mesh_cellCenterCoordinates = 0.0_pReal -do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell - mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el)) -enddo -mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / FE_NcellnodesPerCell(c) + t = mesh_element(2_pInt,el) ! get element type + g = FE_geomtype(t) ! get geometry type + c = FE_celltype(g) ! get cell type + mesh_cellCenterCoordinates = 0.0_pReal + do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell + mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el)) + enddo + mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / FE_NcellnodesPerCell(c) -endfunction mesh_cellCenterCoordinates + end function mesh_cellCenterCoordinates #ifdef Spectral @@ -1189,61 +1147,39 @@ end function mesh_spectral_getHomogenization !-------------------------------------------------------------------------------------------------- !> @brief Count overall number of nodes and elements in mesh and stores them in -!! 'mesh_Nelems' and 'mesh_Nnodes' +!! 'mesh_Nelems', 'mesh_Nnodes' and 'mesh_NcpElems' !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_count_nodesAndElements() +subroutine mesh_spectral_count(myUnit) implicit none - mesh_Nelems = product(res) - mesh_Nnodes = product(res+1_pInt) - -end subroutine mesh_spectral_count_nodesAndElements - - -!-------------------------------------------------------------------------------------------------- -!> @brief Count overall number of CP elements in mesh and stores them in 'mesh_NcpElems' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_count_cpElements - - implicit none + integer(pInt), intent(in) :: myUnit + integer(pInt), dimension(3) :: grid + grid = mesh_spectral_getGrid(myUnit) + mesh_Nelems = product(grid) mesh_NcpElems = mesh_Nelems - -end subroutine mesh_spectral_count_cpElements + mesh_Nnodes = product(grid+1_pInt) + +end subroutine mesh_spectral_count !-------------------------------------------------------------------------------------------------- -!> @brief Maps elements from FE ID to internal (consecutive) representation. -!! Allocates global array 'mesh_mapFEtoCPelem' +!> @brief fake map node from FE ID to internal (consecutive) representation for node and element +!! Allocates global array 'mesh_mapFEtoCPnode' and 'mesh_mapFEtoCPelem' !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_map_elements +subroutine mesh_spectral_mapNodesAndElems + use math, only: & + math_range implicit none - integer(pInt) :: i - allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt + allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes), source = 0_pInt) + allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems), source = 0_pInt) - forall (i = 1_pInt:mesh_NcpElems) & - mesh_mapFEtoCPelem(1:2,i) = i + mesh_mapFEtoCPnode = spread(math_range(mesh_Nnodes),1,2) + mesh_mapFEtoCPelem = spread(math_range(mesh_NcpElems),1,2) -end subroutine mesh_spectral_map_elements - - -!-------------------------------------------------------------------------------------------------- -!> @brief Maps node from FE ID to internal (consecutive) representation. -!! Allocates global array 'mesh_mapFEtoCPnode' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_map_nodes - - implicit none - integer(pInt) :: i - - allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - - forall (i = 1_pInt:mesh_Nnodes) & - mesh_mapFEtoCPnode(1:2,i) = i - -end subroutine mesh_spectral_map_nodes +end subroutine mesh_spectral_mapNodesAndElems !-------------------------------------------------------------------------------------------------- @@ -1272,24 +1208,30 @@ end subroutine mesh_spectral_count_cpSizes !> @brief Store x,y,z coordinates of all nodes in mesh. !! Allocates global arrays 'mesh_node0' and 'mesh_node' !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_build_nodes() +subroutine mesh_spectral_build_nodes(myUnit) implicit none - integer(pInt) :: n + integer(pInt), intent(in) :: myUnit + integer(pInt) :: n + integer(pInt), dimension(3) :: grid + real(pReal), dimension(3) :: geomSize - allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal - allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal + allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal) + allocate (mesh_node (3,mesh_Nnodes), source = 0.0_pReal) + grid = mesh_spectral_getGrid(myUnit) + geomSize = mesh_spectral_getSize(myUnit) + forall (n = 0_pInt:mesh_Nnodes-1_pInt) mesh_node0(1,n+1_pInt) = mesh_unitlength * & - geomdim(1) * real(mod(n,(res(1)+1_pInt) ),pReal) & - / real(res(1),pReal) + geomSize(1)*real(mod(n,(grid(1)+1_pInt) ),pReal) & + / real(grid(1),pReal) mesh_node0(2,n+1_pInt) = mesh_unitlength * & - geomdim(2) * real(mod(n/(res(1)+1_pInt),(res(2)+1_pInt)),pReal) & - / real(res(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 * & - geomdim(3) * real(mod(n/(res(1)+1_pInt)/(res(2)+1_pInt),(res(3)+1_pInt)),pReal) & - / real(res(3),pReal) + geomSize(3)*real(mod(n/(grid(1)+1_pInt)/(grid(2)+1_pInt),(grid(3)+1_pInt)),pReal) & + / real(grid(3),pReal) end forall mesh_node = mesh_node0 @@ -1314,15 +1256,29 @@ subroutine mesh_spectral_build_elements(myUnit) IO_countContinuousIntValues implicit none - integer(pInt), intent(in) :: myUnit - - integer(pInt), dimension (1_pInt+7_pInt*2_pInt) :: myPos - integer(pInt) :: e, i, headerLength = 0_pInt, maxIntCount - integer(pInt), dimension(:), allocatable :: microstructures - integer(pInt), dimension(1,1) :: dummySet = 0_pInt - character(len=65536) :: line,keyword - character(len=64), dimension(1) :: dummyName = '' + integer(pInt), intent(in) :: & + myUnit + integer(pInt), dimension(1_pInt+7_pInt*2_pInt) :: & + myPos + integer(pInt) :: & + e, i, & + headerLength = 0_pInt, & + maxIntCount, & + homog + integer(pInt), dimension(:), allocatable :: & + microstructures + integer(pInt), dimension(3) :: & + grid + integer(pInt), dimension(1,1) :: & + dummySet = 0_pInt + character(len=65536) :: & + line, & + keyword + character(len=64), dimension(1) :: & + dummyName = '' + grid = mesh_spectral_getGrid(myUnit) + homog = mesh_spectral_getHomogenization(myUnit) call IO_checkAndRewind(myUnit) read(myUnit,'(a65536)') line @@ -1364,15 +1320,15 @@ subroutine mesh_spectral_build_elements(myUnit) mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type mesh_element( 3,e) = homog ! homogenization mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure - mesh_element( 5,e) = e + (e-1_pInt)/res(1) + & - ((e-1_pInt)/(res(1)*res(2)))*(res(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) + res(1) + 2_pInt - mesh_element( 8,e) = mesh_element(5,e) + res(1) + 1_pInt - mesh_element( 9,e) = mesh_element(5,e) +(res(1) + 1_pInt) * (res(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) + res(1) + 2_pInt - mesh_element(12,e) = mesh_element(9,e) + res(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 @@ -1388,56 +1344,59 @@ end subroutine mesh_spectral_build_elements !> @brief build neighborhood relations for spectral !> @details assign globals: mesh_ipNeighborhood !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_build_ipNeighborhood +subroutine mesh_spectral_build_ipNeighborhood(myUnit) -implicit none -integer(pInt) x,y,z, & - e - -allocate(mesh_ipNeighborhood(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) -mesh_ipNeighborhood = 0_pInt - - -e = 0_pInt -do z = 0_pInt,res(3)-1_pInt - do y = 0_pInt,res(2)-1_pInt - do x = 0_pInt,res(1)-1_pInt - e = e + 1_pInt - mesh_ipNeighborhood(1,1,1,e) = z * res(1) * res(2) & - + y * res(1) & - + modulo(x+1_pInt,res(1)) & - + 1_pInt - mesh_ipNeighborhood(1,2,1,e) = z * res(1) * res(2) & - + y * res(1) & - + modulo(x-1_pInt,res(1)) & - + 1_pInt - mesh_ipNeighborhood(1,3,1,e) = z * res(1) * res(2) & - + modulo(y+1_pInt,res(2)) * res(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(1,4,1,e) = z * res(1) * res(2) & - + modulo(y-1_pInt,res(2)) * res(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,res(3)) * res(1) * res(2) & - + y * res(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,res(3)) * res(1) * res(2) & - + y * res(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt - mesh_ipNeighborhood(3,1,1,e) = 2_pInt - mesh_ipNeighborhood(3,2,1,e) = 1_pInt - mesh_ipNeighborhood(3,3,1,e) = 4_pInt - mesh_ipNeighborhood(3,4,1,e) = 3_pInt - mesh_ipNeighborhood(3,5,1,e) = 6_pInt - mesh_ipNeighborhood(3,6,1,e) = 5_pInt - enddo - enddo -enddo + implicit none + integer(pInt), intent(in) :: & + myUnit + integer(pInt) :: & + x,y,z, & + e + integer(pInt), dimension(3) :: & + grid + allocate(mesh_ipNeighborhood(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems),source=0_pInt) + + grid = mesh_spectral_getGrid(myUnit) + e = 0_pInt + do z = 0_pInt,grid(1)-1_pInt + do y = 0_pInt,grid(2)-1_pInt + do x = 0_pInt,grid(3)-1_pInt + e = e + 1_pInt + 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 * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x-1_pInt,grid(1)) & + + 1_pInt + 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 * 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,grid(3)) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1_pInt + mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,grid(3)) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1_pInt + mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt + mesh_ipNeighborhood(3,1,1,e) = 2_pInt + mesh_ipNeighborhood(3,2,1,e) = 1_pInt + mesh_ipNeighborhood(3,3,1,e) = 4_pInt + mesh_ipNeighborhood(3,4,1,e) = 3_pInt + mesh_ipNeighborhood(3,5,1,e) = 6_pInt + mesh_ipNeighborhood(3,6,1,e) = 5_pInt + enddo + enddo + enddo end subroutine mesh_spectral_build_ipNeighborhood @@ -1454,6 +1413,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) getSolverJobName, & GeometryFile use IO, only: & + IO_open_file, & IO_read_jobBinaryFile ,& IO_read_jobBinaryIntFile ,& IO_write_jobBinaryFile, & @@ -1467,16 +1427,17 @@ function mesh_regrid(adaptive,resNewInput,minRes) math_mul33x3 implicit none - character(len=1024):: formatString, N_Digits logical, intent(in) :: adaptive ! if true, choose adaptive grid based on resNewInput, otherwise keep it constant integer(pInt), dimension(3), optional, intent(in) :: resNewInput ! f2py cannot handle optional arguments correctly (they are always present) integer(pInt), dimension(3), optional, intent(in) :: minRes - integer(pInt), dimension(3) :: mesh_regrid, ratio + integer(pInt), dimension(3) :: mesh_regrid, ratio, grid + integer(pInt), parameter :: myUnit = 777_pInt integer(pInt), dimension(3,2) :: possibleResNew - integer(pInt):: maxsize, i, j, k, ielem, NpointsNew, spatialDim + integer(pInt):: maxsize, i, j, k, ielem, NpointsNew, spatialDim, Nelems integer(pInt), dimension(3) :: resNew integer(pInt), dimension(:), allocatable :: indices - real(pReal), dimension(3) :: geomdimNew + real(pReal) :: wgt + real(pReal), dimension(3) :: geomSizeNew, geomSize real(pReal), dimension(3,3) :: Favg, Favg_LastInc, & FavgNew, Favg_LastIncNew, & deltaF, deltaF_lastInc @@ -1497,13 +1458,21 @@ function mesh_regrid(adaptive,resNewInput,minRes) F_lastIncNew real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: & dPdF, dPdFNew - + character(len=1024):: formatString, N_Digits integer(pInt), dimension(:,:), allocatable :: & sizeStateHomog integer(pInt), dimension(:,:,:), allocatable :: & material_phase, material_phaseNew, & sizeStateConst - + + call IO_open_file(myUnit,trim(geometryFile)) + grid = mesh_spectral_getGrid(myUnit) + geomSize = mesh_spectral_getsize(myUnit) + close(myUnit) + + Nelems = product(grid) + wgt = 1.0_pReal/real(Nelems,pReal) + write(6,'(a)') 'Regridding geometry' if (adaptive) then write(6,'(a)') 'adaptive resolution determination' @@ -1519,53 +1488,54 @@ function mesh_regrid(adaptive,resNewInput,minRes) endif endif - allocate(coordinates(3,mesh_NcpElems)) + allocate(coordinates(3,Nelems)) + !-------------------------------------------------------------------------------------------------- ! read in deformation gradient to calculate coordinates, shape depend of selected solver select case(myspectralsolver) case('basic') - allocate(spectralF33(3,3,res(1),res(2),res(3))) + allocate(spectralF33(3,3,grid(1),grid(2),grid(3))) call IO_read_jobBinaryFile(777,'F',trim(getSolverJobName()),size(spectralF33)) read (777,rec=1) spectralF33 close (777) Favg = sum(sum(sum(spectralF33,dim=5),dim=4),dim=3) * wgt - coordinates = reshape(mesh_deformedCoordsFFT(geomdim,spectralF33),[3,mesh_NcpElems]) + coordinates = reshape(mesh_deformedCoordsFFT(geomSize,spectralF33),[3,mesh_NcpElems]) case('basicpetsc','al') - allocate(spectralF9(9,res(1),res(2),res(3))) + allocate(spectralF9(9,grid(1),grid(2),grid(3))) call IO_read_jobBinaryFile(777,'F',trim(getSolverJobName()),size(spectralF9)) read (777,rec=1) spectralF9 close (777) Favg = reshape(sum(sum(sum(spectralF9,dim=4),dim=3),dim=2) * wgt, [3,3]) - coordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(spectralF9, & - [3,3,res(1),res(2),res(3)])),[3,mesh_NcpElems]) + coordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(spectralF9, & + [3,3,grid(1),grid(2),grid(3)])),[3,mesh_NcpElems]) end select !-------------------------------------------------------------------------------------------------- ! sanity check 2D/3D case - if (res(3)== 1_pInt) then + if (grid(3)== 1_pInt) then spatialDim = 2_pInt if (present (minRes)) then if (minRes(1) > 0_pInt .or. minRes(2) > 0_pInt) then if (minRes(3) /= 1_pInt .or. & mod(minRes(1),2_pInt) /= 0_pInt .or. & - mod(minRes(2),2_pInt) /= 0_pInt) call IO_error(890_pInt, ext_msg = '2D minRes') ! as f2py has problems with present, use pyf file for initialization to -1 + mod(minRes(2),2_pInt) /= 0_pInt) call IO_error(890_pInt, ext_msg = '2D minRes') ! as f2py has problems with present, use pyf file for initialization to -1 endif; endif else spatialDim = 3_pInt if (present (minRes)) then if (any(minRes > 0_pInt)) then - if (mod(minRes(1),2_pInt) /= 0_pInt.or. & + if (mod(minRes(1),2_pInt) /= 0_pInt .or. & mod(minRes(2),2_pInt) /= 0_pInt .or. & - mod(minRes(3),2_pInt) /= 0_pInt) call IO_error(890_pInt, ext_msg = '3D minRes') ! as f2py has problems with present, use pyf file for initialization to -1 + mod(minRes(3),2_pInt) /= 0_pInt) call IO_error(890_pInt, ext_msg = '3D minRes') ! as f2py has problems with present, use pyf file for initialization to -1 endif; endif endif !-------------------------------------------------------------------------------------------------- ! Automatic detection based on current geom - geomdimNew = math_mul33x3(Favg,geomdim) + geomSizeNew = math_mul33x3(Favg,geomSize) if (adaptive) then - ratio = floor(real(resNewInput,pReal) * (geomdimNew/geomdim), pInt) + ratio = floor(real(resNewInput,pReal) * (geomSizeNew/geomSize), pInt) possibleResNew = 1_pInt do i = 1_pInt, spatialDim @@ -1574,15 +1544,14 @@ function mesh_regrid(adaptive,resNewInput,minRes) else possibleResNew(i,1:2) = [ratio(i)-1_pInt, ratio(i) + 1_pInt] endif - if (.not.present(minRes)) then ! calling from fortran, optional argument not given + if (.not.present(minRes)) then ! calling from fortran, optional argument not given possibleResNew = possibleResNew - else ! optional argument is there + else ! optional argument is there if (any(minRes<1_pInt)) then - possibleResNew = possibleResNew ! f2py calling, but without specification (or choosing invalid values), standard from pyf = -1 - else ! given useful values - do k = 1_pInt,3_pInt; do j = 1_pInt,3_pInt - possibleResNew(j,k) = max(possibleResNew(j,k), minRes(j)) - enddo; enddo + possibleResNew = possibleResNew ! f2py calling, but without specification (or choosing invalid values), standard from pyf = -1 + else ! given useful values + forall(k = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & + possibleResNew(j,k) = max(possibleResNew(j,k), minRes(j)) endif endif enddo @@ -1602,11 +1571,11 @@ function mesh_regrid(adaptive,resNewInput,minRes) endif enddo else - resNew = res + resNew = grid endif mesh_regrid = resNew - NpointsNew = resNew(1)*resNew(2)*resNew(3) + NpointsNew = product(resNew) !-------------------------------------------------------------------------------------------------- ! Calculate regular new coordinates @@ -1614,14 +1583,14 @@ function mesh_regrid(adaptive,resNewInput,minRes) ielem = 0_pInt do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) ielem = ielem + 1_pInt - coordinatesNew(1:3,ielem) = math_mul33x3(Favg, geomdim/real(resNew,pReal)*real([i,j,k],pReal) & - - geomdim/real(2_pInt*resNew,pReal)) + coordinatesNew(1:3,ielem) = math_mul33x3(Favg, geomSize/real(resNew,pReal)*real([i,j,k],pReal) & + - geomSize/real(2_pInt*resNew,pReal)) enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! Nearest neighbour search allocate(indices(NpointsNew)) - indices = math_periodicNearestNeighbor(geomdim, Favg, coordinatesNew, coordinates) + indices = math_periodicNearestNeighbor(geomSize, Favg, coordinatesNew, coordinates) deallocate(coordinates) !-------------------------------------------------------------------------------------------------- @@ -1665,8 +1634,8 @@ function mesh_regrid(adaptive,resNewInput,minRes) formatString = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//',a)' open(777,file=trim(getSolverWorkingDirectoryName())//trim(GeometryFile),status='REPLACE') write(777, '(A)') '3 header' - write(777, '(3(A, I8))') 'resolution a ', resNew(1), ' b ', resNew(2), ' c ', resNew(3) - write(777, '(3(A, g17.10))') 'dimension x ', geomdim(1), ' y ', geomdim(2), ' z ', geomdim(3) + write(777, '(3(A, I8))') 'grid a ', resNew(1), ' b ', resNew(2), ' c ', resNew(3) + write(777, '(3(A, g17.10))') 'size x ', geomSize(1), ' y ', geomSize(2), ' z ', geomSize(3) write(777, '(A)') 'homogenization 1' do i = 1_pInt, NpointsNew write(777,trim(formatString),advance='no') mesh_element(4,indices(i)), ' ' @@ -2763,11 +2732,12 @@ end subroutine mesh_marc_map_nodes !-------------------------------------------------------------------------------------------------- subroutine mesh_marc_build_nodes(myUnit) - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_fixedIntValue, & - IO_fixedNoEFloatValue + use IO, only: & + IO_lc, & + IO_stringValue, & + IO_stringPos, & + IO_fixedIntValue, & + IO_fixedNoEFloatValue implicit none integer(pInt), intent(in) :: myUnit @@ -3464,7 +3434,7 @@ subroutine mesh_abaqus_build_nodes(myUnit) myPos = IO_stringPos(line,maxNchunks) m = mesh_FEasCP('node',IO_intValue(line,myPos,1_pInt)) do j=1_pInt, 3_pInt - mesh_node0(j,m) = mesh_unitlength * IO_floatValue(line,myPos,j+1_pInt) + mesh_node0(j,m) = numerics_unitlength * IO_floatValue(line,myPos,j+1_pInt) enddo enddo endif @@ -3705,9 +3675,10 @@ use IO, only: & endselect endif enddo -#endif 610 FORMAT(A300) +#endif + 620 end subroutine mesh_get_damaskOptions @@ -3715,19 +3686,17 @@ use IO, only: & !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- subroutine mesh_build_ipAreas - use math, only: & math_norm3, & math_vectorproduct - implicit none integer(pInt) :: e,t,g,c,i,f,n,m real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals real(pReal), dimension(3) :: normal - allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal - allocate(mesh_ipAreaNormal(3_pInt,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal + allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)); mesh_ipArea = 0.0_pReal + allocate(mesh_ipAreaNormal(3_pInt,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)); mesh_ipAreaNormal = 0.0_pReal !$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals) do e = 1_pInt,mesh_NcpElems ! loop over cpElems