223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

divergence_correction for basic solver variants has now 3 possibilities: 

0: uncorrected, slope per sidelength (physical dimension) e = res/dim
1: corrected by sidelength, slope per unitlength          e = res/1
2: corrected such that distance between FPs               e = 1

alway regarding the medium length of x,y,z direction
This commit is contained in:
Martin Diehl 2013-01-16 10:40:53 +00:00
parent 20bc97b7eb
commit c6a79d2b3d
5 changed files with 49 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
code/DAMASK_spectral_solverBasic.f90
View File

@ -63,6 +63,7 @@ subroutine basic_init(temperature)
res, & res, &
wgt, & wgt, &
geomdim, & geomdim, &
scaledDim, &
mesh_ipCoordinates, & mesh_ipCoordinates, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_deformedCoordsFFT mesh_deformedCoordsFFT
@ -82,6 +83,7 @@ subroutine basic_init(temperature)
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a,3(f12.5)/)') ' scaledDim x y z:', scaledDim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields

2
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -92,6 +92,7 @@ subroutine basicPETSc_init(temperature)
res, & res, &
wgt, & wgt, &
geomdim, & geomdim, &
scaledDim, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_ipCoordinates, & mesh_ipCoordinates, &
mesh_deformedCoordsFFT mesh_deformedCoordsFFT
@ -117,6 +118,7 @@ subroutine basicPETSc_init(temperature)
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' 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 $' write(6,'(a)') ' $Id: DAMASK_spectral_SolverBasicPETSC.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a,3(f12.5)/)') ' scaledDim x y z:', scaledDim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields

8
code/DAMASK_spectral_utilities.f90
View File

@ -123,7 +123,7 @@ subroutine utilities_init()
use mesh, only: & use mesh, only: &
res, & res, &
res1_red, & res1_red, &
virt_dim scaledDim
use math ! must use the whole module for use of FFTW use math ! must use the whole module for use of FFTW
implicit none 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 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 do i = 1_pInt, res1_red
k_s(1) = i - 1_pInt ! symmetry, junst running from 0,1,...,N/2,N/2+1 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 enddo; enddo; enddo
if(memory_efficient) then ! allocate just single fourth order tensor if(memory_efficient) then ! allocate just single fourth order tensor
@ -307,7 +307,7 @@ end subroutine utilities_updateGamma
subroutine utilities_FFTforward(row,column) subroutine utilities_FFTforward(row,column)
use math use math
use mesh, only : & use mesh, only : &
virt_dim, & scaledDim, &
res, & res, &
res1_red res1_red
@ -325,7 +325,7 @@ subroutine utilities_FFTforward(row,column)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! call function to calculate divergence from math (for post processing) to check results ! call function to calculate divergence from math (for post processing) to check results
if (debugDivergence) & 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 ! doing the FFT

52
code/mesh.f90
View File

@ -34,7 +34,7 @@ module mesh
implicit none implicit none
private private
integer(pInt), public :: & integer(pInt), public, protected :: &
mesh_NcpElems, & !< total number of CP elements in mesh mesh_NcpElems, & !< total number of CP elements in mesh
mesh_NelemSets, & mesh_NelemSets, &
mesh_maxNelemInSet, & mesh_maxNelemInSet, &
@ -46,27 +46,31 @@ module mesh
mesh_maxNsharedElems, & !< max number of CP elements sharing a node mesh_maxNsharedElems, & !< max number of CP elements sharing a node
mesh_maxNsubNodes mesh_maxNsubNodes
integer(pInt), dimension(:,:), allocatable, public :: & integer(pInt), dimension(:,:), allocatable, public, protected :: &
mesh_element, & !< FEid, type(internal representation), material, texture, node indices mesh_element, & !< FEid, type(internal representation), material, texture, node indices
mesh_sharedElem, & !< entryCount and list of elements containing node 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) 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] mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
real(pReal), dimension(:,:), allocatable, public :: & 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!!!) mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:,:), allocatable, public :: & real(pReal), dimension(:,:), allocatable, public, protected :: &
mesh_ipCoordinates, & !< IP x,y,z coordinates (after deformation!) mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_ipArea !< area of interface to neighboring 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!) 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 :: & integer(pInt), private :: &
mesh_Nelems !< total number of elements in mesh mesh_Nelems !< total number of elements in mesh
@ -98,10 +102,16 @@ module mesh
#ifdef Spectral #ifdef Spectral
include 'fftw3.f03' include 'fftw3.f03'
real(pReal), dimension(3), public :: geomdim, virt_dim ! physical dimension of volume element per direction real(pReal), dimension(3), public, protected :: &
integer(pInt), dimension(3), public :: res geomdim, & !< physical dimension of volume element per direction
real(pReal), public :: wgt scaledDim !< scaled dimension of volume element, depending on selected divergence calculation
integer(pInt), public :: res1_red, homog 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 integer(pInt), private :: i
#endif #endif
@ -435,12 +445,20 @@ subroutine mesh_init(ip,element)
wgt = 1.0/real(res(1)*res(2)*res(3),pReal) wgt = 1.0/real(res(1)*res(2)*res(3),pReal)
geomdim = mesh_spectral_getDimension(fileUnit) geomdim = mesh_spectral_getDimension(fileUnit)
homog = mesh_spectral_getHomogenization(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 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 enddo
else else
virt_dim = geomdim scaledDim = geomdim
endif endif
write(6,'(a,3(i12 ))') ' resolution a b c:', res write(6,'(a,3(i12 ))') ' resolution a b c:', res
write(6,'(a,3(f12.5))') ' dimension x y z:', geomdim write(6,'(a,3(f12.5))') ' dimension x y z:', geomdim

10
code/numerics.f90
View File

@ -104,10 +104,10 @@ module numerics
itmax = 20_pInt, & !< maximum number of iterations itmax = 20_pInt, & !< maximum number of iterations
itmin = 2_pInt, & !< minimum number of iterations itmin = 2_pInt, & !< minimum number of iterations
maxCutBack = 3_pInt, & !< max number of cut backs 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 :: & logical, protected , public :: &
memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate 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 update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness
#endif #endif
@ -285,7 +285,7 @@ subroutine numerics_init
case ('rotation_tol') case ('rotation_tol')
rotation_tol = IO_floatValue(line,positions,2_pInt) rotation_tol = IO_floatValue(line,positions,2_pInt)
case ('divergence_correction') case ('divergence_correction')
divergence_correction = IO_intValue(line,positions,2_pInt) > 0_pInt divergence_correction = IO_intValue(line,positions,2_pInt)
case ('update_gamma') case ('update_gamma')
update_gamma = IO_intValue(line,positions,2_pInt) > 0_pInt update_gamma = IO_intValue(line,positions,2_pInt) > 0_pInt
#ifdef PETSc #ifdef PETSc
@ -409,7 +409,7 @@ subroutine numerics_init
write(6,'(a24,1x,a)') ' myfilter: ',trim(myfilter) write(6,'(a24,1x,a)') ' myfilter: ',trim(myfilter)
write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag 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,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 write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma
#ifdef PETSc #ifdef PETSc
write(6,'(a24,1x,es8.1)') ' err_f_tol: ',err_f_tol 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 (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 (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 (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 (maxCutBack <= 1.0_pInt) call IO_error(301_pInt,ext_msg='maxCutBack')
if (update_gamma .and. & if (update_gamma .and. &
.not. memory_efficient) call IO_error(error_ID = 847_pInt) .not. memory_efficient) call IO_error(error_ID = 847_pInt)