diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index 8b1f0b7d8..3dbdad716 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -63,6 +63,7 @@ subroutine basic_init(temperature) res, & wgt, & geomdim, & + scaledDim, & mesh_ipCoordinates, & mesh_NcpElems, & mesh_deformedCoordsFFT @@ -82,6 +83,7 @@ subroutine basic_init(temperature) write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' write(6,'(a)') ' $Id$' #include "compilation_info.f90" + write(6,'(a,3(f12.5)/)') ' scaledDim x y z:', scaledDim !-------------------------------------------------------------------------------------------------- ! allocate global fields diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 4eb8a42ab..2283032a8 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -92,6 +92,7 @@ subroutine basicPETSc_init(temperature) res, & wgt, & geomdim, & + scaledDim, & mesh_NcpElems, & mesh_ipCoordinates, & mesh_deformedCoordsFFT @@ -117,6 +118,7 @@ subroutine basicPETSc_init(temperature) write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' write(6,'(a)') ' $Id: DAMASK_spectral_SolverBasicPETSC.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' #include "compilation_info.f90" + write(6,'(a,3(f12.5)/)') ' scaledDim x y z:', scaledDim !-------------------------------------------------------------------------------------------------- ! allocate global fields diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 2c83a34b6..110717964 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -123,7 +123,7 @@ subroutine utilities_init() use mesh, only: & res, & res1_red, & - virt_dim + scaledDim use math ! must use the whole module for use of FFTW implicit none @@ -230,7 +230,7 @@ subroutine utilities_init() 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 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)/virt_dim ! if divergence_correction is set, frequencies are calculated on unit length + xi(1:3,i,j,k) = real(k_s, pReal)/scaledDim ! if divergence_correction is set, frequencies are calculated on unit length enddo; enddo; enddo if(memory_efficient) then ! allocate just single fourth order tensor @@ -307,7 +307,7 @@ end subroutine utilities_updateGamma subroutine utilities_FFTforward(row,column) use math use mesh, only : & - virt_dim, & + scaledDim, & res, & res1_red @@ -325,7 +325,7 @@ subroutine utilities_FFTforward(row,column) !-------------------------------------------------------------------------------------------------- ! call function to calculate divergence from math (for post processing) to check results if (debugDivergence) & - divergence_post = math_divergenceFFT(virt_dim,field_real(1:res(1),1:res(2),1:res(3),1:3,1:3)) ! some elements are padded + divergence_post = math_divergenceFFT(scaledDim,field_real(1:res(1),1:res(2),1:res(3),1:3,1:3)) ! some elements are padded !-------------------------------------------------------------------------------------------------- ! doing the FFT diff --git a/code/mesh.f90 b/code/mesh.f90 index 80631cb27..a7ac6e3e7 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -34,7 +34,7 @@ module mesh implicit none private - integer(pInt), public :: & + integer(pInt), public, protected :: & mesh_NcpElems, & !< total number of CP elements in mesh mesh_NelemSets, & mesh_maxNelemInSet, & @@ -46,27 +46,31 @@ module mesh mesh_maxNsharedElems, & !< max number of CP elements sharing a node mesh_maxNsubNodes - integer(pInt), dimension(:,:), allocatable, public :: & + integer(pInt), dimension(:,:), allocatable, public, protected :: & mesh_element, & !< FEid, type(internal representation), material, texture, node indices mesh_sharedElem, & !< entryCount and list of elements containing node mesh_nodeTwins !< node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions) - integer(pInt), dimension(:,:,:,:), allocatable, public :: & + integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: & mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] real(pReal), dimension(:,:), allocatable, public :: & - mesh_ipVolume, & !< volume associated with IP (initially!) - mesh_node0, & !< node x,y,z coordinates (initially!) mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) - real(pReal), dimension(:,:,:), allocatable, public :: & - mesh_ipCoordinates, & !< IP x,y,z coordinates (after deformation!) - mesh_ipArea !< area of interface to neighboring IP (initially!) + real(pReal), dimension(:,:), allocatable, public, protected :: & + mesh_ipVolume, & !< volume associated with IP (initially!) + mesh_node0 !< node x,y,z coordinates (initially!) + + real(pReal), dimension(:,:,:), allocatable, public, protected :: & + mesh_ipArea !< area of interface to neighboring IP (initially!) - real(pReal),dimension(:,:,:,:), allocatable, public :: & + real(pReal), dimension(:,:,:), allocatable, public :: & + mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) + + real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) - logical, dimension(3), public :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes) + logical, dimension(3), public, protected :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes) integer(pInt), private :: & mesh_Nelems !< total number of elements in mesh @@ -98,10 +102,16 @@ module mesh #ifdef Spectral include 'fftw3.f03' - real(pReal), dimension(3), public :: geomdim, virt_dim ! physical dimension of volume element per direction - integer(pInt), dimension(3), public :: res - real(pReal), public :: wgt - integer(pInt), public :: res1_red, homog + 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 integer(pInt), private :: i #endif @@ -435,12 +445,20 @@ subroutine mesh_init(ip,element) wgt = 1.0/real(res(1)*res(2)*res(3),pReal) geomdim = mesh_spectral_getDimension(fileUnit) homog = mesh_spectral_getHomogenization(fileUnit) - if (divergence_correction) then + +!-------------------------------------------------------------------------------------------------- +! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso- +! lution-independent divergence + if (divergence_correction == 1_pInt) then do i = 1_pInt, 3_pInt - if (i /= minloc(geomdim,1) .and. i /= maxloc(geomdim,1)) virt_dim = geomdim/geomdim(i) + if (i/=minloc(geomdim,1) .and. i/=maxloc(geomdim,1)) scaledDim=geomdim/geomdim(i) + enddo + elseif (divergence_correction == 2_pInt) then + do i = 1_pInt, 3_pInt + if (i/=minloc(geomdim/res,1) .and. i/=maxloc(geomdim/res,1)) scaledDim=geomdim/geomdim(i)*res(i) enddo else - virt_dim = geomdim + scaledDim = geomdim endif write(6,'(a,3(i12 ))') ' resolution a b c:', res write(6,'(a,3(f12.5))') ' dimension x y z:', geomdim diff --git a/code/numerics.f90 b/code/numerics.f90 index abe20fecb..7ee43264e 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -104,10 +104,10 @@ module numerics itmax = 20_pInt, & !< maximum number of iterations itmin = 2_pInt, & !< minimum number of iterations maxCutBack = 3_pInt, & !< max number of cut backs - regridMode = 0_pInt !< 0: no regrid; 1: regrid if DAMASK doesn't converge; 2: regrid if DAMASK or BVP Solver doesn't converge + regridMode = 0_pInt, & !< 0: no regrid; 1: regrid if DAMASK doesn't converge; 2: regrid if DAMASK or BVP Solver doesn't converge + divergence_correction = 0_pInt !< correct divergence calculation in fourier space 0: no correction, 1: dimension scaled to 1, 2: dimension scaled to Npoints logical, protected , public :: & memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate - divergence_correction = .false., & !< correct divergence calculation in fourier space, Default .false.: no correction update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness #endif @@ -285,7 +285,7 @@ subroutine numerics_init case ('rotation_tol') rotation_tol = IO_floatValue(line,positions,2_pInt) case ('divergence_correction') - divergence_correction = IO_intValue(line,positions,2_pInt) > 0_pInt + divergence_correction = IO_intValue(line,positions,2_pInt) case ('update_gamma') update_gamma = IO_intValue(line,positions,2_pInt) > 0_pInt #ifdef PETSc @@ -409,7 +409,7 @@ subroutine numerics_init write(6,'(a24,1x,a)') ' myfilter: ',trim(myfilter) write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag write(6,'(a24,1x,es8.1)') ' rotation_tol: ',rotation_tol - write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction + write(6,'(a24,1x,i8)') ' divergence_correction: ',divergence_correction write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma #ifdef PETSc write(6,'(a24,1x,es8.1)') ' err_f_tol: ',err_f_tol @@ -470,6 +470,8 @@ subroutine numerics_init if (err_stress_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolabs') if (itmax <= 1.0_pInt) call IO_error(301_pInt,ext_msg='itmax') if (itmin > itmax .or. itmin < 1_pInt) call IO_error(301_pInt,ext_msg='itmin') + if (divergence_correction < 0_pInt .or. & + divergence_correction > 2_pInt) call IO_error(301_pInt,ext_msg='divergence_correction') if (maxCutBack <= 1.0_pInt) call IO_error(301_pInt,ext_msg='maxCutBack') if (update_gamma .and. & .not. memory_efficient) call IO_error(error_ID = 847_pInt)