From 639ca891334287ece846c32e65f395496da4aa2f Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Fri, 6 Jun 2014 00:38:29 +0000 Subject: [PATCH] DAMASK (except nonlocal) now sees and interacts with FEM solver --- code/CPFEM.f90 | 12 +++++ code/FEsolving.f90 | 11 ++++ code/IO.f90 | 4 ++ code/constitutive.f90 | 7 ++- code/crystallite.f90 | 12 ++++- code/libs.f90 | 3 ++ code/numerics.f90 | 121 ++++++++++++++++++++++++++++++++++++++++-- code/prec.f90 | 3 ++ 8 files changed, 168 insertions(+), 5 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index d03be5fea..77b416736 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -71,6 +71,10 @@ subroutine CPFEM_initAll(temperature,el,ip) use IO, only: & IO_init use DAMASK_interface +#ifdef FEM + use FEZoo, only: & + FEZoo_init +#endif implicit none integer(pInt), intent(in) :: el, & ! FE el number @@ -81,9 +85,15 @@ subroutine CPFEM_initAll(temperature,el,ip) if (.not. CPFEM_init_done) then #ifdef Spectral call DAMASK_interface_init() ! Spectral solver is interfacing to commandline +#endif +#ifdef FEM + call DAMASK_interface_init() ! Spectral solver is interfacing to commandline #endif call prec_init call IO_init +#ifdef FEM + call FEZoo_init() ! Spectral solver is interfacing to commandline +#endif call numerics_init call debug_init call math_init @@ -96,7 +106,9 @@ subroutine CPFEM_initAll(temperature,el,ip) call homogenization_init call CPFEM_init #ifndef Spectral +#ifndef FEM call DAMASK_interface_init() ! Spectral solver init is already done +#endif #endif CPFEM_init_done = .true. endif diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index 52e3a21eb..346c6c5dc 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -67,8 +67,10 @@ subroutine FE_init IO_intValue, & IO_lc, & #ifndef Spectral +#ifndef FEM IO_open_inputFile, & IO_open_logFile, & +#endif #endif IO_warning, & IO_timeStamp @@ -76,6 +78,7 @@ subroutine FE_init implicit none #ifndef Spectral +#ifndef FEM integer(pInt), parameter :: & FILEUNIT = 222_pInt, & maxNchunks = 6_pInt @@ -83,6 +86,7 @@ subroutine FE_init character(len=64) :: tag character(len=1024) :: line integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions +#endif #endif write(6,'(/,a)') ' <<<+- FEsolving init -+>>>' @@ -98,6 +102,13 @@ subroutine FE_init restartInc = 1_pInt endif restartRead = restartInc > 1_pInt ! only read in if "true" restart requested +#elif defined FEM + restartInc = FEMRestartInc + if(restartInc <= 0_pInt) then + call IO_warning(warning_ID=34_pInt) + restartInc = 1_pInt + endif + restartRead = restartInc > 1_pInt #else call IO_open_inputFile(FILEUNIT,modelName) rewind(FILEUNIT) diff --git a/code/IO.f90 b/code/IO.f90 index 6e0fd1d85..849163bbb 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -73,10 +73,12 @@ module IO IO_intOut, & IO_timeStamp #ifndef Spectral +#ifndef FEM public :: & IO_open_inputFile, & IO_open_logFile #endif +#endif #ifdef Abaqus public :: & IO_abaqus_hasNoPart @@ -302,6 +304,7 @@ end function IO_open_JobFile_stat #ifndef Spectral +#ifndef FEM !-------------------------------------------------------------------------------------------------- !> @brief opens FEM input file for reading located in current working directory to given unit !-------------------------------------------------------------------------------------------------- @@ -367,6 +370,7 @@ subroutine IO_open_logFile(fileUnit) end subroutine IO_open_logFile #endif +#endif !-------------------------------------------------------------------------------------------------- diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 209baae0f..f1b868fff 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -29,6 +29,9 @@ module constitutive constitutive_aTolState !< pointer array to absolute state tolerance type(p_vec), public, dimension(:,:,:,:), allocatable :: & constitutive_RKCK45dotState !< pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method + real(pReal), public, dimension(:,:,:), allocatable :: & + constitutive_localDamage, & + constitutive_gradientDamage integer(pInt), public, dimension(:,:,:), allocatable :: & constitutive_sizeDotState, & !< size of dotState array constitutive_sizeState, & !< size of state array per grain @@ -240,6 +243,8 @@ subroutine constitutive_init #ifndef NEWSTATE ! lumped into new state allocate(constitutive_state0(cMax,iMax,eMax)) + allocate(constitutive_localDamage(cMax,iMax,eMax)); constitutive_localDamage = 1.0_pReal + allocate(constitutive_gradientDamage(cMax,iMax,eMax)); constitutive_gradientDamage = 1.0_pReal allocate(constitutive_partionedState0(cMax,iMax,eMax)) allocate(constitutive_subState0(cMax,iMax,eMax)) allocate(constitutive_state(cMax,iMax,eMax)) @@ -775,7 +780,7 @@ use math, only : & C = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el)) FeT = math_transpose33(Fe) - T = 0.5_pReal * math_mul3333xx33(C, math_mul33x33(FeT,Fe)-MATH_I3) + T = 0.5_pReal*math_mul3333xx33(C,math_mul33x33(FeT,Fe)-MATH_I3)*constitutive_gradientDamage(ipc,ip,el) dT_dFe = 0.0_pReal forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt, k=1_pInt:3_pInt, l=1_pInt:3_pInt) & diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 084198935..f956f52fe 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -551,7 +551,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) #else mappingConstitutive, & #endif - constitutive_TandItsTangent + constitutive_TandItsTangent, & + constitutive_localDamage, & + constitutive_gradientDamage implicit none logical, intent(in) :: & @@ -1371,6 +1373,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) crystallite_heat(g,i,e) = 0.98_pReal* & abs(math_mul33xx33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), & crystallite_Lp(1:3,1:3,g,i,e))) + constitutive_localDamage(g,i,e) = & + 1.0_pReal* & + sum(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)/constitutive_gradientDamage(g,i,e))* & + (math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)), & + crystallite_Fe(1:3,1:3,g,i,e))-math_I3))/4.0_pReal + & + 0.0_pReal* & + sum(abs(math_mul33x33(math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)), & + crystallite_Fp(1:3,1:3,g,i,e))-math_I3)/2.0_pReal) enddo enddo enddo elementLooping12 diff --git a/code/libs.f90 b/code/libs.f90 index 009b20805..740cbf2c9 100644 --- a/code/libs.f90 +++ b/code/libs.f90 @@ -11,6 +11,9 @@ end module libs #ifdef Spectral #include "../lib/kdtree2.f90" #endif +#ifdef FEM +#include "../lib/kdtree2.f90" +#endif #include "../lib/IR_Precision.f90" #include "../lib/Lib_Base64.f90" #include "../lib/Lib_VTK_IO.f90" diff --git a/code/numerics.f90 b/code/numerics.f90 index c92c9f221..0ce69aac4 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -12,6 +12,9 @@ module numerics implicit none private +#ifdef FEM +#include +#endif character(len=64), parameter, private :: & numerics_CONFIGFILE = 'numerics.config' !< name of configuration file @@ -89,6 +92,7 @@ module numerics fftw_planner_flag = 32_pInt, & !< conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw itmax = 250_pInt, & !< maximum number of iterations itmin = 2_pInt, & !< minimum number of iterations + stagItMax = 25_pInt, & !< maximum number of iterations maxCutBack = 3_pInt, & !< max number of cut backs continueCalculation = 0_pInt, & !< 0: exit if BVP solver does not converge, 1: continue calculation if BVP solver does not converge divergence_correction = 2_pInt !< correct divergence calculation in fourier space 0: no correction, 1: size scaled to 1, 2: size scaled to Npoints @@ -97,6 +101,33 @@ module numerics update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness #endif +!-------------------------------------------------------------------------------------------------- +! FEM parameters: +#ifdef FEM + real(pReal), protected, public :: & + err_struct_tolAbs = 1.0e-10_pReal, & !< absolute tolerance for equilibrium + err_struct_tolRel = 5.0e-4_pReal, & !< relative tolerance for equilibrium + err_thermal_tol = 1.0e-6_pReal, & + err_damage_tol = 1.0e-4_pReal + character(len=1024), protected, public :: & + petsc_optionsFEM = '-snes_type ngmres & + &-snes_ngmres_anderson ' + integer(pInt), protected, public :: & + itmaxFEM = 250_pInt, & !< maximum number of iterations + itminFEM = 2_pInt, & !< minimum number of iterations + maxCutBackFEM = 3_pInt, & !< max number of cut backs + integrationOrder = 1_pInt, & + structOrder = 2_pInt, & + thermalOrder = 2_pInt, & + damageOrder = 2_pInt, & + worldrank = 0_pInt, & + worldsize = 0_pInt + logical, protected, public :: & + structThermalCoupling = .false., & + structDamageCoupling = .false., & + thermalDamageCoupling = .false. +#endif + public :: numerics_init contains @@ -124,15 +155,20 @@ subroutine numerics_init #ifdef Spectral !$ use OMP_LIB, only: omp_set_num_threads ! Use the standard conforming module file for omp if using the spectral solver +#endif +#ifdef FEM +!$ use OMP_LIB, only: omp_set_num_threads ! Use the standard conforming module file for omp if using the spectral solver #endif implicit none #ifndef Spectral +#ifndef FEM !$ include "omp_lib.h" ! use the not F90 standard conforming include file to prevent crashes with some versions of MSC.Marc +#endif #endif integer(pInt), parameter :: FILEUNIT = 300_pInt ,& maxNchunks = 2_pInt !$ integer :: gotDAMASK_NUM_THREADS = 1 - integer :: i ! no pInt + integer :: i, ierr ! no pInt integer(pInt), dimension(1+2*maxNchunks) :: positions character(len=65536) :: & tag ,& @@ -154,6 +190,11 @@ subroutine numerics_init !$ endif !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution +#ifdef FEM + call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr) + call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,ierr);CHKERRQ(ierr) +#endif + !-------------------------------------------------------------------------------------------------- ! try to open the config file fileExists: if(IO_open_file_stat(FILEUNIT,numerics_configFile)) then @@ -274,6 +315,8 @@ subroutine numerics_init itmax = IO_intValue(line,positions,2_pInt) case ('itmin') itmin = IO_intValue(line,positions,2_pInt) + case ('stagitmax') + stagItMax = IO_intValue(line,positions,2_pInt) case ('maxcutback') maxCutBack = IO_intValue(line,positions,2_pInt) case ('continuecalculation') @@ -312,12 +355,53 @@ subroutine numerics_init #endif #ifndef Spectral case ('err_div_tolabs','err_div_tolrel','err_stress_tolrel','err_stress_tolabs',& ! found spectral parameter for FEM build - 'itmax', 'itmin','memory_efficient','fftw_timelimit','fftw_plan_mode', & + 'itmax', 'itmin','stagitmax','memory_efficient','fftw_timelimit','fftw_plan_mode', & 'divergence_correction','update_gamma','myfilter', & 'err_curl_tolabs','err_curl_tolrel', & 'maxcutback','polaralpha','polarbeta') call IO_warning(40_pInt,ext_msg=tag) #endif + +!-------------------------------------------------------------------------------------------------- +! FEM parameters +#ifdef FEM + case ('err_struct_tolabs') + err_struct_tolAbs = IO_floatValue(line,positions,2_pInt) + case ('err_struct_tolrel') + err_struct_tolRel = IO_floatValue(line,positions,2_pInt) + case ('err_thermal_tol') + err_thermal_tol = IO_floatValue(line,positions,2_pInt) + case ('err_damage_tol') + err_damage_tol = IO_floatValue(line,positions,2_pInt) + case ('itmaxfem') + itmaxFEM = IO_intValue(line,positions,2_pInt) + case ('itminfem') + itminFEM = IO_intValue(line,positions,2_pInt) + case ('maxcutbackfem') + maxCutBackFEM = IO_intValue(line,positions,2_pInt) + case ('integrationorder') + integrationorder = IO_intValue(line,positions,2_pInt) + case ('structorder') + structorder = IO_intValue(line,positions,2_pInt) + case ('thermalorder') + thermalorder = IO_intValue(line,positions,2_pInt) + case ('damageorder') + damageorder = IO_intValue(line,positions,2_pInt) + case ('petsc_optionsfem') + petsc_optionsFEM = trim(line(positions(4):)) + case ('structthermalcoupling') + structThermalCoupling = IO_intValue(line,positions,2_pInt) > 0_pInt + case ('structdamagecoupling') + structDamageCoupling = IO_intValue(line,positions,2_pInt) > 0_pInt + case ('thermaldamagecoupling') + thermalDamageCoupling = IO_intValue(line,positions,2_pInt) > 0_pInt +#endif +#ifndef FEM + case ('err_struct_tolabs','err_struct_tolrel','err_thermal_tol','err_damage_tol',& ! found spectral parameter for FEM build + 'itmaxfem', 'itminfem','maxcutbackfem','integrationorder','structorder','thermalorder', & + 'damageorder','petsc_optionsfem','structthermalcoupling','structdamagecoupling','thermaldamagecoupling') + call IO_warning(40_pInt,ext_msg=tag) +#endif case default ! found unknown keyword call IO_error(300_pInt,ext_msg=tag) endselect @@ -405,6 +489,7 @@ subroutine numerics_init #ifdef Spectral write(6,'(a24,1x,i8)') ' itmax: ',itmax write(6,'(a24,1x,i8)') ' itmin: ',itmin + write(6,'(a24,1x,i8)') ' stagItMax: ',stagItMax write(6,'(a24,1x,i8)') ' maxCutBack: ',maxCutBack write(6,'(a24,1x,i8)') ' continueCalculation: ',continueCalculation write(6,'(a24,1x,L8)') ' memory_efficient: ',memory_efficient @@ -432,6 +517,26 @@ subroutine numerics_init #endif #endif +!-------------------------------------------------------------------------------------------------- +! spectral parameters +#ifdef FEM + write(6,'(a24,1x,i8)') ' itmaxFEM: ',itmaxFEM + write(6,'(a24,1x,i8)') ' itminFEM: ',itminFEM + write(6,'(a24,1x,i8)') ' maxCutBackFEM: ',maxCutBackFEM + 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,es8.1)') ' err_struct_tolAbs: ',err_struct_tolAbs + write(6,'(a24,1x,es8.1)') ' err_struct_tolRel: ',err_struct_tolRel + write(6,'(a24,1x,es8.1)') ' err_thermal_tol: ',err_thermal_tol + write(6,'(a24,1x,es8.1)') ' err_damage_tol: ',err_damage_tol + write(6,'(a24,1x,L8)') ' structThermalCoupling: ',structThermalCoupling + write(6,'(a24,1x,L8)') ' structDamageCoupling: ',structDamageCoupling + write(6,'(a24,1x,L8)') ' thermalDamageCoupling: ',thermalDamageCoupling + write(6,'(a24,1x,a)') ' PETSc_optionsFEM: ',trim(petsc_optionsFEM) +#endif + !-------------------------------------------------------------------------------------------------- ! sanity checks if (relevantStrain <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relevantStrain') @@ -473,7 +578,7 @@ subroutine numerics_init if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='volDiscrPw_RGC') #ifdef Spectral if (itmax <= 1_pInt) call IO_error(301_pInt,ext_msg='itmax') - if (itmin > itmax .or. itmin < 1_pInt) call IO_error(301_pInt,ext_msg='itmin') + if (itmin > itmax .or. itmin < 0_pInt) call IO_error(301_pInt,ext_msg='itmin') if (continueCalculation /= 0_pInt .and. & continueCalculation /= 1_pInt) call IO_error(301_pInt,ext_msg='continueCalculation') if (divergence_correction < 0_pInt .or. & @@ -494,6 +599,16 @@ subroutine numerics_init polarBeta > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarBeta') #endif #endif +#ifdef FEM + if (itmaxFEM <= 1_pInt) call IO_error(301_pInt,ext_msg='itmaxFEM') + if (itminFEM > itmaxFEM .or. & + itminFEM < 0_pInt) call IO_error(301_pInt,ext_msg='itminFEM') + if (maxCutBackFEM < 0_pInt) call IO_error(301_pInt,ext_msg='maxCutBackFEM') + if (err_struct_tolRel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_struct_tolRel') + if (err_struct_tolAbs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_struct_tolAbs') + if (err_thermal_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_thermal_tol') + if (err_damage_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_damage_tol') +#endif end subroutine numerics_init diff --git a/code/prec.f90 b/code/prec.f90 index 338820a3c..dca64e2da 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -18,6 +18,9 @@ module prec #if (FLOAT==4) #ifdef Spectral SPECTRAL SOLVER DOES NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION +#endif +#ifdef FEM + SPECTRAL SOLVER DOES NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION #endif integer, parameter, public :: pReal = 4 !< floating point single precition (was selected_real_kind(6,37), number with 6 significant digits, up to 1e+-37) #ifdef __INTEL_COMPILER