Merge branch '38-introduce-rudimentary-PETSc-based-FEM-solver' of magit1.mpie.de:damask/DAMASK into 38-introduce-rudimentary-PETSc-based-FEM-solver

This commit is contained in:
Jaeyong Jung 2018-09-21 10:13:27 +02:00
commit 9150844c96
11 changed files with 6 additions and 218 deletions

View File

@ -57,22 +57,9 @@ program DAMASK_FEM
maxFields, &
nActiveFields, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_SOLUTE_ID, &
FIELD_MGTWIN_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID, &
COMPONENT_DAMAGE_PHI_ID, &
COMPONENT_SOLUTE_CV_ID, &
COMPONENT_SOLUTE_CVPOT_ID, &
COMPONENT_SOLUTE_CH_ID, &
COMPONENT_SOLUTE_CHPOT_ID, &
COMPONENT_SOLUTE_CVaH_ID, &
COMPONENT_SOLUTE_CVaHPOT_ID, &
COMPONENT_MGTWIN_PHI_ID, &
FIELD_MECH_label
use FEM_mech
@ -131,11 +118,8 @@ program DAMASK_FEM
external :: &
MPI_abort, &
DMGetDimension, &
DMGetLabelSize, &
DMGetLabelIdIS, &
ISDestroy, &
quit
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
@ -264,124 +248,6 @@ program DAMASK_FEM
enddo
endif
enddo
case('temp','temperature') ! thermal field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_THERMAL_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_THERMAL_T_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('mgtwin') ! mgtwin field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MGTWIN_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MGTWIN_PHI_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('damage')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_DAMAGE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_DAMAGE_PHI_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cv')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CV_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('ch')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CH_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('chpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CHPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvah')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaH_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('cvahpot')
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaHPOT_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
end select
enddo; enddo
close(FILEUNIT)

View File

@ -88,9 +88,7 @@ subroutine DAMASK_interface_init()
dateAndTime ! type default integer
PetscErrorCode :: ierr
external :: &
quit,&
PETScErrorF, & ! is called in the CHKERRQ macro
PETScInitialize
quit
open(6, encoding='UTF-8') ! for special characters in output

View File

@ -151,7 +151,6 @@ program DAMASK_spectral
external :: &
quit
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
@ -659,9 +658,6 @@ subroutine quit(stop_id)
PetscErrorCode :: ierr = 0
logical :: ErrorInQuit
external :: &
PETScFinalize
call PETScFinalize(ierr)
if (ierr /= 0) write(6,'(a)') ' Error in PETScFinalize'
#ifdef _OPENMP

View File

@ -62,7 +62,6 @@ module FEM_mech
FEM_mech_solution ,&
FEM_mech_forward, &
FEM_mech_destroy
external :: PETScerrorf
contains
!--------------------------------------------------------------------------------------------------

View File

@ -4,7 +4,8 @@
!--------------------------------------------------------------------------------------------------
module FEM_utilities
#include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petsc.h>
#include <petsc/finclude/petscdmda.h>
#include <petsc/finclude/petscis.h>
use prec, only: pReal, pInt
use PETScdmplex
@ -32,8 +33,6 @@ use PETScis
character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical'
integer(pInt), parameter :: structOrder = 2_pInt
enum, bind(c)
enumerator :: FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
@ -122,7 +121,6 @@ use PETScis
COMPONENT_MECH_Z_ID, &
COMPONENT_THERMAL_T_ID
external :: PETScErrorF
contains
!--------------------------------------------------------------------------------------------------
@ -138,6 +136,7 @@ subroutine utilities_init()
IO_timeStamp, &
IO_open_file
use numerics, only: &
structOrder, &
integrationOrder, &
worldsize, &
worldrank, &

View File

@ -127,12 +127,7 @@ module numerics
#ifdef FEM
integer(pInt), protected, public :: &
integrationOrder = 2_pInt, & !< order of quadrature rule required
structOrder = 2_pInt, & !< order of displacement shape functions
thermalOrder = 2_pInt, & !< order of temperature field shape functions
damageOrder = 2_pInt, & !< order of damage field shape functions
vacancyfluxOrder = 2_pInt, & !< order of vacancy concentration and chemical potential field shape functions
porosityOrder = 2_pInt, & !< order of porosity field shape functions
hydrogenfluxOrder = 2_pInt !< order of hydrogen concentration and chemical potential field shape functions
structOrder = 2_pInt !< order of displacement shape functions
logical, protected, public :: &
BBarStabilisation = .false.
character(len=4096), protected, public :: &
@ -197,8 +192,6 @@ subroutine numerics_init
tag ,&
line
!$ character(len=6) DAMASK_NumThreadsString ! environment variable DAMASK_NUM_THREADS
external :: &
PETScErrorF ! is called in the CHKERRQ macro
#ifdef PETSc
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
@ -425,16 +418,6 @@ subroutine numerics_init
integrationorder = IO_intValue(line,chunkPos,2_pInt)
case ('structorder')
structorder = IO_intValue(line,chunkPos,2_pInt)
case ('thermalorder')
thermalorder = IO_intValue(line,chunkPos,2_pInt)
case ('damageorder')
damageorder = IO_intValue(line,chunkPos,2_pInt)
case ('vacancyfluxorder')
vacancyfluxOrder = IO_intValue(line,chunkPos,2_pInt)
case ('porosityorder')
porosityOrder = IO_intValue(line,chunkPos,2_pInt)
case ('hydrogenfluxorder')
hydrogenfluxOrder = IO_intValue(line,chunkPos,2_pInt)
case ('petsc_options')
petsc_options = trim(line(chunkPos(4):))
case ('bbarstabilisation')
@ -587,11 +570,6 @@ subroutine numerics_init
#ifdef FEM
write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder
write(6,'(a24,1x,i8)') ' structOrder: ',structOrder
write(6,'(a24,1x,i8)') ' thermalOrder: ',thermalOrder
write(6,'(a24,1x,i8)') ' damageOrder: ',damageOrder
write(6,'(a24,1x,i8)') ' vacancyfluxOrder: ',vacancyfluxOrder
write(6,'(a24,1x,i8)') ' porosityOrder: ',porosityOrder
write(6,'(a24,1x,i8)') ' hydrogenfluxOrder: ',hydrogenfluxOrder
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_defaultOptions)//' '//trim(petsc_options)
write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation
#endif

View File

@ -50,8 +50,6 @@ module spectral_damage
spectral_damage_init, &
spectral_damage_solution, &
spectral_damage_forward
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -85,11 +83,6 @@ subroutine spectral_damage_init()
Vec :: uBound, lBound
PetscErrorCode :: ierr
character(len=100) :: snes_type
external :: &
SNESSetOptionsPrefix, &
SNESGetType, &
DMDAGetCorners, &
DMDASNESSetFunctionLocal
write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press, '
@ -194,11 +187,6 @@ type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadC
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
VecMin, &
VecMax, &
SNESSolve
spectral_damage_solution%converged =.false.
!--------------------------------------------------------------------------------------------------

View File

@ -66,8 +66,6 @@ module spectral_mech_basic
basic_init, &
basic_solution, &
basic_forward
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -119,11 +117,6 @@ subroutine basic_init
integer(pInt) :: proc
character(len=1024) :: rankStr
external :: &
SNESSetOptionsPrefix, &
SNESSetConvergenceTest, &
DMDASNESSetFunctionLocal
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:3145, 2015'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
@ -246,9 +239,6 @@ type(tSolutionState) function basic_solution(incInfoIn,timeinc,timeinc_old,stres
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
SNESsolve
incInfo = incInfoIn
!--------------------------------------------------------------------------------------------------

View File

@ -73,8 +73,6 @@ module spectral_mech_Polarisation
Polarisation_init, &
Polarisation_solution, &
Polarisation_forward
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -130,11 +128,6 @@ subroutine Polarisation_init
integer(pInt) :: proc
character(len=1024) :: rankStr
external :: &
SNESSetOptionsPrefix, &
SNESSetConvergenceTest, &
DMDASNESsetFunctionLocal
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:3145, 2015'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
@ -272,9 +265,6 @@ type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_ol
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
SNESSolve
incInfo = incInfoIn
!--------------------------------------------------------------------------------------------------

View File

@ -50,8 +50,6 @@ module spectral_thermal
spectral_thermal_init, &
spectral_thermal_solution, &
spectral_thermal_forward
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -88,11 +86,6 @@ subroutine spectral_thermal_init
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
external :: &
SNESsetOptionsPrefix, &
DMDAgetCorners, &
DMDASNESsetFunctionLocal
write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press,'
write(6,'(/,a)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018'
@ -196,11 +189,6 @@ type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,load
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
VecMin, &
VecMax, &
SNESSolve
spectral_thermal_solution%converged =.false.
!--------------------------------------------------------------------------------------------------

View File

@ -146,8 +146,6 @@ module spectral_utilities
FIELD_DAMAGE_ID
private :: &
utilities_getFreqDerivative
external :: &
PETScErrorF ! is called in the CHKERRQ macro
contains
@ -209,8 +207,6 @@ subroutine utilities_init()
scalarSize = 1_C_INTPTR_T, &
vecSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T
external :: &
PetscOptionsInsertString
write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>'
write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity, 46:3753, 2013'