2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
2012-10-24 17:01:40 +05:30
! $Id$
2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Utilities used by the different spectral solver variants
!--------------------------------------------------------------------------------------------------
2012-10-24 17:01:40 +05:30
module DAMASK_spectral_utilities
use , intrinsic :: iso_c_binding
2012-07-24 22:37:10 +05:30
use prec , only : &
pReal , &
pInt
2015-06-03 23:00:31 +05:30
use math , only : &
2015-09-11 18:35:46 +05:30
math_I3
2012-10-24 17:01:40 +05:30
2012-07-19 22:54:56 +05:30
implicit none
2013-03-28 16:07:00 +05:30
private
2012-12-15 23:37:49 +05:30
#ifdef PETSc
2015-08-04 20:34:53 +05:30
#include <petsc/finclude/petscsys.h>
2012-12-15 23:37:49 +05:30
#endif
2015-03-25 21:36:19 +05:30
include 'fftw3-mpi.f03'
2015-10-15 00:06:19 +05:30
2015-03-25 21:36:19 +05:30
logical , public :: cutBack = . false . !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
2013-11-14 00:51:35 +05:30
integer ( pInt ) , public , parameter :: maxPhaseFields = 2_pInt
2015-06-03 23:00:31 +05:30
integer ( pInt ) , public :: nActiveFields = 0_pInt
2015-10-15 00:06:19 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
! field labels information
enum , bind ( c )
enumerator :: FIELD_UNDEFINED_ID , &
FIELD_MECH_ID , &
FIELD_THERMAL_ID , &
FIELD_DAMAGE_ID , &
FIELD_VACANCYDIFFUSION_ID
end enum
2015-10-15 00:06:19 +05:30
2013-05-08 21:22:29 +05:30
!--------------------------------------------------------------------------------------------------
! grid related information information
real ( pReal ) , public :: wgt !< weighting factor 1/Nelems
2015-10-15 00:06:19 +05:30
2012-07-23 15:42:31 +05:30
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
2015-03-25 21:36:19 +05:30
integer ( pInt ) , public :: grid1Red !< grid(1)/2
2015-09-11 14:22:03 +05:30
real ( C_DOUBLE ) , public , dimension ( : , : , : , : , : ) , pointer :: tensorField_real !< real representation (some stress or deformation) of field_fourier
complex ( C_DOUBLE_COMPLEX ) , public , dimension ( : , : , : , : , : ) , pointer :: tensorField_fourier !< field on which the Fourier transform operates
real ( C_DOUBLE ) , public , dimension ( : , : , : , : ) , pointer :: vectorField_real !< vector field real representation for fftw
complex ( C_DOUBLE_COMPLEX ) , public , dimension ( : , : , : , : ) , pointer :: vectorField_fourier !< vector field fourier representation for fftw
real ( C_DOUBLE ) , public , dimension ( : , : , : ) , pointer :: scalarField_real !< scalar field real representation for fftw
complex ( C_DOUBLE_COMPLEX ) , public , dimension ( : , : , : ) , pointer :: scalarField_fourier !< scalar field fourier representation for fftw
2015-12-14 23:42:09 +05:30
complex ( pReal ) , private , dimension ( : , : , : , : , : , : , : ) , allocatable :: gamma_hat !< gamma operator (field) for spectral method
complex ( pReal ) , private , dimension ( : , : , : , : ) , allocatable :: xi1st !< wave vector field for first derivatives
complex ( pReal ) , private , dimension ( : , : , : , : ) , allocatable :: xi2nd !< wave vector field for second derivatives
2015-06-03 23:00:31 +05:30
real ( pReal ) , private , dimension ( 3 , 3 , 3 , 3 ) :: C_ref !< mechanic reference stiffness
real ( pReal ) , protected , public , dimension ( 3 ) :: scaledGeomSize !< scaled geometry size for calculation of divergence (Basic, Basic PETSc)
2015-10-15 00:06:19 +05:30
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! plans for FFTW
2013-04-10 15:49:16 +05:30
type ( C_PTR ) , private :: &
2015-09-11 14:22:03 +05:30
planTensorForth , & !< FFTW MPI plan P(x) to P(k)
planTensorBack , & !< FFTW MPI plan F(k) to F(x)
planVectorForth , & !< FFTW MPI plan v(x) to v(k)
planVectorBack , & !< FFTW MPI plan v(k) to v(x)
planScalarForth , & !< FFTW MPI plan s(x) to s(k)
planScalarBack !< FFTW MPI plan s(k) to s(x)
2012-11-12 19:44:39 +05:30
2012-07-23 15:42:31 +05:30
!--------------------------------------------------------------------------------------------------
2012-10-24 17:01:40 +05:30
! variables controlling debugging
2013-05-08 21:22:29 +05:30
logical , private :: &
2012-10-24 17:01:40 +05:30
debugGeneral , & !< general debugging of spectral solver
2012-12-15 23:37:49 +05:30
debugRotation , & !< also printing out results in lab frame
debugPETSc !< use some in debug defined options for more verbose PETSc solution
2012-07-23 15:42:31 +05:30
2012-07-25 19:31:39 +05:30
!--------------------------------------------------------------------------------------------------
2012-08-03 14:55:48 +05:30
! derived types
2013-04-10 15:49:16 +05:30
type , public :: tSolutionState !< return type of solution from spectral solver variants
2015-10-15 00:06:19 +05:30
logical :: converged = . true .
logical :: regrid = . false .
logical :: stagConverged = . true .
logical :: termIll = . false .
2012-10-24 17:01:40 +05:30
integer ( pInt ) :: iterationsNeeded = 0_pInt
2012-10-02 20:56:56 +05:30
end type tSolutionState
2012-08-03 14:55:48 +05:30
2013-04-10 15:49:16 +05:30
type , public :: tBoundaryCondition !< set of parameters defining a boundary condition
2012-10-24 17:01:40 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: values = 0.0_pReal
real ( pReal ) , dimension ( 3 , 3 ) :: maskFloat = 0.0_pReal
2012-08-03 14:55:48 +05:30
logical , dimension ( 3 , 3 ) :: maskLogical = . false .
2012-10-24 17:01:40 +05:30
character ( len = 64 ) :: myType = 'None'
2012-10-02 20:56:56 +05:30
end type tBoundaryCondition
2013-11-14 00:51:35 +05:30
2015-06-03 23:00:31 +05:30
type , public :: tLoadCase
real ( pReal ) , dimension ( 3 , 3 ) :: rotation = math_I3 !< rotation of BC
type ( tBoundaryCondition ) :: P , & !< stress BC
deformation !< deformation BC (Fdot or L)
2015-07-24 20:27:29 +05:30
real ( pReal ) :: time = 0.0_pReal !< length of increment
2015-06-03 23:00:31 +05:30
integer ( pInt ) :: incs = 0_pInt , & !< number of increments
outputfrequency = 1_pInt , & !< frequency of result writes
restartfrequency = 0_pInt , & !< frequency of restart writes
logscale = 0_pInt !< linear/logarithmic time inc flag
2015-10-15 00:06:19 +05:30
logical :: followFormerTrajectory = . true . !< follow trajectory of former loadcase
2015-06-03 23:00:31 +05:30
integer ( kind ( FIELD_UNDEFINED_ID ) ) , allocatable :: ID ( : )
end type tLoadCase
2015-10-15 00:06:19 +05:30
2015-06-03 23:00:31 +05:30
type , public :: tSolutionParams !< @todo use here the type definition for a full loadcase including mask
real ( pReal ) , dimension ( 3 , 3 ) :: P_BC , rotation_BC
real ( pReal ) :: timeinc
real ( pReal ) :: timeincOld
real ( pReal ) :: density
end type tSolutionParams
2015-10-15 00:06:19 +05:30
2015-06-03 23:00:31 +05:30
type ( tSolutionParams ) , private :: params
2015-10-15 00:06:19 +05:30
2015-03-25 21:36:19 +05:30
type , public :: phaseFieldDataBin !< set of parameters defining a phase field
real ( pReal ) :: diffusion = 0.0_pReal , & !< thermal conductivity
mobility = 0.0_pReal , & !< thermal mobility
2015-10-15 00:06:19 +05:30
phaseField0 = 0.0_pReal !< homogeneous damage field starting condition
2013-11-14 00:51:35 +05:30
logical :: active = . false .
character ( len = 64 ) :: label = ''
end type phaseFieldDataBin
2015-09-11 18:35:46 +05:30
2015-10-15 00:06:19 +05:30
enum , bind ( c )
2015-12-14 23:42:09 +05:30
enumerator :: DERIVATIVE_CONTINUOUS_ID , &
DERIVATIVE_CENTRAL_DIFF_ID , &
DERIVATIVE_FWBW_DIFF_ID
2015-09-11 18:35:46 +05:30
end enum
2015-12-14 23:42:09 +05:30
integer ( kind ( DERIVATIVE_CONTINUOUS_ID ) ) :: &
spectral_derivative_ID
2015-10-15 00:06:19 +05:30
2012-10-24 17:01:40 +05:30
public :: &
utilities_init , &
utilities_updateGamma , &
2015-06-03 23:00:31 +05:30
utilities_FFTtensorForward , &
utilities_FFTtensorBackward , &
utilities_FFTvectorForward , &
utilities_FFTvectorBackward , &
utilities_FFTscalarForward , &
utilities_FFTscalarBackward , &
utilities_fourierGammaConvolution , &
utilities_fourierGreenConvolution , &
2012-10-24 17:01:40 +05:30
utilities_divergenceRMS , &
2013-03-07 01:04:30 +05:30
utilities_curlRMS , &
2015-06-03 23:00:31 +05:30
utilities_fourierScalarGradient , &
utilities_fourierVectorDivergence , &
2015-12-14 23:42:09 +05:30
utilities_fourierVectorGradient , &
utilities_fourierTensorDivergence , &
2012-10-24 17:01:40 +05:30
utilities_maskedCompliance , &
utilities_constitutiveResponse , &
utilities_calculateRate , &
utilities_forwardField , &
2015-03-13 03:58:33 +05:30
utilities_destroy , &
2015-06-03 23:00:31 +05:30
utilities_updateIPcoords , &
FIELD_UNDEFINED_ID , &
FIELD_MECH_ID , &
FIELD_THERMAL_ID , &
FIELD_DAMAGE_ID
2012-10-24 17:01:40 +05:30
private :: &
2015-12-14 23:42:09 +05:30
utilities_getFreqDerivative
2012-08-03 14:55:48 +05:30
2015-10-15 00:06:19 +05:30
contains
2012-07-19 22:54:56 +05:30
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
2012-11-12 19:44:39 +05:30
!> @brief allocates all neccessary fields, sets debug flags, create plans for FFTW
2015-10-15 00:06:19 +05:30
!> @details Sets the debug levels for general, divergence, restart and FFTW from the biwise coding
2012-08-09 18:34:56 +05:30
!> provided by the debug module to logicals.
!> Allocates all fields used by FFTW and create the corresponding plans depending on the debug
!> level chosen.
!> Initializes FFTW.
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
2012-10-24 17:01:40 +05:30
subroutine utilities_init ( )
2012-08-27 13:34:47 +05:30
use , intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
2012-10-24 17:01:40 +05:30
use IO , only : &
2013-01-10 19:03:43 +05:30
IO_error , &
2013-02-25 22:04:59 +05:30
IO_warning , &
2013-05-08 21:22:29 +05:30
IO_timeStamp , &
IO_open_file
2015-06-03 23:00:31 +05:30
use numerics , only : &
2015-12-14 23:42:09 +05:30
spectral_derivative , &
2012-07-20 21:03:13 +05:30
fftw_planner_flag , &
2012-07-30 19:36:22 +05:30
fftw_timelimit , &
2012-12-15 23:37:49 +05:30
memory_efficient , &
2015-06-09 18:58:50 +05:30
petsc_defaultOptions , &
2013-05-08 21:22:29 +05:30
petsc_options , &
2015-03-25 21:36:19 +05:30
divergence_correction , &
worldrank
2012-07-20 21:03:13 +05:30
use debug , only : &
debug_level , &
2013-05-16 14:25:10 +05:30
debug_SPECTRAL , &
debug_LEVELBASIC , &
debug_SPECTRALDIVERGENCE , &
debug_SPECTRALFFTW , &
debug_SPECTRALPETSC , &
debug_SPECTRALROTATION
2012-12-15 23:37:49 +05:30
#ifdef PETSc
2013-01-02 22:32:12 +05:30
use debug , only : &
2013-05-16 14:25:10 +05:30
PETSCDEBUG
2012-12-15 23:37:49 +05:30
#endif
2015-10-15 00:06:19 +05:30
use math
2013-05-08 21:22:29 +05:30
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid , &
grid3 , &
grid3Offset , &
geomSize
2012-12-15 23:37:49 +05:30
2012-07-24 22:37:10 +05:30
implicit none
2012-12-15 23:37:49 +05:30
#ifdef PETSc
2013-02-13 23:24:56 +05:30
external :: &
PETScOptionsClear , &
PETScOptionsInsertString , &
MPI_Abort
2012-12-15 23:37:49 +05:30
PetscErrorCode :: ierr
2015-10-15 00:06:19 +05:30
#endif
2012-10-24 17:01:40 +05:30
integer ( pInt ) :: i , j , k
2012-07-24 22:37:10 +05:30
integer ( pInt ) , dimension ( 3 ) :: k_s
2012-10-24 17:01:40 +05:30
type ( C_PTR ) :: &
2015-09-11 18:35:46 +05:30
tensorField , & !< field containing data for FFTW in real and fourier space (in place)
vectorField , & !< field containing data for FFTW in real space when debugging FFTW (no in place)
scalarField !< field containing data for FFTW in real space when debugging FFTW (no in place)
2015-09-24 23:08:49 +05:30
integer ( C_INTPTR_T ) , dimension ( 3 ) :: gridFFTW
integer ( C_INTPTR_T ) :: alloc_local , local_K , local_K_offset
2015-03-31 02:19:17 +05:30
integer ( C_INTPTR_T ) , parameter :: &
scalarSize = 1_C_INTPTR_T , &
vecSize = 3_C_INTPTR_T , &
tensorSize = 9_C_INTPTR_T
2015-10-15 00:06:19 +05:30
2015-03-25 21:36:19 +05:30
mainProcess : if ( worldrank == 0 ) then
write ( 6 , '(/,a)' ) ' <<<+- DAMASK_spectral_utilities init -+>>>'
write ( 6 , '(a)' ) ' $Id$'
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
2012-07-19 22:54:56 +05:30
#include "compilation_info.f90"
2015-03-25 21:36:19 +05:30
endif mainProcess
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
2012-07-23 15:42:31 +05:30
! set debugging parameters
2013-05-16 14:25:10 +05:30
debugGeneral = iand ( debug_level ( debug_SPECTRAL ) , debug_LEVELBASIC ) / = 0
debugRotation = iand ( debug_level ( debug_SPECTRAL ) , debug_SPECTRALROTATION ) / = 0
debugPETSc = iand ( debug_level ( debug_SPECTRAL ) , debug_SPECTRALPETSC ) / = 0
2015-03-31 02:19:17 +05:30
2015-03-25 21:36:19 +05:30
if ( debugPETSc . and . worldrank == 0_pInt ) write ( 6 , '(3(/,a),/)' ) &
2013-07-08 21:18:13 +05:30
' Initializing PETSc with debug options: ' , &
trim ( PETScDebug ) , &
' add more using the PETSc_Options keyword in numerics.config '
2013-01-10 03:49:32 +05:30
flush ( 6 )
call PetscOptionsClear ( ierr ) ; CHKERRQ ( ierr )
2013-05-16 14:25:10 +05:30
if ( debugPETSc ) call PetscOptionsInsertString ( trim ( PETSCDEBUG ) , ierr ) ; CHKERRQ ( ierr )
2015-06-09 18:58:50 +05:30
call PetscOptionsInsertString ( trim ( petsc_defaultOptions ) , ierr ) ; CHKERRQ ( ierr )
2013-01-10 03:49:32 +05:30
call PetscOptionsInsertString ( trim ( petsc_options ) , ierr ) ; CHKERRQ ( ierr )
2013-03-07 01:04:30 +05:30
2015-09-11 14:22:03 +05:30
grid1Red = grid ( 1 ) / 2_pInt + 1_pInt
wgt = 1.0 / real ( product ( grid ) , pReal )
2013-05-08 21:22:29 +05:30
2015-03-25 21:36:19 +05:30
if ( worldrank == 0 ) then
2015-09-11 14:22:03 +05:30
write ( 6 , '(a,3(i12 ))' ) ' grid a b c: ' , grid
write ( 6 , '(a,3(es12.5))' ) ' size x y z: ' , geomSize
2015-03-25 21:36:19 +05:30
endif
2015-10-15 00:06:19 +05:30
2015-12-14 23:42:09 +05:30
select case ( spectral_derivative )
case ( 'continuous' ) ! default, no weighting
spectral_derivative_ID = DERIVATIVE_CONTINUOUS_ID
case ( 'central_difference' ) ! cosine curve with 1 for avg and zero for highest freq
spectral_derivative_ID = DERIVATIVE_CENTRAL_DIFF_ID
case ( 'fwbw_difference' ) ! gradient, might need grid scaling as for cosine filter
spectral_derivative_ID = DERIVATIVE_FWBW_DIFF_ID
case default
call IO_error ( 892_pInt , ext_msg = trim ( spectral_derivative ) )
end select
2013-05-08 21:22:29 +05:30
!--------------------------------------------------------------------------------------------------
2015-09-11 14:22:03 +05:30
! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and
! resolution-independent divergence
2013-05-08 21:22:29 +05:30
if ( divergence_correction == 1_pInt ) then
do j = 1_pInt , 3_pInt
2015-09-11 14:22:03 +05:30
if ( j / = minloc ( geomSize , 1 ) . and . j / = maxloc ( geomSize , 1 ) ) &
scaledGeomSize = geomSize / geomSize ( j )
2013-05-08 21:22:29 +05:30
enddo
elseif ( divergence_correction == 2_pInt ) then
do j = 1_pInt , 3_pInt
2015-09-11 14:22:03 +05:30
if ( j / = minloc ( geomSize / grid , 1 ) . and . j / = maxloc ( geomSize / grid , 1 ) ) &
scaledGeomSize = geomSize / geomSize ( j ) * grid ( j )
2013-05-08 21:22:29 +05:30
enddo
else
2015-09-11 14:22:03 +05:30
scaledGeomSize = geomSize
2013-05-08 21:22:29 +05:30
endif
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
2015-03-25 21:36:19 +05:30
! MPI allocation
2015-09-11 14:22:03 +05:30
gridFFTW = int ( grid , C_INTPTR_T )
2015-03-25 21:36:19 +05:30
alloc_local = fftw_mpi_local_size_3d ( gridFFTW ( 3 ) , gridFFTW ( 2 ) , gridFFTW ( 1 ) / 2 + 1 , &
MPI_COMM_WORLD , local_K , local_K_offset )
2015-12-14 23:42:09 +05:30
allocate ( xi1st ( 3 , grid1Red , grid ( 2 ) , grid3 ) , source = cmplx ( 0.0_pReal , 0.0_pReal , pReal ) ) ! frequencies, only half the size for first dimension
allocate ( xi2nd ( 3 , grid1Red , grid ( 2 ) , grid3 ) , source = cmplx ( 0.0_pReal , 0.0_pReal , pReal ) ) ! frequencies, only half the size for first dimension
2015-09-11 14:22:03 +05:30
tensorField = fftw_alloc_complex ( tensorSize * alloc_local )
call c_f_pointer ( tensorField , tensorField_real , [ 3_C_INTPTR_T , 3_C_INTPTR_T , &
2015-06-03 23:00:31 +05:30
2_C_INTPTR_T * ( gridFFTW ( 1 ) / 2_C_INTPTR_T + 1_C_INTPTR_T ) , gridFFTW ( 2 ) , local_K ] ) ! place a pointer for a real tensor representation
2015-09-11 14:22:03 +05:30
call c_f_pointer ( tensorField , tensorField_fourier , [ 3_C_INTPTR_T , 3_C_INTPTR_T , &
2015-06-03 23:00:31 +05:30
gridFFTW ( 1 ) / 2_C_INTPTR_T + 1_C_INTPTR_T , gridFFTW ( 2 ) , local_K ] ) ! place a pointer for a fourier tensor representation
2015-09-11 14:22:03 +05:30
vectorField = fftw_alloc_complex ( vecSize * alloc_local )
call c_f_pointer ( vectorField , vectorField_real , [ 3_C_INTPTR_T , &
2015-06-03 23:00:31 +05:30
2_C_INTPTR_T * ( gridFFTW ( 1 ) / 2_C_INTPTR_T + 1_C_INTPTR_T ) , gridFFTW ( 2 ) , local_K ] ) ! place a pointer for a real vector representation
2015-09-11 14:22:03 +05:30
call c_f_pointer ( vectorField , vectorField_fourier , [ 3_C_INTPTR_T , &
2015-06-03 23:00:31 +05:30
gridFFTW ( 1 ) / 2_C_INTPTR_T + 1_C_INTPTR_T , gridFFTW ( 2 ) , local_K ] ) ! place a pointer for a fourier vector representation
2015-09-11 14:22:03 +05:30
scalarField = fftw_alloc_complex ( scalarSize * alloc_local ) ! allocate data for real representation (no in place transform)
call c_f_pointer ( scalarField , scalarField_real , &
2015-06-03 23:00:31 +05:30
[ 2_C_INTPTR_T * ( gridFFTW ( 1 ) / 2_C_INTPTR_T + 1 ) , gridFFTW ( 2 ) , local_K ] ) ! place a pointer for a real scalar representation
2015-09-11 14:22:03 +05:30
call c_f_pointer ( scalarField , scalarField_fourier , &
2015-06-03 23:00:31 +05:30
[ gridFFTW ( 1 ) / 2_C_INTPTR_T + 1 , gridFFTW ( 2 ) , local_K ] ) ! place a pointer for a fourier scarlar representation
2015-10-15 00:06:19 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
! tensor MPI fftw plans
2015-09-11 14:22:03 +05:30
planTensorForth = fftw_mpi_plan_many_dft_r2c ( 3 , [ gridFFTW ( 3 ) , gridFFTW ( 2 ) , gridFFTW ( 1 ) ] , & ! dimension, logical length in each dimension in reversed order
2015-06-03 23:00:31 +05:30
tensorSize , FFTW_MPI_DEFAULT_BLOCK , FFTW_MPI_DEFAULT_BLOCK , & ! no. of transforms, default iblock and oblock
2015-09-11 14:22:03 +05:30
tensorField_real , tensorField_fourier , & ! input data, output data
2015-06-03 23:00:31 +05:30
MPI_COMM_WORLD , fftw_planner_flag ) ! use all processors, planer precision
2015-07-09 19:08:21 +05:30
if ( . not . C_ASSOCIATED ( planTensorForth ) ) call IO_error ( 810 , ext_msg = 'planTensorForth' )
2015-09-11 14:22:03 +05:30
planTensorBack = fftw_mpi_plan_many_dft_c2r ( 3 , [ gridFFTW ( 3 ) , gridFFTW ( 2 ) , gridFFTW ( 1 ) ] , & ! dimension, logical length in each dimension in reversed order
2015-06-03 23:00:31 +05:30
tensorSize , FFTW_MPI_DEFAULT_BLOCK , FFTW_MPI_DEFAULT_BLOCK , & ! no. of transforms, default iblock and oblock
2015-09-11 14:22:03 +05:30
tensorField_fourier , tensorField_real , & ! input data, output data
2015-06-03 23:00:31 +05:30
MPI_COMM_WORLD , fftw_planner_flag ) ! all processors, planer precision
2015-07-09 19:08:21 +05:30
if ( . not . C_ASSOCIATED ( planTensorBack ) ) call IO_error ( 810 , ext_msg = 'planTensorBack' )
2015-10-15 00:06:19 +05:30
2015-03-13 03:58:33 +05:30
!--------------------------------------------------------------------------------------------------
2015-06-03 23:00:31 +05:30
! vector MPI fftw plans
2015-09-11 14:22:03 +05:30
planVectorForth = fftw_mpi_plan_many_dft_r2c ( 3 , [ gridFFTW ( 3 ) , gridFFTW ( 2 ) , gridFFTW ( 1 ) ] , & ! dimension, logical length in each dimension in reversed order
2015-06-03 23:00:31 +05:30
vecSize , FFTW_MPI_DEFAULT_BLOCK , FFTW_MPI_DEFAULT_BLOCK , & ! no. of transforms, default iblock and oblock
2015-09-11 14:22:03 +05:30
vectorField_real , vectorField_fourier , & ! input data, output data
2015-06-03 23:00:31 +05:30
MPI_COMM_WORLD , fftw_planner_flag ) ! use all processors, planer precision
2015-07-09 19:08:21 +05:30
if ( . not . C_ASSOCIATED ( planVectorForth ) ) call IO_error ( 810 , ext_msg = 'planVectorForth' )
2015-09-11 14:22:03 +05:30
planVectorBack = fftw_mpi_plan_many_dft_c2r ( 3 , [ gridFFTW ( 3 ) , gridFFTW ( 2 ) , gridFFTW ( 1 ) ] , & ! dimension, logical length in each dimension in reversed order
2015-06-03 23:00:31 +05:30
vecSize , FFTW_MPI_DEFAULT_BLOCK , FFTW_MPI_DEFAULT_BLOCK , & ! no. of transforms, default iblock and oblock
2015-09-11 14:22:03 +05:30
vectorField_fourier , vectorField_real , & ! input data, output data
2015-06-03 23:00:31 +05:30
MPI_COMM_WORLD , fftw_planner_flag ) ! all processors, planer precision
2015-07-09 19:08:21 +05:30
if ( . not . C_ASSOCIATED ( planVectorBack ) ) call IO_error ( 810 , ext_msg = 'planVectorBack' )
2015-10-15 00:06:19 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
2015-10-15 00:06:19 +05:30
! scalar MPI fftw plans
planScalarForth = fftw_mpi_plan_many_dft_r2c ( 3 , [ gridFFTW ( 3 ) , gridFFTW ( 2 ) , gridFFTW ( 1 ) ] , & ! dimension, logical length in each dimension in reversed order
2015-06-03 23:00:31 +05:30
scalarSize , FFTW_MPI_DEFAULT_BLOCK , FFTW_MPI_DEFAULT_BLOCK , & ! no. of transforms, default iblock and oblock
2015-09-11 14:22:03 +05:30
scalarField_real , scalarField_fourier , & ! input data, output data
2015-06-03 23:00:31 +05:30
MPI_COMM_WORLD , fftw_planner_flag ) ! use all processors, planer precision
2015-07-09 19:08:21 +05:30
if ( . not . C_ASSOCIATED ( planScalarForth ) ) call IO_error ( 810 , ext_msg = 'planScalarForth' )
planScalarBack = fftw_mpi_plan_many_dft_c2r ( 3 , [ gridFFTW ( 3 ) , gridFFTW ( 2 ) , gridFFTW ( 1 ) ] , & ! dimension, logical length in each dimension in reversed order, no. of transforms
2015-06-03 23:00:31 +05:30
scalarSize , FFTW_MPI_DEFAULT_BLOCK , FFTW_MPI_DEFAULT_BLOCK , & ! no. of transforms, default iblock and oblock
2015-09-11 14:22:03 +05:30
scalarField_fourier , scalarField_real , & ! input data, output data
2015-06-03 23:00:31 +05:30
MPI_COMM_WORLD , fftw_planner_flag ) ! use all processors, planer precision
2015-07-09 19:08:21 +05:30
if ( . not . C_ASSOCIATED ( planScalarBack ) ) call IO_error ( 810 , ext_msg = 'planScalarBack' )
2012-07-19 22:54:56 +05:30
2015-03-25 21:36:19 +05:30
!--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details)
if ( pReal / = C_DOUBLE . or . pInt / = C_INT ) call IO_error ( 0_pInt , ext_msg = 'Fortran to C' ) ! check for correct precision in C
call fftw_set_timelimit ( fftw_timelimit ) ! set timelimit for plan creation
2012-07-19 22:54:56 +05:30
2015-03-25 21:36:19 +05:30
if ( debugGeneral . and . worldrank == 0_pInt ) write ( 6 , '(/,a)' ) ' FFTW initialized'
flush ( 6 )
2015-10-15 00:06:19 +05:30
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
2012-07-23 15:42:31 +05:30
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
2015-09-11 14:22:03 +05:30
do k = grid3Offset + 1_pInt , grid3Offset + grid3
2012-07-23 15:42:31 +05:30
k_s ( 3 ) = k - 1_pInt
2015-09-11 14:22:03 +05:30
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 )
2012-07-23 15:42:31 +05:30
k_s ( 2 ) = j - 1_pInt
2015-09-11 14:22:03 +05:30
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
2013-05-08 21:22:29 +05:30
do i = 1_pInt , grid1Red
2012-10-24 17:01:40 +05:30
k_s ( 1 ) = i - 1_pInt ! symmetry, junst running from 0,1,...,N/2,N/2+1
2015-12-14 23:42:09 +05:30
xi2nd ( 1 : 3 , i , j , k - grid3Offset ) = utilities_getFreqDerivative ( k_s ) ! if divergence_correction is set, frequencies are calculated on unit length
where ( mod ( grid , 2 ) == 0 . and . [ i , j , k ] == grid / 2 + 1 . and . &
spectral_derivative_ID == DERIVATIVE_CONTINUOUS_ID ) ! for even grids, set the Nyquist Freq component to 0.0
xi1st ( 1 : 3 , i , j , k - grid3Offset ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2015-09-11 18:35:46 +05:30
elsewhere
xi1st ( 1 : 3 , i , j , k - grid3Offset ) = xi2nd ( 1 : 3 , i , j , k - grid3Offset )
endwhere
2012-07-23 15:42:31 +05:30
enddo ; enddo ; enddo
2015-10-15 00:06:19 +05:30
2012-07-19 22:54:56 +05:30
if ( memory_efficient ) then ! allocate just single fourth order tensor
2015-12-14 23:42:09 +05:30
allocate ( gamma_hat ( 3 , 3 , 3 , 3 , 1 , 1 , 1 ) , source = cmplx ( 0.0_pReal , 0.0_pReal , pReal ) )
2012-07-19 22:54:56 +05:30
else ! precalculation of gamma_hat field
2015-12-14 23:42:09 +05:30
allocate ( gamma_hat ( 3 , 3 , 3 , 3 , grid1Red , grid ( 2 ) , grid3 ) , source = cmplx ( 0.0_pReal , 0.0_pReal , pReal ) )
2012-07-19 22:54:56 +05:30
endif
2012-08-03 14:55:48 +05:30
2012-10-24 17:01:40 +05:30
end subroutine utilities_init
2012-07-26 19:28:47 +05:30
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief updates references stiffness and potentially precalculated gamma operator
2012-08-09 18:34:56 +05:30
!> @details Sets the current reference stiffness to the stiffness given as an argument.
!> If the gamma operator is precalculated, it is calculated with this stiffness.
!> In case of a on-the-fly calculation, only the reference stiffness is updated.
2012-10-24 17:01:40 +05:30
!> Also writes out the current reference stiffness for restart.
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
2012-10-24 17:01:40 +05:30
subroutine utilities_updateGamma ( C , saveReference )
use IO , only : &
2013-09-18 19:37:55 +05:30
IO_write_jobRealFile
2012-08-03 14:55:48 +05:30
use numerics , only : &
2015-03-25 21:36:19 +05:30
memory_efficient , &
worldrank
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid3Offset , &
grid3 , &
grid
2012-10-24 17:01:40 +05:30
use math , only : &
2015-12-14 23:42:09 +05:30
math_det33 , &
math_invert
2012-10-24 17:01:40 +05:30
2012-08-03 14:55:48 +05:30
implicit none
2012-10-24 17:01:40 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 , 3 , 3 ) :: C !< input stiffness to store as reference stiffness
logical , intent ( in ) :: saveReference !< save reference stiffness to file for restart
2015-12-14 23:42:09 +05:30
complex ( pReal ) , dimension ( 3 , 3 ) :: temp33_complex , xiDyad_cmplx
real ( pReal ) , dimension ( 6 , 6 ) :: matA , matInvA
2012-10-24 17:01:40 +05:30
integer ( pInt ) :: &
i , j , k , &
l , m , n , o
2015-12-14 23:42:09 +05:30
logical :: ierr
2015-10-15 00:06:19 +05:30
2012-08-03 14:55:48 +05:30
C_ref = C
2012-10-24 17:01:40 +05:30
if ( saveReference ) then
2015-03-25 21:36:19 +05:30
if ( worldrank == 0_pInt ) then
write ( 6 , '(/,a)' ) ' writing reference stiffness to file'
flush ( 6 )
call IO_write_jobRealFile ( 777 , 'C_ref' , size ( C_ref ) )
write ( 777 , rec = 1 ) C_ref
close ( 777 )
endif
2012-10-24 17:01:40 +05:30
endif
2015-10-15 00:06:19 +05:30
if ( . not . memory_efficient ) then
2015-09-11 14:22:03 +05:30
do k = grid3Offset + 1_pInt , grid3Offset + grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid1Red
2015-03-25 21:36:19 +05:30
if ( any ( [ i , j , k ] / = 1_pInt ) ) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
2012-08-03 14:55:48 +05:30
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
2015-12-14 23:42:09 +05:30
xiDyad_cmplx ( l , m ) = conjg ( - xi1st ( l , i , j , k - grid3Offset ) ) * xi1st ( m , i , j , k - grid3Offset )
2012-08-03 14:55:48 +05:30
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
2015-12-14 23:42:09 +05:30
temp33_complex ( l , m ) = sum ( C_ref ( l , 1 : 3 , m , 1 : 3 ) * xiDyad_cmplx )
matA ( 1 : 3 , 1 : 3 ) = real ( temp33_complex ) ; matA ( 4 : 6 , 4 : 6 ) = real ( temp33_complex )
matA ( 1 : 3 , 4 : 6 ) = aimag ( temp33_complex ) ; matA ( 4 : 6 , 1 : 3 ) = - aimag ( temp33_complex )
if ( abs ( math_det33 ( matA ( 1 : 3 , 1 : 3 ) ) ) > 1e-16 ) then
call math_invert ( 6_pInt , matA , matInvA , ierr )
temp33_complex = cmplx ( matInvA ( 1 : 3 , 1 : 3 ) , matInvA ( 1 : 3 , 4 : 6 ) , pReal )
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt , n = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
gamma_hat ( l , m , n , o , i , j , k - grid3Offset ) = temp33_complex ( l , n ) * &
conjg ( - xi1st ( o , i , j , k - grid3Offset ) ) * xi1st ( m , i , j , k - grid3Offset )
else
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt , n = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
gamma_hat ( l , m , n , o , 1 , 1 , 1 ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
endif
2015-10-15 00:06:19 +05:30
endif
2012-08-03 14:55:48 +05:30
enddo ; enddo ; enddo
2015-10-15 00:06:19 +05:30
endif
2012-10-24 17:01:40 +05:30
end subroutine utilities_updateGamma
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
2015-09-11 18:35:46 +05:30
!> @brief forward FFT of data in field_real to field_fourier
!> @details Does an unweighted filtered FFT transform from real to complex
2012-08-03 14:55:48 +05:30
!--------------------------------------------------------------------------------------------------
2015-06-03 23:00:31 +05:30
subroutine utilities_FFTtensorForward ( )
2015-03-25 21:36:19 +05:30
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid3 , &
grid
2012-10-24 17:01:40 +05:30
implicit none
2013-03-07 01:04:30 +05:30
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
2015-09-11 18:35:46 +05:30
! doing the tensor FFT
2015-09-11 14:22:03 +05:30
call fftw_mpi_execute_dft_r2c ( planTensorForth , tensorField_real , tensorField_fourier )
2015-10-15 00:06:19 +05:30
2015-06-03 23:00:31 +05:30
end subroutine utilities_FFTtensorForward
2012-07-26 19:28:47 +05:30
2012-08-09 18:34:56 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief backward FFT of data in field_fourier to field_real
2015-09-11 14:22:03 +05:30
!> @details Does an weighted inverse FFT transform from complex to real
2012-08-09 18:34:56 +05:30
!--------------------------------------------------------------------------------------------------
2015-06-03 23:00:31 +05:30
subroutine utilities_FFTtensorBackward ( )
2012-10-24 17:01:40 +05:30
implicit none
2013-04-10 15:49:16 +05:30
2015-09-11 18:35:46 +05:30
call fftw_mpi_execute_dft_c2r ( planTensorBack , tensorField_fourier , tensorField_real )
2015-09-11 14:22:03 +05:30
tensorField_real = tensorField_real * wgt ! normalize the result by number of elements
2012-08-03 14:55:48 +05:30
2015-06-03 23:00:31 +05:30
end subroutine utilities_FFTtensorBackward
2012-07-26 19:28:47 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
2015-09-11 14:22:03 +05:30
!> @brief forward FFT of data in scalarField_real to scalarField_fourier
2015-09-11 18:35:46 +05:30
!> @details Does an unweighted filtered FFT transform from real to complex
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarForward ( )
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid3 , &
grid
2015-10-15 00:06:19 +05:30
2015-09-11 14:22:03 +05:30
implicit none
2015-10-15 00:06:19 +05:30
2015-09-11 14:22:03 +05:30
!--------------------------------------------------------------------------------------------------
2015-06-03 23:00:31 +05:30
! doing the scalar FFT
2015-09-11 14:22:03 +05:30
call fftw_mpi_execute_dft_r2c ( planScalarForth , scalarField_real , scalarField_fourier )
2015-06-03 23:00:31 +05:30
end subroutine utilities_FFTscalarForward
2015-09-24 23:08:49 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
2015-09-11 14:22:03 +05:30
!> @brief backward FFT of data in scalarField_fourier to scalarField_real
!> @details Does an weighted inverse FFT transform from complex to real
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarBackward ( )
2015-09-11 14:22:03 +05:30
implicit none
call fftw_mpi_execute_dft_c2r ( planScalarBack , scalarField_fourier , scalarField_real )
scalarField_real = scalarField_real * wgt ! normalize the result by number of elements
2015-06-03 23:00:31 +05:30
end subroutine utilities_FFTscalarBackward
2015-09-24 23:08:49 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed
2015-09-11 14:22:03 +05:30
!> @details Does an unweighted filtered FFT transform from real to complex.
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorForward ( )
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid3 , &
grid
2015-10-15 00:06:19 +05:30
2015-09-11 14:22:03 +05:30
implicit none
2015-09-11 18:35:46 +05:30
!--------------------------------------------------------------------------------------------------
2015-09-11 14:22:03 +05:30
! doing the vector FFT
call fftw_mpi_execute_dft_r2c ( planVectorForth , vectorField_real , vectorField_fourier )
2015-06-03 23:00:31 +05:30
end subroutine utilities_FFTvectorForward
2015-09-24 23:08:49 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief backward FFT of data in field_fourier to field_real
2015-09-11 14:22:03 +05:30
!> @details Does an weighted inverse FFT transform from complex to real
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorBackward ( )
2015-09-11 14:22:03 +05:30
implicit none
call fftw_mpi_execute_dft_c2r ( planVectorBack , vectorField_fourier , vectorField_real )
vectorField_real = vectorField_real * wgt ! normalize the result by number of elements
2015-06-03 23:00:31 +05:30
end subroutine utilities_FFTvectorBackward
2012-07-26 19:28:47 +05:30
2013-07-26 21:55:37 +05:30
2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
2012-10-24 17:01:40 +05:30
!> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim
2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
2015-06-03 23:00:31 +05:30
subroutine utilities_fourierGammaConvolution ( fieldAim )
2012-07-30 19:36:22 +05:30
use numerics , only : &
memory_efficient
2012-10-24 17:01:40 +05:30
use math , only : &
2015-12-14 23:42:09 +05:30
math_det33 , &
math_invert
2015-03-25 21:36:19 +05:30
use numerics , only : &
2015-10-15 00:06:19 +05:30
worldrank
2015-03-25 21:36:19 +05:30
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid3 , &
grid , &
grid3Offset
2012-07-26 19:28:47 +05:30
2015-10-15 00:06:19 +05:30
implicit none
2012-10-24 17:01:40 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: fieldAim !< desired average value of the field after convolution
2015-12-14 23:42:09 +05:30
complex ( pReal ) , dimension ( 3 , 3 ) :: temp33_complex , xiDyad_cmplx
real ( pReal ) :: matA ( 6 , 6 ) , matInvA ( 6 , 6 )
2015-10-15 00:06:19 +05:30
2012-10-24 17:01:40 +05:30
integer ( pInt ) :: &
i , j , k , &
l , m , n , o
2015-12-14 23:42:09 +05:30
logical :: ierr
2012-10-02 20:56:56 +05:30
2015-10-15 00:06:19 +05:30
if ( worldrank == 0_pInt ) then
2015-09-24 23:08:49 +05:30
write ( 6 , '(/,a)' ) ' ... doing gamma convolution ...............................................'
2015-03-25 21:36:19 +05:30
flush ( 6 )
endif
2015-10-15 00:06:19 +05:30
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
2013-04-18 22:10:49 +05:30
! do the actual spectral method calculation (mechanical equilibrium)
2015-09-11 14:22:03 +05:30
memoryEfficient : if ( memory_efficient ) then
2015-12-14 23:42:09 +05:30
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid1Red
if ( any ( [ i , j , k + grid3Offset ] / = 1_pInt ) ) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
2012-11-29 00:16:07 +05:30
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
2015-12-14 23:42:09 +05:30
xiDyad_cmplx ( l , m ) = conjg ( - xi1st ( l , i , j , k ) ) * xi1st ( m , i , j , k )
2012-11-29 00:16:07 +05:30
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
2015-12-14 23:42:09 +05:30
temp33_complex ( l , m ) = sum ( C_ref ( l , 1 : 3 , m , 1 : 3 ) * xiDyad_cmplx )
matA ( 1 : 3 , 1 : 3 ) = real ( temp33_complex ) ; matA ( 4 : 6 , 4 : 6 ) = real ( temp33_complex )
matA ( 1 : 3 , 4 : 6 ) = aimag ( temp33_complex ) ; matA ( 4 : 6 , 1 : 3 ) = - aimag ( temp33_complex )
if ( abs ( math_det33 ( matA ( 1 : 3 , 1 : 3 ) ) ) > 1e-16 ) then
call math_invert ( 6_pInt , matA , matInvA , ierr )
temp33_complex = cmplx ( matInvA ( 1 : 3 , 1 : 3 ) , matInvA ( 1 : 3 , 4 : 6 ) , pReal )
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt , n = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
gamma_hat ( l , m , n , o , 1 , 1 , 1 ) = temp33_complex ( l , n ) * conjg ( - xi1st ( o , i , j , k ) ) * xi1st ( m , i , j , k )
else
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt , n = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
gamma_hat ( l , m , n , o , 1 , 1 , 1 ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
endif
2012-11-29 00:16:07 +05:30
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
2015-12-14 23:42:09 +05:30
temp33_Complex ( l , m ) = sum ( gamma_hat ( l , m , 1 : 3 , 1 : 3 , 1 , 1 , 1 ) * tensorField_fourier ( 1 : 3 , 1 : 3 , i , j , k ) )
2015-10-15 00:06:19 +05:30
tensorField_fourier ( 1 : 3 , 1 : 3 , i , j , k ) = temp33_Complex
endif
2012-07-30 19:36:22 +05:30
enddo ; enddo ; enddo
2015-09-11 14:22:03 +05:30
else memoryEfficient
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid1Red
2013-03-01 15:06:45 +05:30
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
2015-06-03 23:00:31 +05:30
temp33_Complex ( l , m ) = sum ( gamma_hat ( l , m , 1 : 3 , 1 : 3 , i , j , k ) * &
2015-12-14 23:42:09 +05:30
tensorField_fourier ( 1 : 3 , 1 : 3 , i , j , k ) )
2015-09-11 14:22:03 +05:30
tensorField_fourier ( 1 : 3 , 1 : 3 , i , j , k ) = temp33_Complex
2012-07-30 19:36:22 +05:30
enddo ; enddo ; enddo
2015-09-11 14:22:03 +05:30
endif memoryEfficient
2015-10-15 00:06:19 +05:30
2015-09-11 14:22:03 +05:30
if ( grid3Offset == 0_pInt ) &
2015-12-14 23:42:09 +05:30
tensorField_fourier ( 1 : 3 , 1 : 3 , 1 , 1 , 1 ) = cmplx ( fieldAim / wgt , 0.0_pReal , pReal )
2012-07-30 19:36:22 +05:30
2015-06-03 23:00:31 +05:30
end subroutine utilities_fourierGammaConvolution
2015-10-15 00:06:19 +05:30
2015-09-24 23:08:49 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief doing convolution DamageGreenOp_hat * field_real
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierGreenConvolution ( D_ref , mobility_ref , deltaT )
use math , only : &
math_mul33x3 , &
PI
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid , &
grid3 , &
geomSize
2015-06-03 23:00:31 +05:30
2015-10-15 00:06:19 +05:30
implicit none
2015-06-03 23:00:31 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: D_ref !< desired average value of the field after convolution
real ( pReal ) , intent ( in ) :: mobility_ref , deltaT !< desired average value of the field after convolution
2015-12-14 23:42:09 +05:30
complex ( pReal ) :: GreenOp_hat
2015-06-03 23:00:31 +05:30
integer ( pInt ) :: i , j , k
2015-10-15 00:06:19 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation
2015-09-11 14:22:03 +05:30
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid1Red
2015-12-14 23:42:09 +05:30
GreenOp_hat = cmplx ( 1.0_pReal , 0.0_pReal , pReal ) / &
( cmplx ( mobility_ref , 0.0_pReal , pReal ) + &
deltaT * sum ( conjg ( - xi1st ( 1 : 3 , i , j , k ) ) * matmul ( D_ref , xi1st ( 1 : 3 , i , j , k ) ) ) )
2015-10-15 00:06:19 +05:30
scalarField_fourier ( i , j , k ) = scalarField_fourier ( i , j , k ) * GreenOp_hat
2015-06-03 23:00:31 +05:30
enddo ; enddo ; enddo
end subroutine utilities_fourierGreenConvolution
2012-07-30 19:36:22 +05:30
2015-09-24 23:08:49 +05:30
2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate root mean square of divergence of field_fourier
!--------------------------------------------------------------------------------------------------
2015-10-15 00:06:19 +05:30
real ( pReal ) function utilities_divergenceRMS ( )
2015-03-25 21:36:19 +05:30
use numerics , only : &
2015-10-15 00:06:19 +05:30
worldrank
2015-03-25 21:36:19 +05:30
use mesh , only : &
2015-12-14 23:42:09 +05:30
geomSize , &
2015-09-11 14:22:03 +05:30
grid , &
grid3
2012-10-24 17:01:40 +05:30
implicit none
2015-10-15 00:06:19 +05:30
integer ( pInt ) :: i , j , k
2015-03-25 21:36:19 +05:30
PetscErrorCode :: ierr
2012-07-30 19:36:22 +05:30
2015-10-15 00:06:19 +05:30
if ( worldrank == 0_pInt ) then
2015-03-25 21:36:19 +05:30
write ( 6 , '(/,a)' ) ' ... calculating divergence ................................................'
flush ( 6 )
endif
2013-04-10 15:49:16 +05:30
2012-07-19 22:54:56 +05:30
!--------------------------------------------------------------------------------------------------
2012-07-26 19:28:47 +05:30
! calculating RMS divergence criterion in Fourier space
2012-10-24 17:01:40 +05:30
utilities_divergenceRMS = 0.0_pReal
2015-09-11 14:22:03 +05:30
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 )
2013-05-08 21:22:29 +05:30
do i = 2_pInt , grid1Red - 1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
2012-10-24 17:01:40 +05:30
utilities_divergenceRMS = utilities_divergenceRMS &
2015-12-14 23:42:09 +05:30
+ 2.0_pReal * ( sum ( real ( matmul ( tensorField_fourier ( 1 : 3 , 1 : 3 , i , j , k ) , & ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
conjg ( - xi1st ( 1 : 3 , i , j , k ) ) * geomSize / scaledGeomSize ) ) ** 2.0_pReal ) & ! --> sum squared L_2 norm of vector
+ sum ( aimag ( matmul ( tensorField_fourier ( 1 : 3 , 1 : 3 , i , j , k ) , &
conjg ( - xi1st ( 1 : 3 , i , j , k ) ) * geomSize / scaledGeomSize ) ) ** 2.0_pReal ) )
2012-10-24 17:01:40 +05:30
enddo
2013-05-08 21:22:29 +05:30
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1)
2015-12-14 23:42:09 +05:30
+ sum ( real ( matmul ( tensorField_fourier ( 1 : 3 , 1 : 3 , 1 , j , k ) , &
conjg ( - xi1st ( 1 : 3 , 1 , j , k ) ) * geomSize / scaledGeomSize ) ) ** 2.0_pReal ) &
+ sum ( aimag ( matmul ( tensorField_fourier ( 1 : 3 , 1 : 3 , 1 , j , k ) , &
conjg ( - xi1st ( 1 : 3 , 1 , j , k ) ) * geomSize / scaledGeomSize ) ) ** 2.0_pReal ) &
+ sum ( real ( matmul ( tensorField_fourier ( 1 : 3 , 1 : 3 , grid1Red , j , k ) , &
conjg ( - xi1st ( 1 : 3 , grid1Red , j , k ) ) * geomSize / scaledGeomSize ) ) ** 2.0_pReal ) &
+ sum ( aimag ( matmul ( tensorField_fourier ( 1 : 3 , 1 : 3 , grid1Red , j , k ) , &
conjg ( - xi1st ( 1 : 3 , grid1Red , j , k ) ) * geomSize / scaledGeomSize ) ) ** 2.0_pReal )
2012-10-24 17:01:40 +05:30
enddo ; enddo
2015-09-11 14:22:03 +05:30
if ( grid ( 1 ) == 1_pInt ) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
2015-03-25 21:36:19 +05:30
call MPI_Allreduce ( MPI_IN_PLACE , utilities_divergenceRMS , 1 , MPI_DOUBLE , MPI_SUM , PETSC_COMM_WORLD , ierr )
2013-04-10 15:49:16 +05:30
utilities_divergenceRMS = sqrt ( utilities_divergenceRMS ) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
2012-10-24 17:01:40 +05:30
end function utilities_divergenceRMS
2015-10-15 00:06:19 +05:30
2013-03-06 20:01:13 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate max of curl of field_fourier
!--------------------------------------------------------------------------------------------------
real ( pReal ) function utilities_curlRMS ( )
2015-03-25 21:36:19 +05:30
use numerics , only : &
2015-10-15 00:06:19 +05:30
worldrank
2015-03-25 21:36:19 +05:30
use mesh , only : &
2015-12-14 23:42:09 +05:30
geomSize , &
2015-09-11 14:22:03 +05:30
grid , &
grid3
2013-03-06 20:01:13 +05:30
implicit none
2015-10-15 00:06:19 +05:30
integer ( pInt ) :: i , j , k , l
2013-03-06 20:01:13 +05:30
complex ( pReal ) , dimension ( 3 , 3 ) :: curl_fourier
2015-03-25 21:36:19 +05:30
PetscErrorCode :: ierr
2013-03-06 20:01:13 +05:30
2015-10-15 00:06:19 +05:30
if ( worldrank == 0_pInt ) then
2015-03-25 21:36:19 +05:30
write ( 6 , '(/,a)' ) ' ... calculating curl ......................................................'
flush ( 6 )
endif
!--------------------------------------------------------------------------------------------------
2013-03-06 20:01:13 +05:30
! calculating max curl criterion in Fourier space
utilities_curlRMS = 0.0_pReal
2015-10-15 00:06:19 +05:30
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ;
do i = 2_pInt , grid1Red - 1_pInt
2015-09-11 14:49:13 +05:30
do l = 1_pInt , 3_pInt
2015-12-14 23:42:09 +05:30
curl_fourier ( l , 1 ) = ( + tensorField_fourier ( l , 3 , i , j , k ) * xi1st ( 2 , i , j , k ) * geomSize ( 2 ) / scaledGeomSize ( 2 ) &
- tensorField_fourier ( l , 2 , i , j , k ) * xi1st ( 3 , i , j , k ) * geomSize ( 3 ) / scaledGeomSize ( 3 ) )
curl_fourier ( l , 2 ) = ( + tensorField_fourier ( l , 1 , i , j , k ) * xi1st ( 3 , i , j , k ) * geomSize ( 3 ) / scaledGeomSize ( 3 ) &
- tensorField_fourier ( l , 3 , i , j , k ) * xi1st ( 1 , i , j , k ) * geomSize ( 1 ) / scaledGeomSize ( 1 ) )
curl_fourier ( l , 3 ) = ( + tensorField_fourier ( l , 2 , i , j , k ) * xi1st ( 1 , i , j , k ) * geomSize ( 1 ) / scaledGeomSize ( 1 ) &
- tensorField_fourier ( l , 1 , i , j , k ) * xi1st ( 2 , i , j , k ) * geomSize ( 2 ) / scaledGeomSize ( 2 ) )
2015-09-11 14:49:13 +05:30
enddo
utilities_curlRMS = utilities_curlRMS + &
2.0_pReal * sum ( real ( curl_fourier ) ** 2.0_pReal + aimag ( curl_fourier ) ** 2.0_pReal ) ! Has somewhere a conj. complex counterpart. Therefore count it twice.
2015-10-15 00:06:19 +05:30
enddo
do l = 1_pInt , 3_pInt
2015-12-14 23:42:09 +05:30
curl_fourier = ( + tensorField_fourier ( l , 3 , 1 , j , k ) * xi1st ( 2 , 1 , j , k ) * geomSize ( 2 ) / scaledGeomSize ( 2 ) &
- tensorField_fourier ( l , 2 , 1 , j , k ) * xi1st ( 3 , 1 , j , k ) * geomSize ( 3 ) / scaledGeomSize ( 3 ) )
curl_fourier = ( + tensorField_fourier ( l , 1 , 1 , j , k ) * xi1st ( 3 , 1 , j , k ) * geomSize ( 3 ) / scaledGeomSize ( 3 ) &
- tensorField_fourier ( l , 3 , 1 , j , k ) * xi1st ( 1 , 1 , j , k ) * geomSize ( 1 ) / scaledGeomSize ( 1 ) )
curl_fourier = ( + tensorField_fourier ( l , 2 , 1 , j , k ) * xi1st ( 1 , 1 , j , k ) * geomSize ( 1 ) / scaledGeomSize ( 1 ) &
- tensorField_fourier ( l , 1 , 1 , j , k ) * xi1st ( 2 , 1 , j , k ) * geomSize ( 2 ) / scaledGeomSize ( 2 ) )
2013-03-06 20:01:13 +05:30
enddo
utilities_curlRMS = utilities_curlRMS + &
2015-09-11 14:49:13 +05:30
sum ( real ( curl_fourier ) ** 2.0_pReal + aimag ( curl_fourier ) ** 2.0_pReal ) ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1)
2015-10-15 00:06:19 +05:30
do l = 1_pInt , 3_pInt
2015-12-14 23:42:09 +05:30
curl_fourier = ( + tensorField_fourier ( l , 3 , grid1Red , j , k ) * xi1st ( 2 , grid1Red , j , k ) * geomSize ( 2 ) / scaledGeomSize ( 2 ) &
- tensorField_fourier ( l , 2 , grid1Red , j , k ) * xi1st ( 3 , grid1Red , j , k ) * geomSize ( 3 ) / scaledGeomSize ( 3 ) )
curl_fourier = ( + tensorField_fourier ( l , 1 , grid1Red , j , k ) * xi1st ( 3 , grid1Red , j , k ) * geomSize ( 3 ) / scaledGeomSize ( 3 ) &
- tensorField_fourier ( l , 3 , grid1Red , j , k ) * xi1st ( 1 , grid1Red , j , k ) * geomSize ( 1 ) / scaledGeomSize ( 1 ) )
curl_fourier = ( + tensorField_fourier ( l , 2 , grid1Red , j , k ) * xi1st ( 1 , grid1Red , j , k ) * geomSize ( 1 ) / scaledGeomSize ( 1 ) &
- tensorField_fourier ( l , 1 , grid1Red , j , k ) * xi1st ( 2 , grid1Red , j , k ) * geomSize ( 2 ) / scaledGeomSize ( 2 ) )
2015-09-11 14:49:13 +05:30
enddo
utilities_curlRMS = utilities_curlRMS + &
2015-10-15 00:06:19 +05:30
sum ( real ( curl_fourier ) ** 2.0_pReal + aimag ( curl_fourier ) ** 2.0_pReal ) ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1)
2013-03-06 20:01:13 +05:30
enddo ; enddo
2015-09-11 14:49:13 +05:30
2015-03-25 21:36:19 +05:30
call MPI_Allreduce ( MPI_IN_PLACE , utilities_curlRMS , 1 , MPI_DOUBLE , MPI_SUM , PETSC_COMM_WORLD , ierr )
2013-09-20 21:47:25 +05:30
utilities_curlRMS = sqrt ( utilities_curlRMS ) * wgt
2015-09-11 14:49:13 +05:30
if ( grid ( 1 ) == 1_pInt ) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
2013-03-06 20:01:13 +05:30
end function utilities_curlRMS
2012-07-24 22:37:10 +05:30
2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
2013-04-18 22:10:49 +05:30
!> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC
2012-07-30 19:36:22 +05:30
!--------------------------------------------------------------------------------------------------
2012-10-24 17:01:40 +05:30
function utilities_maskedCompliance ( rot_BC , mask_stress , C )
use IO , only : &
IO_error
2015-03-25 21:36:19 +05:30
use numerics , only : &
worldrank
2012-10-24 17:01:40 +05:30
use math , only : &
math_Plain3333to99 , &
math_plain99to3333 , &
math_rotate_forward3333 , &
math_rotate_forward33 , &
math_invert
2012-08-28 22:29:45 +05:30
implicit none
2015-03-25 21:36:19 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: utilities_maskedCompliance !< masked compliance
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 , 3 , 3 ) :: C !< current average stiffness
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: rot_BC !< rotation of load frame
logical , intent ( in ) , dimension ( 3 , 3 ) :: mask_stress !< mask of stress BC
2015-10-15 00:06:19 +05:30
integer ( pInt ) :: j , k , m , n
2012-10-19 14:14:21 +05:30
logical , dimension ( 9 ) :: mask_stressVector
2015-10-15 00:06:19 +05:30
real ( pReal ) , dimension ( 9 , 9 ) :: temp99_Real
integer ( pInt ) :: size_reduced = 0_pInt
2012-10-24 17:01:40 +05:30
real ( pReal ) , dimension ( : , : ) , allocatable :: &
s_reduced , & !< reduced compliance matrix (depending on number of stress BC)
2015-10-15 00:06:19 +05:30
c_reduced , & !< reduced stiffness (depending on number of stress BC)
2012-10-24 17:01:40 +05:30
sTimesC !< temp variable to check inversion
2012-07-30 19:36:22 +05:30
logical :: errmatinv
2012-08-28 22:29:45 +05:30
character ( len = 1024 ) :: formatString
2015-10-15 00:06:19 +05:30
2012-10-19 14:14:21 +05:30
mask_stressVector = reshape ( transpose ( mask_stress ) , [ 9 ] )
size_reduced = int ( count ( mask_stressVector ) , pInt )
2012-07-30 19:36:22 +05:30
if ( size_reduced > 0_pInt ) then
allocate ( c_reduced ( size_reduced , size_reduced ) , source = 0.0_pReal )
allocate ( s_reduced ( size_reduced , size_reduced ) , source = 0.0_pReal )
2012-08-28 22:29:45 +05:30
allocate ( sTimesC ( size_reduced , size_reduced ) , source = 0.0_pReal )
2012-07-30 19:36:22 +05:30
2012-10-24 17:01:40 +05:30
temp99_Real = math_Plain3333to99 ( math_rotate_forward3333 ( C , rot_BC ) )
2013-01-10 03:49:32 +05:30
2015-10-15 00:06:19 +05:30
if ( debugGeneral . and . worldrank == 0_pInt ) then
2013-01-10 03:49:32 +05:30
write ( 6 , '(/,a)' ) ' ... updating masked compliance ............................................'
2013-04-18 22:10:49 +05:30
write ( 6 , '(/,a,/,9(9(2x,f12.7,1x)/))' , advance = 'no' ) ' Stiffness C (load) / GPa =' , &
2012-10-24 17:01:40 +05:30
transpose ( temp99_Real ) / 1.e9_pReal
2013-01-10 03:49:32 +05:30
flush ( 6 )
endif
2012-10-24 17:01:40 +05:30
k = 0_pInt ! calculate reduced stiffness
2012-07-26 19:28:47 +05:30
do n = 1_pInt , 9_pInt
if ( mask_stressVector ( n ) ) then
k = k + 1_pInt
j = 0_pInt
do m = 1_pInt , 9_pInt
if ( mask_stressVector ( m ) ) then
j = j + 1_pInt
2012-07-30 19:36:22 +05:30
c_reduced ( k , j ) = temp99_Real ( n , m )
endif ; enddo ; endif ; enddo
2015-10-15 00:06:19 +05:30
2012-08-28 22:29:45 +05:30
call math_invert ( size_reduced , c_reduced , s_reduced , errmatinv ) ! invert reduced stiffness
2012-12-14 23:00:22 +05:30
if ( errmatinv ) call IO_error ( error_ID = 400_pInt , ext_msg = 'utilities_maskedCompliance' )
2012-10-24 17:01:40 +05:30
temp99_Real = 0.0_pReal ! fill up compliance with zeros
2012-07-30 19:36:22 +05:30
k = 0_pInt
do n = 1_pInt , 9_pInt
if ( mask_stressVector ( n ) ) then
k = k + 1_pInt
j = 0_pInt
do m = 1_pInt , 9_pInt
if ( mask_stressVector ( m ) ) then
j = j + 1_pInt
temp99_Real ( n , m ) = s_reduced ( k , j )
endif ; enddo ; endif ; enddo
2015-10-15 00:06:19 +05:30
2012-10-24 17:01:40 +05:30
!--------------------------------------------------------------------------------------------------
2013-03-07 01:04:30 +05:30
! check if inversion was successful
2012-08-28 22:29:45 +05:30
sTimesC = matmul ( c_reduced , s_reduced )
do m = 1_pInt , size_reduced
do n = 1_pInt , size_reduced
2012-10-24 17:01:40 +05:30
if ( m == n . and . abs ( sTimesC ( m , n ) ) > ( 1.0_pReal + 1 0.0e-12_pReal ) ) errmatinv = . true . ! diagonal elements of S*C should be 1
if ( m / = n . and . abs ( sTimesC ( m , n ) ) > ( 0.0_pReal + 1 0.0e-12_pReal ) ) errmatinv = . true . ! off diagonal elements of S*C should be 0
2012-08-28 22:29:45 +05:30
enddo
enddo
2015-03-25 21:36:19 +05:30
if ( ( debugGeneral . or . errmatinv ) . and . ( worldrank == 0_pInt ) ) then ! report
2012-08-28 22:29:45 +05:30
write ( formatString , '(I16.16)' ) size_reduced
2013-01-10 03:49:32 +05:30
formatString = '(/,a,/,' / / trim ( formatString ) / / '(' / / trim ( formatString ) / / '(2x,es9.2,1x)/))'
2013-04-18 22:10:49 +05:30
write ( 6 , trim ( formatString ) , advance = 'no' ) ' C * S (load) ' , &
transpose ( matmul ( c_reduced , s_reduced ) )
write ( 6 , trim ( formatString ) , advance = 'no' ) ' S (load) ' , transpose ( s_reduced )
2012-08-28 22:29:45 +05:30
endif
2012-12-14 23:00:22 +05:30
if ( errmatinv ) call IO_error ( error_ID = 400_pInt , ext_msg = 'utilities_maskedCompliance' )
2012-07-30 19:36:22 +05:30
deallocate ( c_reduced )
deallocate ( s_reduced )
2012-08-28 22:29:45 +05:30
deallocate ( sTimesC )
2012-07-30 19:36:22 +05:30
else
temp99_real = 0.0_pReal
endif
2015-03-25 21:36:19 +05:30
if ( debugGeneral . and . worldrank == 0_pInt ) & ! report
2013-04-18 22:10:49 +05:30
write ( 6 , '(/,a,/,9(9(2x,f12.7,1x)/),/)' , advance = 'no' ) ' Masked Compliance (load) * GPa =' , &
2013-01-10 03:49:32 +05:30
transpose ( temp99_Real * 1.e9_pReal )
flush ( 6 )
2012-10-24 17:01:40 +05:30
utilities_maskedCompliance = math_Plain99to3333 ( temp99_Real )
2015-10-15 00:06:19 +05:30
end function utilities_maskedCompliance
2012-10-24 17:01:40 +05:30
2015-06-03 23:00:31 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate scalar gradient in fourier field
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierScalarGradient ( )
use math , only : &
PI
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid3 , &
grid , &
geomSize
implicit none
2015-06-03 23:00:31 +05:30
integer ( pInt ) :: i , j , k
2015-09-11 14:22:03 +05:30
vectorField_fourier = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid1Red
2015-12-14 23:42:09 +05:30
vectorField_fourier ( 1 : 3 , i , j , k ) = scalarField_fourier ( i , j , k ) * xi1st ( 1 : 3 , i , j , k )
2015-06-03 23:00:31 +05:30
enddo ; enddo ; enddo
end subroutine utilities_fourierScalarGradient
!--------------------------------------------------------------------------------------------------
!> @brief calculate vector divergence in fourier field
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierVectorDivergence ( )
use math , only : &
PI
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid3 , &
grid , &
geomSize
implicit none
2015-06-03 23:00:31 +05:30
integer ( pInt ) :: i , j , k , m
2015-09-11 14:22:03 +05:30
scalarField_fourier = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid1Red
2015-06-03 23:00:31 +05:30
do m = 1_pInt , 3_pInt
2015-09-11 14:22:03 +05:30
scalarField_fourier ( i , j , k ) = &
scalarField_fourier ( i , j , k ) + &
2015-12-14 23:42:09 +05:30
vectorField_fourier ( m , i , j , k ) * conjg ( - xi1st ( m , i , j , k ) )
2015-06-03 23:00:31 +05:30
enddo
enddo ; enddo ; enddo
end subroutine utilities_fourierVectorDivergence
2015-09-11 18:35:46 +05:30
2015-12-14 23:42:09 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate vector gradient in fourier field
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierVectorGradient ( )
use mesh , only : &
grid3 , &
grid , &
geomSize
implicit none
integer ( pInt ) :: i , j , k , m , n
tensorField_fourier = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid1Red
do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt
tensorField_fourier ( m , n , i , j , k ) = vectorField_fourier ( m , i , j , k ) * xi1st ( n , i , j , k )
enddo ; enddo
enddo ; enddo ; enddo
end subroutine utilities_fourierVectorGradient
!--------------------------------------------------------------------------------------------------
!> @brief calculate tensor divergence in fourier field
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierTensorDivergence ( )
use mesh , only : &
grid3 , &
grid , &
geomSize
implicit none
integer ( pInt ) :: i , j , k , m , n
vectorField_fourier = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid1Red
do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt
vectorField_fourier ( m , i , j , k ) = &
vectorField_fourier ( m , i , j , k ) + &
tensorField_fourier ( m , n , i , j , k ) * conjg ( - xi1st ( n , i , j , k ) )
enddo ; enddo
enddo ; enddo ; enddo
end subroutine utilities_fourierTensorDivergence
2012-10-24 17:01:40 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates constitutive response
!--------------------------------------------------------------------------------------------------
2015-10-15 00:06:19 +05:30
subroutine utilities_constitutiveResponse ( F_lastInc , F , timeinc , &
2013-03-06 20:01:13 +05:30
P , C_volAvg , C_minmaxAvg , P_av , forwardData , rotation_BC )
2012-10-24 17:01:40 +05:30
use debug , only : &
debug_reset , &
debug_info
2015-03-25 21:36:19 +05:30
use numerics , only : &
worldrank
2012-10-24 17:01:40 +05:30
use math , only : &
math_transpose33 , &
2013-03-07 01:04:30 +05:30
math_rotate_forward33 , &
math_det33
2015-03-25 21:36:19 +05:30
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid , &
grid3
2012-10-24 17:01:40 +05:30
use FEsolving , only : &
restartWrite
use CPFEM , only : &
2013-03-01 17:18:29 +05:30
CPFEM_general , &
CPFEM_COLLECT , &
CPFEM_CALCRESULTS , &
2014-07-24 17:56:01 +05:30
CPFEM_AGERESULTS
2013-01-24 01:26:45 +05:30
use homogenization , only : &
materialpoint_F0 , &
materialpoint_F , &
materialpoint_P , &
materialpoint_dPdF
2015-10-15 00:06:19 +05:30
2012-10-24 17:01:40 +05:30
implicit none
2015-09-11 14:22:03 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ) :: &
2012-10-24 17:01:40 +05:30
F_lastInc , & !< target deformation gradient
F !< previous deformation gradient
real ( pReal ) , intent ( in ) :: timeinc !< loading time
logical , intent ( in ) :: forwardData !< age results
2013-05-08 21:22:29 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: rotation_BC !< rotation of load frame
2015-10-15 00:06:19 +05:30
2013-05-08 21:22:29 +05:30
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
2015-09-11 14:22:03 +05:30
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ) :: P !< PK stress
2015-10-15 00:06:19 +05:30
2012-10-24 17:01:40 +05:30
integer ( pInt ) :: &
calcMode , & !< CPFEM mode for calculation
2013-04-10 15:49:16 +05:30
j , k
2013-03-06 20:01:13 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: max_dPdF , min_dPdF
2013-04-10 15:49:16 +05:30
real ( pReal ) :: max_dPdF_norm , min_dPdF_norm , defgradDetMin , defgradDetMax , defgradDet
2015-03-25 21:36:19 +05:30
PetscErrorCode :: ierr
2015-10-15 00:06:19 +05:30
2015-09-11 14:22:03 +05:30
external :: &
MPI_Allreduce
2012-10-24 17:01:40 +05:30
2015-10-15 00:06:19 +05:30
if ( worldrank == 0_pInt ) then
2015-03-25 21:36:19 +05:30
write ( 6 , '(/,a)' ) ' ... evaluating constitutive response ......................................'
flush ( 6 )
endif
2013-03-01 17:18:29 +05:30
calcMode = CPFEM_CALCRESULTS
2014-08-04 23:20:01 +05:30
2012-10-24 17:01:40 +05:30
if ( forwardData ) then ! aging results
2015-10-15 00:06:19 +05:30
calcMode = ior ( calcMode , CPFEM_AGERESULTS )
2015-09-11 14:22:03 +05:30
materialpoint_F0 = reshape ( F_lastInc , [ 3 , 3 , 1 , product ( grid ( 1 : 2 ) ) * grid3 ] )
2012-10-24 17:01:40 +05:30
endif
if ( cutBack ) then ! restore saved variables
2015-10-15 00:06:19 +05:30
calcMode = iand ( calcMode , not ( CPFEM_AGERESULTS ) )
2012-10-24 17:01:40 +05:30
endif
2013-03-07 01:04:30 +05:30
2014-08-04 23:20:01 +05:30
call CPFEM_general ( CPFEM_COLLECT , F_lastInc ( 1 : 3 , 1 : 3 , 1 , 1 , 1 ) , F ( 1 : 3 , 1 : 3 , 1 , 1 , 1 ) , &
2015-07-24 20:27:29 +05:30
timeinc , 1_pInt , 1_pInt )
2013-04-09 15:38:00 +05:30
2015-09-11 14:22:03 +05:30
materialpoint_F = reshape ( F , [ 3 , 3 , 1 , product ( grid ( 1 : 2 ) ) * grid3 ] )
2013-04-09 15:38:00 +05:30
call debug_reset ( )
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! calculate bounds of det(F) and report
2013-03-07 01:04:30 +05:30
if ( debugGeneral ) then
defgradDetMax = - huge ( 1.0_pReal )
defgradDetMin = + huge ( 1.0_pReal )
2015-09-11 14:22:03 +05:30
do j = 1_pInt , product ( grid ( 1 : 2 ) ) * grid3
2013-04-09 15:38:00 +05:30
defgradDet = math_det33 ( materialpoint_F ( 1 : 3 , 1 : 3 , 1 , j ) )
2013-03-07 01:04:30 +05:30
defgradDetMax = max ( defgradDetMax , defgradDet )
2015-10-15 00:06:19 +05:30
defgradDetMin = min ( defgradDetMin , defgradDet )
2013-04-09 15:38:00 +05:30
end do
2015-03-25 21:36:19 +05:30
call MPI_reduce ( MPI_IN_PLACE , defgradDetMax , 1 , MPI_DOUBLE , MPI_MAX , 0 , PETSC_COMM_WORLD , ierr )
call MPI_reduce ( MPI_IN_PLACE , defgradDetMin , 1 , MPI_DOUBLE , MPI_MIN , 0 , PETSC_COMM_WORLD , ierr )
2015-10-15 00:06:19 +05:30
if ( worldrank == 0_pInt ) then
2015-03-25 21:36:19 +05:30
write ( 6 , '(a,1x,es11.4)' ) ' max determinant of deformation =' , defgradDetMax
write ( 6 , '(a,1x,es11.4)' ) ' min determinant of deformation =' , defgradDetMin
flush ( 6 )
endif
2013-03-07 01:04:30 +05:30
endif
2015-06-03 23:00:31 +05:30
2015-03-25 21:36:19 +05:30
call CPFEM_general ( calcMode , F_lastInc ( 1 : 3 , 1 : 3 , 1 , 1 , 1 ) , F ( 1 : 3 , 1 : 3 , 1 , 1 , 1 ) , & ! first call calculates everything
2015-07-24 20:27:29 +05:30
timeinc , 1_pInt , 1_pInt )
2014-06-11 13:49:07 +05:30
2013-03-06 20:01:13 +05:30
max_dPdF = 0.0_pReal
max_dPdF_norm = 0.0_pReal
min_dPdF = huge ( 1.0_pReal )
min_dPdF_norm = huge ( 1.0_pReal )
2015-09-11 14:22:03 +05:30
do k = 1_pInt , product ( grid ( 1 : 2 ) ) * grid3
2013-03-06 20:01:13 +05:30
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 )
2015-10-15 00:06:19 +05:30
endif
2013-03-06 20:01:13 +05:30
if ( min_dPdF_norm > sum ( materialpoint_dPdF ( 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 , 1 , k ) ** 2.0_pReal ) ) then
min_dPdF = materialpoint_dPdF ( 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 , 1 , k )
min_dPdF_norm = sum ( materialpoint_dPdF ( 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 , 1 , k ) ** 2.0_pReal )
2015-10-15 00:06:19 +05:30
endif
2013-04-09 15:38:00 +05:30
end do
2015-10-15 00:06:19 +05:30
2015-03-25 21:36:19 +05:30
call MPI_Allreduce ( MPI_IN_PLACE , max_dPdF , 81 , MPI_DOUBLE , MPI_MAX , PETSC_COMM_WORLD , ierr )
call MPI_Allreduce ( MPI_IN_PLACE , min_dPdF , 81 , MPI_DOUBLE , MPI_MIN , PETSC_COMM_WORLD , ierr )
2013-01-24 01:26:45 +05:30
2015-10-15 00:06:19 +05:30
C_minmaxAvg = 0.5_pReal * ( max_dPdF + min_dPdF )
2013-03-06 20:01:13 +05:30
C_volAvg = sum ( sum ( materialpoint_dPdF , dim = 6 ) , dim = 5 ) * wgt
2015-10-15 00:06:19 +05:30
2015-03-25 21:36:19 +05:30
call MPI_Allreduce ( MPI_IN_PLACE , C_volAvg , 81 , MPI_DOUBLE , MPI_SUM , PETSC_COMM_WORLD , ierr )
2015-10-15 00:06:19 +05:30
2012-10-24 17:01:40 +05:30
call debug_info ( )
2015-10-15 00:06:19 +05:30
2012-10-24 17:01:40 +05:30
restartWrite = . false . ! reset restartWrite status
cutBack = . false . ! reset cutBack status
2015-10-15 00:06:19 +05:30
2015-09-11 14:22:03 +05:30
P = reshape ( materialpoint_P , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ] )
2015-10-15 00:06:19 +05:30
P_av = sum ( sum ( sum ( P , dim = 5 ) , dim = 4 ) , dim = 3 ) * wgt ! average of P
2015-03-25 21:36:19 +05:30
call MPI_Allreduce ( MPI_IN_PLACE , P_av , 9 , MPI_DOUBLE , MPI_SUM , PETSC_COMM_WORLD , ierr )
if ( debugRotation . and . worldrank == 0_pInt ) &
2013-07-08 21:18:13 +05:30
write ( 6 , '(/,a,/,3(3(2x,f12.4,1x)/))' , advance = 'no' ) ' Piola--Kirchhoff stress (lab) / MPa =' , &
2013-08-07 22:50:05 +05:30
math_transpose33 ( P_av ) * 1.e-6_pReal
2012-11-28 20:34:05 +05:30
P_av = math_rotate_forward33 ( P_av , rotation_BC )
2015-03-25 21:36:19 +05:30
if ( worldrank == 0_pInt ) then
write ( 6 , '(/,a,/,3(3(2x,f12.4,1x)/))' , advance = 'no' ) ' Piola--Kirchhoff stress / MPa =' , &
2013-08-07 22:50:05 +05:30
math_transpose33 ( P_av ) * 1.e-6_pReal
2015-03-25 21:36:19 +05:30
flush ( 6 )
endif
2013-01-10 03:49:32 +05:30
2012-10-24 17:01:40 +05:30
end subroutine utilities_constitutiveResponse
2012-07-30 19:36:22 +05:30
2012-10-24 17:01:40 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates forward rate, either guessing or just add delta/timeinc
!--------------------------------------------------------------------------------------------------
2012-11-28 20:34:05 +05:30
pure function utilities_calculateRate ( avRate , timeinc_old , guess , field_lastInc , field )
2015-03-25 21:36:19 +05:30
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid3 , &
grid
2015-10-15 00:06:19 +05:30
2012-08-28 22:29:45 +05:30
implicit none
2012-11-28 20:34:05 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: avRate !< homogeneous addon
2012-10-24 17:01:40 +05:30
real ( pReal ) , intent ( in ) :: &
2012-11-09 01:02:00 +05:30
timeinc_old !< timeinc of last step
logical , intent ( in ) :: &
guess !< guess along former trajectory
2015-09-11 14:22:03 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ) :: &
2012-10-24 17:01:40 +05:30
field_lastInc , & !< data of previous step
field !< data of current step
2015-09-11 14:22:03 +05:30
real ( pReal ) , dimension ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ) :: &
2015-03-25 21:36:19 +05:30
utilities_calculateRate
2015-10-15 00:06:19 +05:30
2015-03-25 21:36:19 +05:30
if ( guess ) then
2012-10-24 17:01:40 +05:30
utilities_calculateRate = ( field - field_lastInc ) / timeinc_old
2012-08-09 18:34:56 +05:30
else
2015-09-11 14:22:03 +05:30
utilities_calculateRate = spread ( spread ( spread ( avRate , 3 , grid ( 1 ) ) , 4 , grid ( 2 ) ) , 5 , grid3 )
2012-08-09 18:34:56 +05:30
endif
2012-10-24 17:01:40 +05:30
end function utilities_calculateRate
2012-10-02 20:56:56 +05:30
2012-10-24 17:01:40 +05:30
!--------------------------------------------------------------------------------------------------
2015-10-15 00:06:19 +05:30
!> @brief forwards a field with a pointwise given rate, if aim is given,
2013-01-08 15:42:03 +05:30
!> ensures that the average matches the aim
2012-10-24 17:01:40 +05:30
!--------------------------------------------------------------------------------------------------
2015-03-25 21:36:19 +05:30
function utilities_forwardField ( timeinc , field_lastInc , rate , aim )
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid3 , &
grid
2015-10-15 00:06:19 +05:30
2012-10-02 20:56:56 +05:30
implicit none
2015-10-15 00:06:19 +05:30
real ( pReal ) , intent ( in ) :: &
2013-01-08 15:42:03 +05:30
timeinc !< timeinc of current step
2015-09-11 14:22:03 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ) :: &
2013-01-08 15:42:03 +05:30
field_lastInc , & !< initial field
2012-10-24 17:01:40 +05:30
rate !< rate by which to forward
2013-01-08 15:42:03 +05:30
real ( pReal ) , intent ( in ) , optional , dimension ( 3 , 3 ) :: &
aim !< average field value aim
2015-09-11 14:22:03 +05:30
real ( pReal ) , dimension ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ) :: &
2015-03-25 21:36:19 +05:30
utilities_forwardField
2015-09-11 18:35:46 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: fieldDiff !< <a + adot*t> - aim
2015-03-25 21:36:19 +05:30
PetscErrorCode :: ierr
2015-10-15 00:06:19 +05:30
2015-09-11 14:22:03 +05:30
external :: &
MPI_Allreduce
2015-10-15 00:06:19 +05:30
2012-10-24 17:01:40 +05:30
utilities_forwardField = field_lastInc + rate * timeinc
2013-01-08 15:42:03 +05:30
if ( present ( aim ) ) then !< correct to match average
2015-03-25 21:36:19 +05:30
fieldDiff = sum ( sum ( sum ( utilities_forwardField , dim = 5 ) , dim = 4 ) , dim = 3 ) * wgt
call MPI_Allreduce ( MPI_IN_PLACE , fieldDiff , 9 , MPI_DOUBLE , MPI_SUM , PETSC_COMM_WORLD , ierr )
fieldDiff = fieldDiff - aim
2013-01-08 15:42:03 +05:30
utilities_forwardField = utilities_forwardField - &
2015-09-11 14:22:03 +05:30
spread ( spread ( spread ( fieldDiff , 3 , grid ( 1 ) ) , 4 , grid ( 2 ) ) , 5 , grid3 )
2013-01-08 15:42:03 +05:30
endif
2012-10-02 20:56:56 +05:30
2012-10-24 17:01:40 +05:30
end function utilities_forwardField
2012-10-02 20:56:56 +05:30
2012-08-07 22:53:13 +05:30
2012-10-24 17:01:40 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates filter for fourier convolution depending on type given in numerics.config
!--------------------------------------------------------------------------------------------------
2015-12-14 23:42:09 +05:30
pure function utilities_getFreqDerivative ( k_s )
2012-10-24 17:01:40 +05:30
use math , only : &
PI
2015-03-25 21:36:19 +05:30
use mesh , only : &
2015-12-14 23:42:09 +05:30
geomSize , &
2015-09-11 14:22:03 +05:30
grid
2012-10-24 17:01:40 +05:30
implicit none
2015-12-14 23:42:09 +05:30
integer ( pInt ) , intent ( in ) , dimension ( 3 ) :: k_s !< indices of frequency
complex ( pReal ) , dimension ( 3 ) :: utilities_getFreqDerivative
select case ( spectral_derivative_ID )
case ( DERIVATIVE_CONTINUOUS_ID )
utilities_getFreqDerivative = cmplx ( 0.0_pReal , 2.0_pReal * PI * real ( k_s , pReal ) / geomSize , pReal )
case ( DERIVATIVE_CENTRAL_DIFF_ID )
utilities_getFreqDerivative = cmplx ( 0.0_pReal , sin ( 2.0_pReal * PI * real ( k_s , pReal ) / real ( grid , pReal ) ) , pReal ) / &
cmplx ( 2.0_pReal * geomSize / real ( grid , pReal ) , 0.0_pReal , pReal )
case ( DERIVATIVE_FWBW_DIFF_ID )
utilities_getFreqDerivative ( 1 ) = &
cmplx ( cos ( 2.0_pReal * PI * real ( k_s ( 1 ) , pReal ) / real ( grid ( 1 ) , pReal ) ) - 1.0_pReal , &
sin ( 2.0_pReal * PI * real ( k_s ( 1 ) , pReal ) / real ( grid ( 1 ) , pReal ) ) , pReal ) * &
cmplx ( cos ( 2.0_pReal * PI * real ( k_s ( 2 ) , pReal ) / real ( grid ( 2 ) , pReal ) ) + 1.0_pReal , &
sin ( 2.0_pReal * PI * real ( k_s ( 2 ) , pReal ) / real ( grid ( 2 ) , pReal ) ) , pReal ) * &
cmplx ( cos ( 2.0_pReal * PI * real ( k_s ( 3 ) , pReal ) / real ( grid ( 3 ) , pReal ) ) + 1.0_pReal , &
sin ( 2.0_pReal * PI * real ( k_s ( 3 ) , pReal ) / real ( grid ( 3 ) , pReal ) ) , pReal ) / &
cmplx ( 4.0_pReal * geomSize ( 1 ) / real ( grid ( 1 ) , pReal ) , 0.0_pReal , pReal )
utilities_getFreqDerivative ( 2 ) = &
cmplx ( cos ( 2.0_pReal * PI * real ( k_s ( 1 ) , pReal ) / real ( grid ( 1 ) , pReal ) ) + 1.0_pReal , &
sin ( 2.0_pReal * PI * real ( k_s ( 1 ) , pReal ) / real ( grid ( 1 ) , pReal ) ) , pReal ) * &
cmplx ( cos ( 2.0_pReal * PI * real ( k_s ( 2 ) , pReal ) / real ( grid ( 2 ) , pReal ) ) - 1.0_pReal , &
sin ( 2.0_pReal * PI * real ( k_s ( 2 ) , pReal ) / real ( grid ( 2 ) , pReal ) ) , pReal ) * &
cmplx ( cos ( 2.0_pReal * PI * real ( k_s ( 3 ) , pReal ) / real ( grid ( 3 ) , pReal ) ) + 1.0_pReal , &
sin ( 2.0_pReal * PI * real ( k_s ( 3 ) , pReal ) / real ( grid ( 3 ) , pReal ) ) , pReal ) / &
cmplx ( 4.0_pReal * geomSize ( 2 ) / real ( grid ( 2 ) , pReal ) , 0.0_pReal , pReal )
utilities_getFreqDerivative ( 3 ) = &
cmplx ( cos ( 2.0_pReal * PI * real ( k_s ( 1 ) , pReal ) / real ( grid ( 1 ) , pReal ) ) + 1.0_pReal , &
sin ( 2.0_pReal * PI * real ( k_s ( 1 ) , pReal ) / real ( grid ( 1 ) , pReal ) ) , pReal ) * &
cmplx ( cos ( 2.0_pReal * PI * real ( k_s ( 2 ) , pReal ) / real ( grid ( 2 ) , pReal ) ) + 1.0_pReal , &
sin ( 2.0_pReal * PI * real ( k_s ( 2 ) , pReal ) / real ( grid ( 2 ) , pReal ) ) , pReal ) * &
cmplx ( cos ( 2.0_pReal * PI * real ( k_s ( 3 ) , pReal ) / real ( grid ( 3 ) , pReal ) ) - 1.0_pReal , &
sin ( 2.0_pReal * PI * real ( k_s ( 3 ) , pReal ) / real ( grid ( 3 ) , pReal ) ) , pReal ) / &
cmplx ( 4.0_pReal * geomSize ( 3 ) / real ( grid ( 3 ) , pReal ) , 0.0_pReal , pReal )
2015-09-11 18:35:46 +05:30
end select
2012-07-25 19:31:39 +05:30
2015-12-14 23:42:09 +05:30
end function utilities_getFreqDerivative
2012-10-02 20:56:56 +05:30
2015-03-13 03:58:33 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate coordinates in current configuration for given defgrad field
! using integration in Fourier space. Similar as in mesh.f90, but using data already defined for
! convolution
!--------------------------------------------------------------------------------------------------
subroutine utilities_updateIPcoords ( F )
2015-09-11 14:22:03 +05:30
use math , only : &
math_mul33x3
2015-03-13 03:58:33 +05:30
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid , &
grid3 , &
grid3Offset , &
geomSize , &
2015-10-15 00:06:19 +05:30
mesh_ipCoordinates
2015-03-13 03:58:33 +05:30
implicit none
2015-09-11 14:22:03 +05:30
real ( pReal ) , dimension ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ) , intent ( in ) :: F
2015-03-13 03:58:33 +05:30
integer ( pInt ) :: i , j , k , m
2015-12-14 23:42:09 +05:30
real ( pReal ) , dimension ( 3 ) :: step , offset_coords
2015-03-13 03:58:33 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: Favg
2015-03-25 21:36:19 +05:30
PetscErrorCode :: ierr
2015-09-11 14:22:03 +05:30
external &
MPI_Bcast
2015-03-13 03:58:33 +05:30
2015-12-14 23:42:09 +05:30
!--------------------------------------------------------------------------------------------------
! integration in Fourier space
2015-09-11 14:22:03 +05:30
tensorField_real = 0.0_pReal
tensorField_real ( 1 : 3 , 1 : 3 , 1 : grid ( 1 ) , 1 : grid ( 2 ) , 1 : grid3 ) = F
2015-06-03 23:00:31 +05:30
call utilities_FFTtensorForward ( )
2015-12-14 23:42:09 +05:30
call utilities_fourierTensorDivergence ( )
2015-03-13 03:58:33 +05:30
2015-12-14 23:42:09 +05:30
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid1Red
if ( any ( abs ( xi1st ( 1 : 3 , i , j , k ) ) > tiny ( 0.0_pReal ) ) ) &
vectorField_fourier ( 1 : 3 , i , j , k ) = vectorField_fourier ( 1 : 3 , i , j , k ) / &
sum ( conjg ( - xi1st ( 1 : 3 , i , j , k ) ) * xi1st ( 1 : 3 , i , j , k ) )
enddo ; enddo ; enddo
call fftw_mpi_execute_dft_c2r ( planVectorBack , vectorField_fourier , vectorField_real )
2015-10-15 00:06:19 +05:30
2015-03-13 03:58:33 +05:30
!--------------------------------------------------------------------------------------------------
! average F
2015-09-11 14:22:03 +05:30
if ( grid3Offset == 0_pInt ) Favg = real ( tensorField_fourier ( 1 : 3 , 1 : 3 , 1 , 1 , 1 ) , pReal ) * wgt
2015-03-25 21:36:19 +05:30
call MPI_Bcast ( Favg , 9 , MPI_DOUBLE , 0 , PETSC_COMM_WORLD , ierr )
2015-03-13 03:58:33 +05:30
!--------------------------------------------------------------------------------------------------
! add average to fluctuation and put (0,0,0) on (0,0,0)
2015-12-14 23:42:09 +05:30
step = geomSize / real ( grid , pReal )
2015-09-11 14:22:03 +05:30
if ( grid3Offset == 0_pInt ) offset_coords = vectorField_real ( 1 : 3 , 1 , 1 , 1 )
2015-03-25 21:36:19 +05:30
call MPI_Bcast ( offset_coords , 3 , MPI_DOUBLE , 0 , PETSC_COMM_WORLD , ierr )
offset_coords = math_mul33x3 ( Favg , step / 2.0_pReal ) - offset_coords
2015-03-13 03:58:33 +05:30
m = 1_pInt
2015-09-11 14:22:03 +05:30
do k = 1_pInt , grid3 ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid ( 1 )
mesh_ipCoordinates ( 1 : 3 , 1 , m ) = vectorField_real ( 1 : 3 , i , j , k ) &
2015-03-13 03:58:33 +05:30
+ offset_coords &
2015-09-11 14:22:03 +05:30
+ math_mul33x3 ( Favg , step * real ( [ i , j , k + grid3Offset ] - 1_pInt , pReal ) )
2015-03-13 03:58:33 +05:30
m = m + 1_pInt
enddo ; enddo ; enddo
end subroutine utilities_updateIPcoords
2015-09-11 18:35:46 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief cleans up
!--------------------------------------------------------------------------------------------------
subroutine utilities_destroy ( )
implicit none
call fftw_destroy_plan ( planTensorForth )
call fftw_destroy_plan ( planTensorBack )
call fftw_destroy_plan ( planVectorForth )
call fftw_destroy_plan ( planVectorBack )
call fftw_destroy_plan ( planScalarForth )
call fftw_destroy_plan ( planScalarBack )
end subroutine utilities_destroy
2012-10-24 17:01:40 +05:30
end module DAMASK_spectral_utilities