diff --git a/src/material.f90 b/src/material.f90 index 2b83c9967..d71fbb37a 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -692,7 +692,7 @@ subroutine material_parseMicrostructure implicit none character(len=65536), dimension(:), allocatable :: & - str + strings integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: e, m, c, i character(len=65536) :: & @@ -719,21 +719,22 @@ subroutine material_parseMicrostructure allocate(microstructure_texture (microstructure_maxNconstituents,size(config_microstructure)),source=0_pInt) allocate(microstructure_fraction(microstructure_maxNconstituents,size(config_microstructure)),source=0.0_pReal) + allocate(strings(1)) ! Intel 16.0 Bug do m=1_pInt, size(config_microstructure) - str = config_microstructure(m)%getStrings('(constituent)',raw=.true.) - do c = 1_pInt, size(str) - chunkPos = IO_stringPos(str(c)) + strings = config_microstructure(m)%getStrings('(constituent)',raw=.true.) + do c = 1_pInt, size(strings) + chunkPos = IO_stringPos(strings(c)) do i = 1_pInt,5_pInt,2_pInt - tag = IO_stringValue(str(c),chunkPos,i) + tag = IO_stringValue(strings(c),chunkPos,i) select case (tag) case('phase') - microstructure_phase(c,m) = IO_intValue(str(c),chunkPos,i+1_pInt) + microstructure_phase(c,m) = IO_intValue(strings(c),chunkPos,i+1_pInt) case('texture') - microstructure_texture(c,m) = IO_intValue(str(c),chunkPos,i+1_pInt) + microstructure_texture(c,m) = IO_intValue(strings(c),chunkPos,i+1_pInt) case('fraction') - microstructure_fraction(c,m) = IO_floatValue(str(c),chunkPos,i+1_pInt) + microstructure_fraction(c,m) = IO_floatValue(strings(c),chunkPos,i+1_pInt) end select enddo diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index b2b5e8173..9911eef21 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -62,8 +62,6 @@ module spectral_mech_basic integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment - real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal - public :: & basic_init, & basic_solution, & @@ -234,12 +232,12 @@ type(tSolutionState) function basic_solution(incInfoIn,timeinc,timeinc_old,stres !-------------------------------------------------------------------------------------------------- ! input data for solution - character(len=*), intent(in) :: & + character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & timeinc, & !< increment time for current solution timeinc_old !< increment time of last successful increment - type(tBoundaryCondition), intent(in) :: & + type(tBoundaryCondition), intent(in) :: & stress_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC @@ -261,7 +259,7 @@ type(tSolutionState) function basic_solution(incInfoIn,timeinc,timeinc_old,stres !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - mask_stress = stress_BC%maskFloat + params%stress_mask = stress_BC%maskFloat params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC params%timeinc = timeinc @@ -287,7 +285,11 @@ end function basic_solution !-------------------------------------------------------------------------------------------------- !> @brief forms the basic residual vector !-------------------------------------------------------------------------------------------------- -subroutine Basic_formResidual(in,x_scal,f_scal,dummy,ierr) +subroutine Basic_formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work) + F, & ! defgrad field on grid + residuum, & ! residuum field on grid + dummy, & + ierr) use numerics, only: & itmax, & itmin @@ -316,9 +318,9 @@ subroutine Basic_formResidual(in,x_scal,f_scal,dummy,ierr) implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in PetscScalar, & - dimension(3,3, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal !< what is this? + dimension(3,3, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: F PetscScalar, & - dimension(3,3, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: f_scal !< what is this? + dimension(3,3, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: residuum PetscInt :: & PETScIter, & nfuncs @@ -347,28 +349,29 @@ subroutine Basic_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(f_scal,P_av,C_volAvg,C_minMaxAvg, & - x_scal,params%timeinc, params%rotation_BC) + call Utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) + P_av,C_volAvg,C_minMaxAvg, & + F,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = f_scal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform call utilities_FFTtensorForward() ! FFT forward of global "tensorField_real" - err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier + err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg call utilities_FFTtensorBackward() ! FFT backward of global tensorField_fourier !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too + residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too end subroutine Basic_formResidual @@ -391,9 +394,9 @@ subroutine Basic_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,i SNES :: snes_local PetscInt :: PETScIter PetscReal :: & - xnorm, & - snorm, & - fnorm + xnorm, & ! not used + snorm, & ! not used + fnorm ! not used SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr @@ -457,16 +460,16 @@ subroutine Basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,s restartWrite implicit none - logical, intent(in) :: & + logical, intent(in) :: & guess - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & timeinc_old, & timeinc, & loadCaseTime !< remaining time of current load case type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - real(pReal), dimension(3,3), intent(in) ::& + real(pReal), dimension(3,3), intent(in) :: & rotation_BC PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: F diff --git a/src/spectral_mech_Polarisation.f90 b/src/spectral_mech_Polarisation.f90 index 9e567f0c9..9518452fa 100644 --- a/src/spectral_mech_Polarisation.f90 +++ b/src/spectral_mech_Polarisation.f90 @@ -27,7 +27,6 @@ module spectral_mech_Polarisation !-------------------------------------------------------------------------------------------------- ! derived types type(tSolutionParams), private :: params - real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -289,7 +288,7 @@ type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_ol !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - mask_stress = stress_BC%maskFloat + params%stress_mask = stress_BC%maskFloat params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC params%timeinc = timeinc @@ -315,7 +314,11 @@ end function Polarisation_solution !-------------------------------------------------------------------------------------------------- !> @brief forms the Polarisation residual vector !-------------------------------------------------------------------------------------------------- -subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) +subroutine Polarisation_formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work) + FandF_tau, & ! defgrad fields on grid + residuum, & ! residuum fields on grid + dummy, & + ierr) use numerics, only: & itmax, & itmin, & @@ -352,9 +355,9 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in PetscScalar, & - target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal + target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: FandF_tau PetscScalar, & - target, dimension(3,3,2, X_RANGE, Y_RANGE, Z_RANGE), intent(out) :: f_scal + target, dimension(3,3,2, X_RANGE, Y_RANGE, Z_RANGE), intent(out) :: residuum PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & F_tau, & @@ -368,20 +371,20 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) integer(pInt) :: & i, j, k, e - F => x_scal(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) - F_tau => x_scal(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) - residual_F => f_scal(1:3,1:3,1,& - X_RANGE, Y_RANGE, Z_RANGE) - residual_F_tau => f_scal(1:3,1:3,2,& - X_RANGE, Y_RANGE, Z_RANGE) + F => FandF_tau(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_tau => FandF_tau(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => residuum(1:3,1:3,1,& + X_RANGE, Y_RANGE, Z_RANGE) + residual_F_tau => residuum(1:3,1:3,2,& + X_RANGE, Y_RANGE, Z_RANGE) F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment !-------------------------------------------------------------------------------------------------- @@ -494,8 +497,8 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((1.0_pReal-mask_stress) * math_mul3333xx33(C_scale,F_aim-F_av) + & - mask_stress * (P_av-params%stress_BC))) ! mask = 0.0 for no bc + err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim-F_av) + & + params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc !-------------------------------------------------------------------------------------------------- ! error calculation diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index b209ab2ea..f9902ae41 100644 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -96,11 +96,10 @@ module spectral_utilities integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) end type tLoadCase - type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase including mask - real(pReal), dimension(3,3) :: stress_BC, rotation_BC + type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase + real(pReal), dimension(3,3) :: stress_mask, stress_BC, rotation_BC real(pReal) :: timeinc real(pReal) :: timeincOld - real(pReal) :: density end type tSolutionParams type, public :: phaseFieldDataBin !< set of parameters defining a phase field