2013-03-22 23:05:05 +05:30
! Copyright 2011-13 Max-Planck-Institut für Eisenforschung GmbH
!
! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit.
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
2013-07-24 18:36:16 +05:30
! $Id: DAMASK_spectral_solverAL.f90 2487 2013-06-14 09:49:33Z MPIE\f.roters $
2012-08-06 14:23:12 +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 AL scheme solver
!--------------------------------------------------------------------------------------------------
2012-11-12 19:44:39 +05:30
module DAMASK_spectral_solverAL
2012-08-06 14:23:12 +05:30
use prec , only : &
pInt , &
pReal
use math , only : &
math_I3
2012-11-12 19:44:39 +05:30
use DAMASK_spectral_utilities , only : &
2012-10-02 20:56:56 +05:30
tSolutionState
2012-07-25 19:31:39 +05:30
2012-08-06 14:23:12 +05:30
implicit none
2013-03-28 16:07:00 +05:30
private
2012-08-06 18:13:05 +05:30
#include <finclude/petscsys.h>
#include <finclude/petscdmda.h>
#include <finclude/petscsnes.h>
2012-08-06 14:23:12 +05:30
character ( len = * ) , parameter , public :: &
2013-07-24 18:36:16 +05:30
DAMASK_spectral_solverAL_label = 'al'
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
2013-05-13 15:14:23 +05:30
! derived types
2013-06-11 22:05:04 +05:30
type tSolutionParams !< @todo use here the type definition for a full loadcase including mask
2012-08-06 14:23:12 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: P_BC , rotation_BC
real ( pReal ) :: timeinc
2013-07-26 21:55:37 +05:30
real ( pReal ) :: timeincOld
2013-01-09 23:38:08 +05:30
real ( pReal ) :: temperature
2013-07-26 21:55:37 +05:30
real ( pReal ) :: density
2012-10-02 20:56:56 +05:30
end type tSolutionParams
2012-07-26 19:28:47 +05:30
2012-10-02 20:56:56 +05:30
type ( tSolutionParams ) , private :: params
2012-11-12 19:44:39 +05:30
real ( pReal ) , private , dimension ( 3 , 3 ) :: mask_stress = 0.0_pReal
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
! PETSc data
2012-11-06 21:30:51 +05:30
DM , private :: da
2012-10-02 20:56:56 +05:30
SNES , private :: snes
2012-11-06 21:30:51 +05:30
Vec , private :: solution_vec
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
! common pointwise data
2012-11-12 19:44:39 +05:30
real ( pReal ) , private , dimension ( : , : , : , : , : ) , allocatable :: &
F_lastInc , & !< field of previous compatible deformation gradients
2013-07-26 21:55:37 +05:30
F_lastInc2 , & !< field of 2nd previous compatible deformation gradients
F_lambda_lastInc , & !< field of previous incompatible deformation gradient
2012-11-12 19:44:39 +05:30
Fdot , & !< field of assumed rate of compatible deformation gradient
2013-07-26 21:55:37 +05:30
F_lambdaDot !< field of assumed rate of incopatible deformation gradient
complex ( pReal ) , private , dimension ( : , : , : , : , : ) , allocatable :: inertiaField_fourier
2013-01-09 23:38:08 +05:30
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
2012-10-02 20:56:56 +05:30
real ( pReal ) , private , dimension ( 3 , 3 ) :: &
2012-11-12 19:44:39 +05:30
F_aimDot , & !< assumed rate of average deformation gradient
F_aim = math_I3 , & !< current prescribed deformation gradient
F_aim_lastInc = math_I3 , & !< previous average deformation gradient
2013-07-30 21:02:55 +05:30
P_av = 0.0_pReal , & !< average 1st Piola--Kirchhoff stress
P_avLastEval = 0.0_pReal !< average 1st Piola--Kirchhoff stress last call of CPFEM_general
2012-11-12 19:44:39 +05:30
character ( len = 1024 ) , private :: incInfo !< time and increment information
2012-10-02 20:56:56 +05:30
real ( pReal ) , private , dimension ( 3 , 3 , 3 , 3 ) :: &
2013-07-23 23:12:15 +05:30
C_volAvg = 0.0_pReal , & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal , & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal , & !< current (min+max)/2 stiffness
2012-11-12 19:44:39 +05:30
S = 0.0_pReal , & !< current compliance (filled up with zeros)
C_scale = 0.0_pReal , &
2012-10-02 20:56:56 +05:30
S_scale = 0.0_pReal
2012-07-25 19:31:39 +05:30
2012-11-12 19:44:39 +05:30
real ( pReal ) , private :: &
err_stress , & !< deviation from stress BC
err_f , & !< difference between compatible and incompatible deformation gradient
err_p !< difference of stress resulting from compatible and incompatible F
logical , private :: ForwardData
2013-07-23 23:12:15 +05:30
integer ( pInt ) , private :: &
totalIter = 0_pInt !< total iteration in current increment
2013-03-28 13:10:30 +05:30
public :: &
AL_init , &
AL_solution , &
AL_destroy
2013-03-27 17:58:55 +05:30
external :: &
2013-05-08 21:22:29 +05:30
VecDestroy , &
DMDestroy , &
DMDACreate3D , &
DMCreateGlobalVector , &
DMDASetLocalFunction , &
PETScFinalize , &
SNESDestroy , &
SNESGetNumberFunctionEvals , &
SNESGetIterationNumber , &
SNESSolve , &
SNESSetDM , &
SNESGetConvergedReason , &
SNESSetConvergenceTest , &
SNESSetFromOptions , &
SNESCreate , &
MPI_Abort
2013-03-27 17:58:55 +05:30
2012-10-02 20:56:56 +05:30
contains
2013-03-27 17:58:55 +05:30
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
2013-06-11 22:05:04 +05:30
!> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc)
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-09 23:38:08 +05:30
subroutine AL_init ( temperature )
2012-11-12 19:44:39 +05:30
use , intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO , only : &
2013-05-08 21:22:29 +05:30
IO_intOut , &
2012-11-12 19:44:39 +05:30
IO_read_JobBinaryFile , &
2013-02-25 22:04:59 +05:30
IO_write_JobBinaryFile , &
IO_timeStamp
2013-05-08 21:22:29 +05:30
use debug , only : &
debug_level , &
debug_spectral , &
debug_spectralRestart
2012-11-12 19:44:39 +05:30
use FEsolving , only : &
restartInc
use DAMASK_interface , only : &
getSolverJobName
use DAMASK_spectral_Utilities , only : &
Utilities_init , &
Utilities_constitutiveResponse , &
Utilities_updateGamma , &
2013-05-08 21:22:29 +05:30
grid , &
2013-07-26 21:55:37 +05:30
grid1Red , &
2013-05-08 21:22:29 +05:30
geomSize , &
wgt
2012-11-12 19:44:39 +05:30
use mesh , only : &
2013-02-05 20:33:36 +05:30
mesh_ipCoordinates , &
mesh_deformedCoordsFFT
2012-11-12 19:44:39 +05:30
use math , only : &
math_invSym3333
2012-07-25 19:31:39 +05:30
2012-11-12 19:44:39 +05:30
implicit none
2013-01-09 23:38:08 +05:30
real ( pReal ) , intent ( inout ) :: &
temperature
2012-11-07 18:41:41 +05:30
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90>
2013-05-08 21:22:29 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: P
2012-12-14 23:00:22 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: &
temp33_Real = 0.0_pReal
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: &
temp3333_Real = 0.0_pReal , &
temp3333_Real2 = 0.0_pReal
2012-11-12 19:44:39 +05:30
PetscErrorCode :: ierr
PetscObject :: dummy
2013-07-24 18:36:16 +05:30
PetscScalar , pointer , dimension ( : , : , : , : ) :: xx_psc , F , F_lambda
2012-11-12 19:44:39 +05:30
call Utilities_init ( )
write ( 6 , '(/,a)' ) ' <<<+- DAMASK_spectral_solverAL init -+>>>'
2013-07-24 18:36:16 +05:30
write ( 6 , '(a)' ) ' $Id: DAMASK_spectral_solverAL.f90 2487 2013-06-14 09:49:33Z MPIE\f.roters $'
2013-02-25 22:04:59 +05:30
write ( 6 , '(a16,a)' ) ' Current time : ' , IO_timeStamp ( )
2012-08-06 18:13:05 +05:30
#include "compilation_info.f90"
2013-05-08 21:22:29 +05:30
allocate ( P ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
!--------------------------------------------------------------------------------------------------
! allocate global fields
2013-06-11 22:05:04 +05:30
allocate ( F_lastInc ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
2013-07-26 21:55:37 +05:30
allocate ( F_lastInc2 ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
2013-05-08 21:22:29 +05:30
allocate ( Fdot ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
2013-07-24 18:36:16 +05:30
allocate ( F_lambda_lastInc ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
allocate ( F_lambdaDot ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
2013-07-26 21:55:37 +05:30
allocate ( inertiaField_fourier ( grid1Red , grid ( 2 ) , grid ( 3 ) , 3 , 3 ) , source = cmplx ( 0.0_pReal , 0.0_pReal , pReal ) )
2012-08-09 18:34:56 +05:30
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! PETSc Init
2012-12-17 15:48:39 +05:30
call SNESCreate ( PETSC_COMM_WORLD , snes , ierr ) ; CHKERRQ ( ierr )
2012-11-12 19:44:39 +05:30
call DMDACreate3d ( PETSC_COMM_WORLD , &
DMDA_BOUNDARY_NONE , DMDA_BOUNDARY_NONE , DMDA_BOUNDARY_NONE , &
2013-05-08 21:22:29 +05:30
DMDA_STENCIL_BOX , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) , PETSC_DECIDE , PETSC_DECIDE , PETSC_DECIDE , &
2012-11-12 19:44:39 +05:30
18 , 1 , PETSC_NULL_INTEGER , PETSC_NULL_INTEGER , PETSC_NULL_INTEGER , da , ierr )
CHKERRQ ( ierr )
2012-12-17 15:48:39 +05:30
call DMCreateGlobalVector ( da , solution_vec , ierr ) ; CHKERRQ ( ierr )
call DMDASetLocalFunction ( da , AL_formResidual , ierr ) ; CHKERRQ ( ierr )
call SNESSetDM ( snes , da , ierr ) ; CHKERRQ ( ierr )
2012-11-12 19:44:39 +05:30
call SNESSetConvergenceTest ( snes , AL_converged , dummy , PETSC_NULL_FUNCTION , ierr )
CHKERRQ ( ierr )
2013-01-09 03:24:25 +05:30
call SNESSetFromOptions ( snes , ierr ) ; CHKERRQ ( ierr )
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! init fields
2012-12-17 15:48:39 +05:30
call DMDAVecGetArrayF90 ( da , solution_vec , xx_psc , ierr ) ; CHKERRQ ( ierr ) ! places pointer xx_psc on PETSc data
2012-11-12 19:44:39 +05:30
F = > xx_psc ( 0 : 8 , : , : , : )
2013-07-24 18:36:16 +05:30
F_lambda = > xx_psc ( 9 : 17 , : , : , : )
2012-11-12 19:44:39 +05:30
if ( restartInc == 1_pInt ) then ! no deformation (no restart)
2013-05-08 21:22:29 +05:30
F_lastInc = spread ( spread ( spread ( math_I3 , 3 , grid ( 1 ) ) , 4 , grid ( 2 ) ) , 5 , grid ( 3 ) ) ! initialize to identity
2013-07-26 21:55:37 +05:30
F_lastInc2 = F_lastInc
2013-07-24 18:36:16 +05:30
F_lambda_lastInc = F_lastInc
2013-05-08 21:22:29 +05:30
F = reshape ( F_lastInc , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
2013-07-24 18:36:16 +05:30
F_lambda = F
2013-05-08 21:22:29 +05:30
elseif ( restartInc > 1_pInt ) then
if ( iand ( debug_level ( debug_spectral ) , debug_spectralRestart ) / = 0 ) &
write ( 6 , '(/,a,' / / IO_intOut ( restartInc - 1_pInt ) / / ',a)' ) &
'reading values of increment' , restartInc - 1_pInt , 'from file'
2012-12-14 23:00:22 +05:30
flush ( 6 )
call IO_read_jobBinaryFile ( 777 , 'F' , &
trim ( getSolverJobName ( ) ) , size ( F ) )
2012-11-12 19:44:39 +05:30
read ( 777 , rec = 1 ) F
close ( 777 )
2012-12-14 23:00:22 +05:30
call IO_read_jobBinaryFile ( 777 , 'F_lastInc' , &
2012-11-12 19:44:39 +05:30
trim ( getSolverJobName ( ) ) , size ( F_lastInc ) )
read ( 777 , rec = 1 ) F_lastInc
close ( 777 )
2013-07-26 21:55:37 +05:30
call IO_read_jobBinaryFile ( 777 , 'F_lastInc2' , &
trim ( getSolverJobName ( ) ) , size ( F_lastInc2 ) )
read ( 777 , rec = 1 ) F_lastInc2
close ( 777 )
2012-12-14 23:00:22 +05:30
F_aim = reshape ( sum ( sum ( sum ( F , dim = 4 ) , dim = 3 ) , dim = 2 ) * wgt , [ 3 , 3 ] ) ! average of F
F_aim_lastInc = sum ( sum ( sum ( F_lastInc , dim = 5 ) , dim = 4 ) , dim = 3 ) * wgt ! average of F_lastInc
2013-07-24 18:36:16 +05:30
call IO_read_jobBinaryFile ( 777 , 'F_lambda' , &
trim ( getSolverJobName ( ) ) , size ( F_lambda ) )
read ( 777 , rec = 1 ) F_lambda
2012-11-12 19:44:39 +05:30
close ( 777 )
2013-07-24 18:36:16 +05:30
call IO_read_jobBinaryFile ( 777 , 'F_lambda_lastInc' , &
trim ( getSolverJobName ( ) ) , size ( F_lambda_lastInc ) )
read ( 777 , rec = 1 ) F_lambda_lastInc
2012-11-12 19:44:39 +05:30
close ( 777 )
2013-01-24 00:03:46 +05:30
call IO_read_jobBinaryFile ( 777 , 'F_aimDot' , trim ( getSolverJobName ( ) ) , size ( f_aimDot ) )
read ( 777 , rec = 1 ) f_aimDot
2012-12-14 23:00:22 +05:30
close ( 777 )
2013-07-23 23:12:15 +05:30
call IO_read_jobBinaryFile ( 777 , 'C_volAvg' , trim ( getSolverJobName ( ) ) , size ( C_volAvg ) )
read ( 777 , rec = 1 ) C_volAvg
2012-12-14 23:00:22 +05:30
close ( 777 )
2013-07-23 23:12:15 +05:30
call IO_read_jobBinaryFile ( 777 , 'C_volAvgLastInc' , trim ( getSolverJobName ( ) ) , size ( C_volAvgLastInc ) )
read ( 777 , rec = 1 ) C_volAvgLastInc
2012-11-12 19:44:39 +05:30
close ( 777 )
2012-12-14 23:00:22 +05:30
call IO_read_jobBinaryFile ( 777 , 'C_ref' , trim ( getSolverJobName ( ) ) , size ( temp3333_Real ) )
2013-07-23 23:12:15 +05:30
read ( 777 , rec = 1 ) C_minMaxAvg
2012-11-12 19:44:39 +05:30
close ( 777 )
endif
2013-05-08 21:22:29 +05:30
mesh_ipCoordinates = reshape ( mesh_deformedCoordsFFT ( geomSize , reshape ( &
F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) ) , [ 3 , 1 , product ( grid ) ] )
2013-03-06 20:01:13 +05:30
call Utilities_constitutiveResponse ( F , F , temperature , 0.0_pReal , P , temp3333_Real , temp3333_Real2 , &
2012-12-14 23:00:22 +05:30
temp33_Real , . false . , math_I3 )
2012-12-17 15:48:39 +05:30
call DMDAVecRestoreArrayF90 ( da , solution_vec , xx_psc , ierr ) ; CHKERRQ ( ierr )
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! reference stiffness
2012-12-14 23:00:22 +05:30
if ( restartInc == 1_pInt ) then ! use initial stiffness as reference stiffness
2013-07-23 23:12:15 +05:30
C_minMaxAvg = temp3333_Real2
C_volAvg = temp3333_Real
2012-12-14 23:00:22 +05:30
endif
2012-11-12 19:44:39 +05:30
2013-03-06 20:01:13 +05:30
call Utilities_updateGamma ( temp3333_Real2 , . True . )
C_scale = temp3333_Real2
S_scale = math_invSym3333 ( temp3333_Real2 )
2012-07-25 19:31:39 +05:30
2012-11-12 19:44:39 +05:30
end subroutine AL_init
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief solution for the AL scheme with internal iterations
!--------------------------------------------------------------------------------------------------
2012-11-12 19:44:39 +05:30
type ( tSolutionState ) function &
2013-07-26 21:55:37 +05:30
AL_solution ( incInfoIn , guess , timeinc , timeinc_old , loadCaseTime , P_BC , F_BC , temperature_bc , rotation_BC , density )
2012-11-12 19:44:39 +05:30
use numerics , only : &
2013-07-30 21:02:55 +05:30
update_gamma
2012-11-12 19:44:39 +05:30
use math , only : &
math_mul33x33 , &
2013-07-23 23:12:15 +05:30
math_mul3333xx33 , &
2013-03-06 20:01:13 +05:30
math_rotate_backward33 , &
math_invSym3333
2013-07-24 18:36:16 +05:30
use mesh , only : &
2013-02-05 20:33:36 +05:30
mesh_ipCoordinates , &
mesh_deformedCoordsFFT
2012-11-12 19:44:39 +05:30
use IO , only : &
IO_write_JobBinaryFile
use DAMASK_spectral_Utilities , only : &
2013-05-08 21:22:29 +05:30
grid , &
geomSize , &
2012-11-12 19:44:39 +05:30
tBoundaryCondition , &
Utilities_forwardField , &
Utilities_calculateRate , &
Utilities_maskedCompliance , &
Utilities_updateGamma , &
cutBack
use FEsolving , only : &
restartWrite , &
terminallyIll
implicit none
2012-11-07 18:41:41 +05:30
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90>
2013-05-13 15:14:23 +05:30
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
! input data for solution
2013-05-13 15:14:23 +05:30
real ( pReal ) , intent ( in ) :: &
timeinc , & !< increment in time for current solution
timeinc_old , & !< increment in time of last increment
loadCaseTime , & !< remaining time of current load case
2013-07-30 21:02:55 +05:30
temperature_bc , &
density
2012-12-17 15:48:39 +05:30
logical , intent ( in ) :: &
guess
type ( tBoundaryCondition ) , intent ( in ) :: &
P_BC , &
F_BC
character ( len = * ) , intent ( in ) :: &
incInfoIn
2012-11-12 19:44:39 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: rotation_BC
2013-07-30 21:02:55 +05:30
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
2012-11-12 19:44:39 +05:30
! PETSc Data
2013-07-24 18:36:16 +05:30
PetscScalar , dimension ( : , : , : , : ) , pointer :: xx_psc , F , F_lambda
2012-11-12 19:44:39 +05:30
PetscErrorCode :: ierr
2013-07-24 18:36:16 +05:30
SNESConvergedReason :: reason
2012-08-10 22:31:58 +05:30
2013-07-24 18:36:16 +05:30
incInfo = incInfoIn
2012-12-14 23:00:22 +05:30
call DMDAVecGetArrayF90 ( da , solution_vec , xx_psc , ierr )
F = > xx_psc ( 0 : 8 , : , : , : )
2013-07-24 18:36:16 +05:30
F_lambda = > xx_psc ( 9 : 17 , : , : , : )
2012-12-14 23:00:22 +05:30
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
2012-11-12 19:44:39 +05:30
if ( restartWrite ) then
2013-01-10 03:49:32 +05:30
write ( 6 , '(/,a)' ) ' writing converged results for restart'
flush ( 6 )
2012-12-14 23:00:22 +05:30
call IO_write_jobBinaryFile ( 777 , 'F' , size ( F ) ) ! writing deformation gradient field to file
write ( 777 , rec = 1 ) F
close ( 777 )
call IO_write_jobBinaryFile ( 777 , 'F_lastInc' , size ( F_lastInc ) ) ! writing F_lastInc field to file
write ( 777 , rec = 1 ) F_lastInc
close ( 777 )
2013-07-24 18:36:16 +05:30
call IO_write_jobBinaryFile ( 777 , 'F_lambda' , size ( F_lambda ) ) ! writing deformation gradient field to file
write ( 777 , rec = 1 ) F_lambda
2012-12-14 23:00:22 +05:30
close ( 777 )
2013-07-24 18:36:16 +05:30
call IO_write_jobBinaryFile ( 777 , 'F_lambda_lastInc' , size ( F_lambda_lastInc ) ) ! writing F_lastInc field to file
write ( 777 , rec = 1 ) F_lambda_lastInc
2012-11-12 19:44:39 +05:30
close ( 777 )
2013-01-18 17:00:52 +05:30
call IO_write_jobBinaryFile ( 777 , 'F_aimDot' , size ( F_aimDot ) )
write ( 777 , rec = 1 ) F_aimDot
close ( 777 )
2013-07-23 23:12:15 +05:30
call IO_write_jobBinaryFile ( 777 , 'C_volAvg' , size ( C_volAvg ) )
write ( 777 , rec = 1 ) C_volAvg
2012-11-12 19:44:39 +05:30
close ( 777 )
2013-07-23 23:12:15 +05:30
call IO_write_jobBinaryFile ( 777 , 'C_volAvgLastInc' , size ( C_volAvgLastInc ) )
write ( 777 , rec = 1 ) C_volAvgLastInc
2012-12-14 23:00:22 +05:30
close ( 777 )
2012-11-12 19:44:39 +05:30
endif
AL_solution % converged = . false .
2013-05-13 15:14:23 +05:30
if ( cutBack ) then
2012-11-12 19:44:39 +05:30
F_aim = F_aim_lastInc
2013-07-24 18:36:16 +05:30
F_lambda = reshape ( F_lambda_lastInc , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
2013-05-08 21:22:29 +05:30
F = reshape ( F_lastInc , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
2013-07-23 23:12:15 +05:30
C_volAvg = C_volAvgLastInc
2012-11-12 19:44:39 +05:30
else
2013-07-23 23:12:15 +05:30
C_volAvgLastInc = C_volAvg
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
2012-10-02 20:56:56 +05:30
! calculate rate for aim
if ( F_BC % myType == 'l' ) then ! calculate f_aimDot from given L and current F
f_aimDot = F_BC % maskFloat * math_mul33x33 ( F_BC % values , F_aim )
2013-05-13 15:14:23 +05:30
elseif ( F_BC % myType == 'fdot' ) then ! f_aimDot is prescribed
2012-10-02 20:56:56 +05:30
f_aimDot = F_BC % maskFloat * F_BC % values
2013-05-13 15:14:23 +05:30
elseif ( F_BC % myType == 'f' ) then ! aim at end of load case is prescribed
f_aimDot = F_BC % maskFloat * ( F_BC % values - F_aim ) / loadCaseTime
2012-08-06 14:23:12 +05:30
endif
2012-11-12 19:44:39 +05:30
if ( guess ) f_aimDot = f_aimDot + P_BC % maskFloat * ( F_aim - F_aim_lastInc ) / timeinc_old
2012-10-02 20:56:56 +05:30
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
2013-05-08 21:22:29 +05:30
mesh_ipCoordinates = reshape ( mesh_deformedCoordsFFT ( geomSize , reshape ( &
F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) ) , [ 3 , 1 , product ( grid ) ] )
2012-10-02 20:56:56 +05:30
Fdot = Utilities_calculateRate ( math_rotate_backward33 ( f_aimDot , rotation_BC ) , &
2013-05-08 21:22:29 +05:30
timeinc_old , guess , F_lastInc , reshape ( F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) )
2013-07-24 18:36:16 +05:30
F_lambdaDot = Utilities_calculateRate ( math_rotate_backward33 ( f_aimDot , rotation_BC ) , &
timeinc_old , guess , F_lambda_lastInc , reshape ( F_lambda , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) )
2013-07-26 21:55:37 +05:30
F_lastInc2 = F_lastInc
2013-05-08 21:22:29 +05:30
F_lastInc = reshape ( F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
2013-07-24 18:36:16 +05:30
F_lambda_lastInc = reshape ( F_lambda , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
2012-10-02 20:56:56 +05:30
endif
F_aim = F_aim + f_aimDot * timeinc
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
2012-12-17 15:48:39 +05:30
! update local deformation gradient
2013-03-31 18:36:49 +05:30
F = reshape ( Utilities_forwardField ( timeinc , F_lastInc , Fdot , & ! ensure that it matches rotated F_aim
2013-05-08 21:22:29 +05:30
math_rotate_backward33 ( F_aim , rotation_BC ) ) , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
2013-07-24 18:36:16 +05:30
F_lambda = reshape ( Utilities_forwardField ( timeinc , F_lambda_lastInc , F_lambdadot ) , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) ! does not have any average value as boundary condition
2012-11-12 19:44:39 +05:30
call DMDAVecRestoreArrayF90 ( da , solution_vec , xx_psc , ierr )
CHKERRQ ( ierr )
2012-10-02 20:56:56 +05:30
2012-08-06 14:23:12 +05:30
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
2013-07-23 23:12:15 +05:30
S = Utilities_maskedCompliance ( rotation_BC , P_BC % maskLogical , C_volAvg )
2013-03-06 20:01:13 +05:30
if ( update_gamma ) then
2013-07-23 23:12:15 +05:30
call Utilities_updateGamma ( C_minMaxAvg , restartWrite )
C_scale = C_minMaxAvg
S_scale = math_invSym3333 ( C_minMaxAvg )
2013-03-06 20:01:13 +05:30
endif
2012-11-12 19:44:39 +05:30
ForwardData = . True .
mask_stress = P_BC % maskFloat
params % P_BC = P_BC % values
params % rotation_BC = rotation_BC
params % timeinc = timeinc
2013-07-26 21:55:37 +05:30
params % timeincOld = timeinc_old
2013-01-09 23:38:08 +05:30
params % temperature = temperature_bc
2013-07-26 21:55:37 +05:30
params % density = density
2012-11-12 19:44:39 +05:30
2013-07-23 23:12:15 +05:30
!--------------------------------------------------------------------------------------------------
! solve BVP
2013-07-24 18:36:16 +05:30
call SNESSolve ( snes , PETSC_NULL_OBJECT , solution_vec , ierr )
CHKERRQ ( ierr )
2013-07-23 23:12:15 +05:30
!--------------------------------------------------------------------------------------------------
! check convergence
2013-07-24 18:36:16 +05:30
call SNESGetConvergedReason ( snes , reason , ierr )
CHKERRQ ( ierr )
AL_solution % termIll = terminallyIll
terminallyIll = . false .
2013-07-30 21:02:55 +05:30
AL_solution % converged = . true .
if ( reason < 1 ) AL_solution % converged = . false .
AL_solution % iterationsNeeded = totalIter
2012-11-12 19:44:39 +05:30
end function AL_solution
2012-10-02 20:56:56 +05:30
2012-08-03 14:55:48 +05:30
2012-08-06 22:57:53 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief forms the AL residual vector
!--------------------------------------------------------------------------------------------------
2012-11-12 19:44:39 +05:30
subroutine AL_formResidual ( in , x_scal , f_scal , dummy , ierr )
use numerics , only : &
itmax , &
2013-02-28 23:05:02 +05:30
itmin , &
polarAlpha , &
2013-07-30 21:02:55 +05:30
polarBeta , &
err_stress_tolrel , &
err_stress_tolabs , &
err_f_tol , &
err_p_tol , &
err_stress_tolabs
2013-05-08 21:22:29 +05:30
use IO , only : &
IO_intOut
2012-11-12 19:44:39 +05:30
use math , only : &
math_rotate_backward33 , &
math_transpose33 , &
2013-03-04 15:19:40 +05:30
math_mul3333xx33 , &
2013-07-30 21:02:55 +05:30
math_invSym3333 , &
PI
2012-11-12 19:44:39 +05:30
use DAMASK_spectral_Utilities , only : &
2013-05-08 21:22:29 +05:30
grid , &
2013-07-26 21:55:37 +05:30
geomSize , &
2013-05-08 21:22:29 +05:30
wgt , &
2012-11-12 19:44:39 +05:30
field_real , &
2013-07-26 21:55:37 +05:30
field_fourier , &
2012-11-12 19:44:39 +05:30
Utilities_FFTforward , &
Utilities_fourierConvolution , &
2013-07-26 21:55:37 +05:30
Utilities_inverseLaplace , &
2012-11-12 19:44:39 +05:30
Utilities_FFTbackward , &
2013-05-08 21:22:29 +05:30
Utilities_constitutiveResponse
use debug , only : &
debug_level , &
debug_spectral , &
debug_spectralRotation
2013-03-04 15:19:40 +05:30
use homogenization , only : &
materialpoint_P , &
materialpoint_dPdF
2012-11-12 19:44:39 +05:30
implicit none
2013-01-10 21:06:55 +05:30
!--------------------------------------------------------------------------------------------------
! strange syntax in the next line because otherwise macros expand beyond 132 character limit
DMDALocalInfo , dimension ( &
DMDA_LOCAL_INFO_SIZE ) :: &
2013-01-10 19:03:43 +05:30
in
2013-01-10 21:06:55 +05:30
PetscScalar , target , dimension ( 3 , 3 , 2 , &
XG_RANGE , YG_RANGE , ZG_RANGE ) :: &
2013-01-10 19:03:43 +05:30
x_scal
2013-01-10 21:06:55 +05:30
PetscScalar , target , dimension ( 3 , 3 , 2 , &
X_RANGE , Y_RANGE , Z_RANGE ) :: &
2013-01-10 19:03:43 +05:30
f_scal
2013-01-10 21:06:55 +05:30
PetscScalar , pointer , dimension ( : , : , : , : , : ) :: &
2013-01-10 19:03:43 +05:30
F , &
2013-07-24 18:36:16 +05:30
F_lambda , &
2013-01-10 21:06:55 +05:30
residual_F , &
2013-07-24 18:36:16 +05:30
residual_F_lambda
2013-01-10 19:03:43 +05:30
PetscInt :: &
2013-07-23 23:12:15 +05:30
PETScIter , &
2013-01-10 19:03:43 +05:30
nfuncs
2012-11-12 19:44:39 +05:30
PetscObject :: dummy
PetscErrorCode :: ierr
2013-01-11 00:20:14 +05:30
integer ( pInt ) :: &
2013-07-24 18:36:16 +05:30
i , j , k
2013-07-30 21:02:55 +05:30
real ( pReal ) :: correctionFactor
2012-11-12 19:44:39 +05:30
2013-03-31 18:36:49 +05:30
F = > x_scal ( 1 : 3 , 1 : 3 , 1 , &
2013-01-10 21:06:55 +05:30
XG_RANGE , YG_RANGE , ZG_RANGE )
2013-07-24 18:36:16 +05:30
F_lambda = > x_scal ( 1 : 3 , 1 : 3 , 2 , &
2013-01-10 21:06:55 +05:30
XG_RANGE , YG_RANGE , ZG_RANGE )
2013-03-31 18:36:49 +05:30
residual_F = > f_scal ( 1 : 3 , 1 : 3 , 1 , &
2013-01-10 21:06:55 +05:30
X_RANGE , Y_RANGE , Z_RANGE )
2013-07-24 18:36:16 +05:30
residual_F_lambda = > f_scal ( 1 : 3 , 1 : 3 , 2 , &
2013-01-10 21:06:55 +05:30
X_RANGE , Y_RANGE , Z_RANGE )
2012-11-12 19:44:39 +05:30
2013-01-10 03:49:32 +05:30
call SNESGetNumberFunctionEvals ( snes , nfuncs , ierr ) ; CHKERRQ ( ierr )
2013-07-23 23:12:15 +05:30
call SNESGetIterationNumber ( snes , PETScIter , ierr ) ; CHKERRQ ( ierr )
2012-11-12 19:44:39 +05:30
2013-07-30 21:02:55 +05:30
if ( nfuncs == 0 . and . PETScIter == 0 ) totalIter = - 1_pInt ! new increment
if ( totalIter < = PETScIter ) then ! new iteration
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
2013-07-23 23:12:15 +05:30
totalIter = totalIter + 1_pInt
2013-01-10 03:49:32 +05:30
write ( 6 , '(1x,a,3(a,' / / IO_intOut ( itmax ) / / '))' ) trim ( incInfo ) , &
2013-07-23 23:12:15 +05:30
' @ Iteration ' , itmin , '≤' , totalIter , '≤' , itmax
2013-05-08 21:22:29 +05:30
if ( iand ( debug_level ( debug_spectral ) , debug_spectralRotation ) / = 0 ) &
2013-01-10 03:49:32 +05:30
write ( 6 , '(/,a,/,3(3(f12.7,1x)/))' , advance = 'no' ) ' deformation gradient aim (lab) =' , &
2012-12-17 15:48:39 +05:30
math_transpose33 ( math_rotate_backward33 ( F_aim , params % rotation_BC ) )
2013-01-10 03:49:32 +05:30
write ( 6 , '(/,a,/,3(3(f12.7,1x)/))' , advance = 'no' ) ' deformation gradient aim =' , &
2012-12-17 15:48:39 +05:30
math_transpose33 ( F_aim )
2013-01-10 03:49:32 +05:30
flush ( 6 )
2012-10-02 20:56:56 +05:30
endif
2013-07-26 21:55:37 +05:30
!--------------------------------------------------------------------------------------------------
! evaluate inertia
2013-07-30 21:02:55 +05:30
dynamic : if ( params % density > 0.0_pReal ) then
residual_F = ( ( F - F_lastInc ) / params % timeinc - ( F_lastInc - F_lastInc2 ) / params % timeincOld ) / &
( ( params % timeinc + params % timeincOld ) / 2.0_pReal )
2013-07-26 21:55:37 +05:30
residual_F = params % density * product ( geomSize / grid ) * residual_F
field_real = 0.0_pReal
field_real ( 1 : grid ( 1 ) , 1 : grid ( 2 ) , 1 : grid ( 3 ) , 1 : 3 , 1 : 3 ) = reshape ( residual_F , [ grid ( 1 ) , grid ( 2 ) , grid ( 3 ) , 3 , 3 ] , &
order = [ 4 , 5 , 1 , 2 , 3 ] ) ! field real has a different order
call Utilities_FFTforward ( )
call Utilities_inverseLaplace ( )
inertiaField_fourier = field_fourier
2013-07-30 21:02:55 +05:30
else dynamic
2013-07-26 21:55:37 +05:30
inertiaField_fourier = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2013-07-30 21:02:55 +05:30
endif dynamic
2013-07-26 21:55:37 +05:30
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
!
field_real = 0.0_pReal
2013-05-08 21:22:29 +05:30
do k = 1_pInt , grid ( 3 ) ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid ( 1 )
2013-07-24 18:36:16 +05:30
field_real ( i , j , k , 1 : 3 , 1 : 3 ) = math_mul3333xx33 ( C_scale , polarBeta * F ( 1 : 3 , 1 : 3 , i , j , k ) - &
polarAlpha * F_lambda ( 1 : 3 , 1 : 3 , i , j , k ) )
2012-11-12 19:44:39 +05:30
enddo ; enddo ; enddo
!--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space
call Utilities_FFTforward ( )
2013-07-26 21:55:37 +05:30
field_fourier = field_fourier + polarAlpha * inertiaField_fourier
2013-02-28 23:05:02 +05:30
call Utilities_fourierConvolution ( math_rotate_backward33 ( polarBeta * F_aim , params % rotation_BC ) )
2012-11-12 19:44:39 +05:30
call Utilities_FFTbackward ( )
!--------------------------------------------------------------------------------------------------
! constructing residual
2013-07-24 18:36:16 +05:30
residual_F_lambda = polarBeta * F - reshape ( field_real ( 1 : grid ( 1 ) , 1 : grid ( 2 ) , 1 : grid ( 3 ) , 1 : 3 , 1 : 3 ) , &
2013-05-08 21:22:29 +05:30
[ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] , order = [ 3 , 4 , 5 , 1 , 2 ] )
2013-02-28 23:05:02 +05:30
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
2013-07-31 02:34:41 +05:30
P_avLastEval = P_av
2013-07-24 18:36:16 +05:30
call Utilities_constitutiveResponse ( F_lastInc , F - residual_F_lambda / polarBeta , params % temperature , params % timeinc , &
2013-07-23 23:12:15 +05:30
residual_F , C_volAvg , C_minMaxAvg , P_av , ForwardData , params % rotation_BC )
2013-02-28 23:05:02 +05:30
ForwardData = . False .
2013-07-30 21:02:55 +05:30
!--------------------------------------------------------------------------------------------------
! stress BC handling
write ( 6 , '(/,a)' ) ' ... correcting F to fullfill stress BC ....................................'
2013-07-31 02:34:41 +05:30
correctionFactor = ( cos ( ( 1.0 - 1000 0.0_pReal ** ( - sum ( ( P_av - P_avLastEval ) ** 2.0_pReal ) / & ! only correct when averages stress of last two calls doesn't strongly deviate
2013-07-30 21:02:55 +05:30
sum ( P_av ** 2.0_pReal ) ) ) * PI ) + 1.0 ) / 2.0_pReal
2013-07-31 02:34:41 +05:30
write ( 6 , '(/,a,f10.4)' ) ' stress BC correction factor = ' , correctionFactor
2013-07-30 21:02:55 +05:30
F_aim = F_aim - correctionFactor * math_mul3333xx33 ( S , ( ( P_av - params % P_BC ) ) ) ! S = 0.0 for no bc
err_stress = maxval ( abs ( mask_stress * ( P_av - params % P_BC ) ) ) ! mask = 0.0 for no bc
2013-07-24 18:36:16 +05:30
2013-02-28 23:05:02 +05:30
!--------------------------------------------------------------------------------------------------
! constructing residual
2013-03-04 15:19:40 +05:30
err_p = 0.0_pReal
2013-05-08 21:22:29 +05:30
do k = 1_pInt , grid ( 3 ) ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid ( 1 )
2013-07-30 21:02:55 +05:30
err_p = err_p + sum ( ( math_mul3333xx33 ( C_scale , F_lambda ( 1 : 3 , 1 : 3 , i , j , k ) - math_I3 ) - residual_F ( 1 : 3 , 1 : 3 , i , j , k ) ) ** 2.0_pReal )
2013-07-24 18:36:16 +05:30
residual_F ( 1 : 3 , 1 : 3 , i , j , k ) = math_I3 + math_mul3333xx33 ( S_scale , residual_F ( 1 : 3 , 1 : 3 , i , j , k ) ) - &
F_lambda ( 1 : 3 , 1 : 3 , i , j , k ) &
+ residual_F_lambda ( 1 : 3 , 1 : 3 , i , j , k )
2013-02-28 23:05:02 +05:30
enddo ; enddo ; enddo
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! calculating errors
2013-07-24 18:36:16 +05:30
err_f = wgt * sqrt ( sum ( residual_F_lambda ** 2.0_pReal ) ) / polarBeta
2013-03-04 15:19:40 +05:30
err_p = wgt * sqrt ( err_p )
2012-08-09 18:34:56 +05:30
2012-11-12 19:44:39 +05:30
end subroutine AL_formResidual
2012-08-06 14:23:12 +05:30
2012-08-06 22:57:53 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
2013-07-23 23:12:15 +05:30
subroutine AL_converged ( snes_local , PETScIter , xnorm , snorm , fnorm , reason , dummy , ierr )
2013-07-30 21:02:55 +05:30
use numerics , only : &
2012-11-12 19:44:39 +05:30
itmax , &
itmin , &
err_f_tol , &
err_p_tol , &
2013-07-30 21:02:55 +05:30
err_stress_tolabs , &
err_stress_tolrel
2013-03-06 20:01:13 +05:30
use FEsolving , only : &
terminallyIll
2012-11-12 19:44:39 +05:30
2013-01-10 03:49:32 +05:30
implicit none
2012-11-12 19:44:39 +05:30
SNES :: snes_local
2013-07-23 23:12:15 +05:30
PetscInt :: PETScIter
2013-01-10 03:49:32 +05:30
PetscReal :: &
xnorm , &
snorm , &
fnorm
2012-11-12 19:44:39 +05:30
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
2013-07-30 21:02:55 +05:30
real ( pReal ) :: &
mismatch_f_tol , &
mismatch_p_tol , &
stressBC_tol
mismatch_f_tol = max ( maxval ( abs ( F_aim - math_I3 ) ) * err_stress_tolrel , err_f_tol )
mismatch_p_tol = max ( maxval ( abs ( P_av ) ) * err_stress_tolrel , err_p_tol )
stressBC_tol = max ( maxval ( abs ( P_av ) ) * err_stress_tolrel , err_stress_tolabs )
2013-07-23 23:12:15 +05:30
write ( 6 , '(1/,a)' ) ' ... reporting .............................................................'
2013-01-10 03:49:32 +05:30
write ( 6 , '(/,a,f8.2,a,es11.5,a,es11.4,a)' ) ' mismatch F = ' , &
2013-07-30 21:02:55 +05:30
err_f / mismatch_f_tol , &
' (' , err_f , ' -, tol =' , mismatch_f_tol , ')'
2013-01-10 03:49:32 +05:30
write ( 6 , '(a,f8.2,a,es11.5,a,es11.4,a)' ) ' mismatch P = ' , &
2013-07-30 21:02:55 +05:30
err_p / mismatch_p_tol , &
' (' , err_p , ' -, tol =' , mismatch_p_tol , ')'
write ( 6 , '(a,f8.2,a,es11.5,a,es11.4,a)' ) ' error stress BC = ' , &
err_stress / stressBC_tol , ' (' , err_stress , ' Pa, tol =' , stressBC_tol , ')'
write ( 6 , '(/,a)' ) ' ==========================================================================='
flush ( 6 )
2013-07-24 18:36:16 +05:30
converged : if ( ( totalIter > = itmin . and . &
2013-07-30 21:02:55 +05:30
all ( [ err_f / mismatch_f_tol , &
err_p / mismatch_p_tol , &
err_stress / stressBC_tol ] < 1.0_pReal ) ) &
2013-07-24 18:36:16 +05:30
. or . terminallyIll ) then
reason = 1
2013-07-30 21:02:55 +05:30
elseif ( totalIter > = itmax ) then converged
2013-07-24 18:36:16 +05:30
reason = - 1
else converged
reason = 0
endif converged
2013-01-10 03:49:32 +05:30
2012-11-12 19:44:39 +05:30
end subroutine AL_converged
2012-08-03 14:55:48 +05:30
2012-08-06 22:57:53 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
2012-11-12 19:44:39 +05:30
subroutine AL_destroy ( )
use DAMASK_spectral_Utilities , only : &
Utilities_destroy
2013-01-10 03:49:32 +05:30
2012-11-12 19:44:39 +05:30
implicit none
PetscErrorCode :: ierr
2013-01-10 03:49:32 +05:30
call VecDestroy ( solution_vec , ierr ) ; CHKERRQ ( ierr )
call SNESDestroy ( snes , ierr ) ; CHKERRQ ( ierr )
call DMDestroy ( da , ierr ) ; CHKERRQ ( ierr )
call PetscFinalize ( ierr ) ; CHKERRQ ( ierr )
2012-11-12 19:44:39 +05:30
call Utilities_destroy ( )
2012-07-26 19:28:47 +05:30
2012-11-12 19:44:39 +05:30
end subroutine AL_destroy
2012-07-26 19:28:47 +05:30
2012-11-12 19:44:39 +05:30
end module DAMASK_spectral_SolverAL