diff --git a/Makefile_bk b/Makefile_bk new file mode 100755 index 000000000..8be738090 --- /dev/null +++ b/Makefile_bk @@ -0,0 +1,42 @@ +SHELL = /bin/sh +######################################################################################## +# Makefile for the installation of DAMASK +######################################################################################## +.PHONY: all +all: spectral marc processing + +.PHONY: spectral +spectral: + $(MAKE) DAMASK_spectral.exe -C code + +.PHONY: FEM +FEM: + $(MAKE) DAMASK_FEM.exe -C code + +.PHONY: marc +marc: + @./installation/mods_MarcMentat/apply_DAMASK_modifications.sh ${MAKEFLAGS} + +.PHONY: processing +processing: + @if hash cython 2>/dev/null; then \ + cd ./lib/damask; \ + CC=gcc python setup_corientation.py build_ext --inplace; \ + rm -rv build; \ + rm *.c; \ + fi + @./installation/compile_CoreModule.py ${MAKEFLAGS} + +.PHONY: tidy +tidy: + @$(MAKE) tidy -C code >/dev/null + +.PHONY: clean +clean: + @$(MAKE) cleanDAMASK -C code >/dev/null + +.PHONY: install +install: + @./installation/symlink_Code.py ${MAKEFLAGS} + @./installation/symlink_Processing.py ${MAKEFLAGS} + diff --git a/code/DAMASK_marc2011.f90 b/code/DAMASK_marc2011.f90 new file mode 120000 index 000000000..2c5bec706 --- /dev/null +++ b/code/DAMASK_marc2011.f90 @@ -0,0 +1 @@ +DAMASK_marc.f90 \ No newline at end of file diff --git a/code/DAMASK_marc2012.f90 b/code/DAMASK_marc2012.f90 new file mode 120000 index 000000000..2c5bec706 --- /dev/null +++ b/code/DAMASK_marc2012.f90 @@ -0,0 +1 @@ +DAMASK_marc.f90 \ No newline at end of file diff --git a/code/DAMASK_marc2013.1.f90 b/code/DAMASK_marc2013.1.f90 new file mode 120000 index 000000000..2c5bec706 --- /dev/null +++ b/code/DAMASK_marc2013.1.f90 @@ -0,0 +1 @@ +DAMASK_marc.f90 \ No newline at end of file diff --git a/code/DAMASK_marc2013.f90 b/code/DAMASK_marc2013.f90 new file mode 120000 index 000000000..2c5bec706 --- /dev/null +++ b/code/DAMASK_marc2013.f90 @@ -0,0 +1 @@ +DAMASK_marc.f90 \ No newline at end of file diff --git a/code/DAMASK_marc2014.2.f90 b/code/DAMASK_marc2014.2.f90 new file mode 120000 index 000000000..2c5bec706 --- /dev/null +++ b/code/DAMASK_marc2014.2.f90 @@ -0,0 +1 @@ +DAMASK_marc.f90 \ No newline at end of file diff --git a/code/DAMASK_marc2014.f90 b/code/DAMASK_marc2014.f90 new file mode 120000 index 000000000..2c5bec706 --- /dev/null +++ b/code/DAMASK_marc2014.f90 @@ -0,0 +1 @@ +DAMASK_marc.f90 \ No newline at end of file diff --git a/code/DAMASK_marc2015.f90 b/code/DAMASK_marc2015.f90 new file mode 120000 index 000000000..2c5bec706 --- /dev/null +++ b/code/DAMASK_marc2015.f90 @@ -0,0 +1 @@ +DAMASK_marc.f90 \ No newline at end of file diff --git a/code/Makefile b/code/Makefile index c539cc31e..a18cbebac 100644 --- a/code/Makefile +++ b/code/Makefile @@ -30,13 +30,13 @@ LINKERNAME ?= $(FLINKER) # # setting up for HDF5 support (hard link for now) # 1. Location of HDF5 binaries (with include/ and lib/ underneath) -HDF5 = /mnt/research/CMM/opt/hdf5-1.8.16 +HDF5 = /mnt/research/CMM/opt/hdf5 # 2. Location of External Libraries (missing in the 1.8.12 version) -LIBZ = $(HDF5)/lib/libz.a -LIBSZ = $(HDF5)/lib/libsz.a +LIBZ = /mnt/research/CMM/opt/hdf5/lib/libz.a +LIBSZ = /mnt/research/CMM/opt/hdf5/lib/libsz.a # 3. Set libraries for HDF5 (LIBS: shared lib, LIBZ: external lib) -HDFLIBS = -I$(HDF5)/include $(HDF5)/lib/libhdf5_fortran.a $(HDF5)/lib/libhdf5.a -HDFLIBZ = $(LIBZ) $(LIBSZ) -lm +HDFLIBS = -I$(HDF5)/include -L$(HDF5)/lib +HDFLIBZ = -L$(LIBZ) -L$(LIBSZ) # MPI compiler wrappers will tell if they are pointing to ifort or gfortran COMPILEROUT :=$(shell $(FC) -show) @@ -306,8 +306,8 @@ PRECISION_gfortran :=-fdefault-real-8 -fdefault-double-8 -DFLOAT=8 -DINT=4 #-fdefault-integer-8: Use it to set precision to 8 bytes for integer, don't use it for the standard case of pInt=4 (there is no -fdefault-integer-4) ################################################################################################### -COMPILE =$(OPENMP_FLAG_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) -COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) +COMPILE =$(OPENMP_FLAG_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) +COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) ################################################################################################### SOURCE_FILES = \ source_thermal_dissipation.o source_thermal_externalheat.o \ @@ -342,9 +342,6 @@ HYDROGENFLUX_FILES = \ HOMOGENIZATION_FILES = \ homogenization_RGC.o homogenization_isostrain.o homogenization_none.o -HDF5_FILES = \ - damask_hdf5.o - ##################### # Spectral Solver ##################### @@ -364,7 +361,7 @@ DAMASK_spectral.o: INTERFACENAME := spectral_interface.f90 SPECTRAL_SOLVER_FILES = spectral_mech_AL.o spectral_mech_Basic.o spectral_mech_Polarisation.o \ spectral_thermal.o spectral_damage.o -SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o \ +SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o damask_hdf5.o \ FEsolving.o mesh.o material.o lattice.o \ $(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \ crystallite.o \ @@ -372,10 +369,10 @@ SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o $(HOMOGENIZATION_FILES) homogenization.o \ CPFEM2.o \ spectral_utilities.o \ - $(HDF5_FILES) \ $(SPECTRAL_SOLVER_FILES) -DAMASK_spectral.exe: DAMASK_spectral.o +DAMASK_spectral.exe: DAMASK_spectral.o \ + $(SPECTRAL_FILES) $(PREFIX) $(LINKERNAME) $(OPENMP_FLAG_$(F90)) $(LINK_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) \ -o DAMASK_spectral.exe DAMASK_spectral.o \ $(SPECTRAL_FILES) $(LIBRARIES) $(HDFLIBS) $(HDFLIBZ) $(SUFFIX) @@ -662,10 +659,12 @@ DAMASK_interface.o: spectral_interface.f90 \ # -diag-disable 7410 should disable warning about directory statement in inquire function, but does not work. hence the other 2 statements prec.o: prec.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 $(SUFFIX) -endif -damask_hdf5.o: damask_hdf5.f90 - $(PREFIX) $(COMPILERNAME) $(COMPILE) -c damask_hdf5.f90 $(HDFLIBS) $(HDFLIBZ) +damask_hdf5.o: damask_hdf5.f90 \ + prec.o \ + IO.o + $(PREFIX) $(COMPILERNAME) $(HDFLIBS) $(HDFLIBZ) -c damask_hdf5.f90 $(SUFFIX) -lm +endif %.o : %.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) -c $< $(SUFFIX) diff --git a/code/Makefile_bk b/code/Makefile_bk new file mode 100644 index 000000000..a18cbebac --- /dev/null +++ b/code/Makefile_bk @@ -0,0 +1,699 @@ +SHELL = /bin/sh +######################################################################################## +# Makefile to compile the Material subroutine for BVP solution using spectral method +######################################################################################## +# Be sure to remove all files compiled with different options by using "make clean" +######################################################################################## +# OPTIONS = standard (alternative): meaning +#------------------------------------------------------------- +# F90 = ifort (gfortran): compiler type, choose Intel or GNU +# COMPILERNAME = name of the compiler executable (if not the same as the ype), e.g. using mpich-g90 instead of ifort +# PORTABLE = TRUE (FALSE): decision, if executable is optimized for the machine on which it was built. +# OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE,ULTRA): Optimization mode: O2, O0, O3 + further options for most files, O3 + further options for all files +# OPENMP = TRUE (FALSE): OpenMP multiprocessor support +# PREFIX = arbitrary prefix (before compilername) +# OPTION = arbitrary option (just before file to compile) +# SUFFIX = arbitrary suffix (after file to compile) +# STANDARD_CHECK = checking for Fortran 2008, compiler dependend +######################################################################################## +# including PETSc files. PETSC_ARCH is loaded from these files. +DAMASKVERSION :=$(shell cat ../VERSION) + +include ${PETSC_DIR}/lib/petsc/conf/variables +include ${PETSC_DIR}/lib/petsc/conf/rules + +INCLUDE_DIRS := $(PETSC_FC_INCLUDES) -DPETSc -I../lib +LIBRARIES := $(PETSC_WITH_EXTERNAL_LIB) +COMPILERNAME ?= $(FC) +LINKERNAME ?= $(FLINKER) + +# +# setting up for HDF5 support (hard link for now) +# 1. Location of HDF5 binaries (with include/ and lib/ underneath) +HDF5 = /mnt/research/CMM/opt/hdf5 +# 2. Location of External Libraries (missing in the 1.8.12 version) +LIBZ = /mnt/research/CMM/opt/hdf5/lib/libz.a +LIBSZ = /mnt/research/CMM/opt/hdf5/lib/libsz.a +# 3. Set libraries for HDF5 (LIBS: shared lib, LIBZ: external lib) +HDFLIBS = -I$(HDF5)/include -L$(HDF5)/lib +HDFLIBZ = -L$(LIBZ) -L$(LIBSZ) + +# MPI compiler wrappers will tell if they are pointing to ifort or gfortran +COMPILEROUT :=$(shell $(FC) -show) +# search in FC or COMPILEROUT for gfortran/ifort if not defined +ifeq ($(strip $(F90)),) + F90 :=$(findstring gfortran,$(FC) $(COMPILEROUT)) +endif +ifeq ($(strip $(F90)),) + F90 :=$(findstring ifort,$(FC) $(COMPILEROUT)) +endif + +OPENMP ?= ON +OPTIMIZATION ?= DEFENSIVE + +ifeq "$(OPTIMIZATION)" "OFF" +OPTI := OFF +MAXOPTI := OFF +endif +ifeq "$(OPTIMIZATION)" "DEFENSIVE" +OPTI := DEFENSIVE +MAXOPTI := DEFENSIVE +endif +ifeq "$(OPTIMIZATION)" "AGGRESSIVE" +OPTI := AGGRESSIVE +MAXOPTI := DEFENSIVE +endif +ifeq "$(OPTIMIZATION)" "ULTRA" +OPTI := AGGRESSIVE +MAXOPTI := AGGRESSIVE +endif + +ifndef OPTI +OPTI := DEFENSIVE +MAXOPTI := DEFENSIVE +endif + +# settings for shared memory multicore support +ifeq "$(OPENMP)" "ON" +OPENMP_FLAG_ifort =-openmp -openmp-report0 -parallel +OPENMP_FLAG_gfortran =-fopenmp +endif + +ifdef STANDARD_CHECK +STANDARD_CHECK_ifort =$(STANDARD_CHECK) +STANDARD_CHECK_gfortran =$(STANDARD_CHECK) +endif + +STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics +STANDARD_CHECK_gfortran ?=-std=f2008ts -pedantic-errors + +#-pedantic: more strict on standard, enables some warnings +# -pedantic-errors: like pedantic, but errors instead of warnings +OPTIMIZATION_OFF_ifort :=-O0 -no-ip +OPTIMIZATION_OFF_gfortran :=-O0 +OPTIMIZATION_DEFENSIVE_ifort :=-O2 +OPTIMIZATION_DEFENSIVE_gfortran :=-O2 +OPTIMIZATION_AGGRESSIVE_ifort :=-ipo -O3 -no-prec-div -fp-model fast=2 -xHost #-fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost +OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 -ffast-math -funroll-loops -ftree-vectorize + + +LINK_OPTIONS_ifort :=-shared-intel +COMPILE_OPTIONS_ifort :=-DDAMASKVERSION=\"${DAMASKVERSION}\"\ + -fpp\ + -ftz\ + -assume byterecl,fpe_summary\ + -diag-disable 5268\ + -warn declarations\ + -warn general\ + -warn usage\ + -warn interfaces\ + -warn ignore_loc\ + -warn alignments\ + -warn unused + +################################################################################################### +#COMPILE SWITCHES +#-shared-intel: Link against shared Intel libraries instead of static ones +#-fpp: preprocessor +#-ftz: flush unterflow to zero, automatically set if O<0,1,2,3> >0 +#-assume byterecl record length is given in bytes (also set by -standard-semantics) +# fpe_summary print list of floating point exceptions occured during execution +#-fimplicit-none: assume "implicit-none" even if not present in source +#-diag-disable: disables warnings, where +# warning ID 5268: the text exceeds right hand column allowed on the line (we have only comments there) +#-warn: enables warnings, where +# declarations: any undeclared names (alternative name: -implicitnone) +# general: warning messages and informational messages are issued by the compiler +# usage: questionable programming practices +# interfaces: checks the interfaces of all SUBROUTINEs called and FUNCTIONs invoked in your compilation against an external set of interface blocks +# ignore_loc: %LOC is stripped from an actual argument +# alignments: data that is not naturally aligned +# unused: declared variables that are never used +# stderrors: warnings about Fortran standard violations are changed to errors (STANDARD_CHECK) +# +################################################################################################### +#MORE OPTIONS FOR DEBUGGING DURING COMPILATION +#-warn: enables warnings, where +# truncated_source: Determines whether warnings occur when source exceeds the maximum column width in fixed-format files. (too many warnings because we have comments beyond character 132) +# uncalled: Determines whether warnings occur when a statement function is never called +# all: +# -name as_is: case sensitive Fortran! + +DEBUG_OPTIONS_ifort :=-g\ + -traceback\ + -gen-interfaces\ + -fp-stack-check\ + -fp-model strict\ + -check bounds,format,output_conversion,pointers,uninit\ + -ftrapuv\ + -fpe-all0\ + -warn errors\ + -warn stderrors\ + -debug-parameters all + +################################################################################################### +#COMPILE SWITCHES FOR RUNTIME DEBUGGING +#-g: Generate symbolic debugging information in the object file +#-traceback: Generate extra information in the object file to provide source file traceback information when a severe error occurs at run time. +#-gen-interfaces: Generate an interface block for each routine. http://software.intel.com/en-us/blogs/2012/01/05/doctor-fortran-gets-explicit-again/ +#-fp-stack-check: Generate extra code after every function call to ensure that the floating-point (FP) stack is in the expected state. +#-ftrapuv Trap uninitalized variables +#-check: checks at runtime, where +# bounds: check if an array index is too small (<1) or too large! +# format: Checking for the data type of an item being formatted for output. +# output_conversion: Checking for the fit of data items within a designated format descriptor field. +# pointers: Checking for certain disassociated or uninitialized pointers or unallocated allocatable objects. +# uninit: Checking for uninitialized variables. +#-fpe-all0 capture all floating-point exceptions, sets -ftz automatically +#-warn: enables warnings, where +# errors: warnings are changed to errors +# stderrors: warnings about Fortran standard violations are changed to errors +# information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/ +################################################################################################### +#MORE OPTIONS FOR RUNTIME DEBUGGING +#-heap-arrays: should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits +#-check: checks at runtime, where +# arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays (performance?) +# stack: +LINK_OPTIONS_gfortran :=-Wl,-undefined,dynamic_lookup +COMPILE_OPTIONS_gfortran :=-DDAMASKVERSION=\"${DAMASKVERSION}\"\ + -xf95-cpp-input\ + -ffree-line-length-132\ + -fimplicit-none\ + -fmodule-private\ + -Wall\ + -Wextra\ + -Wcharacter-truncation\ + -Wunderflow\ + -Wsuggest-attribute=pure\ + -Wsuggest-attribute=noreturn\ + -Wconversion-extra\ + -Wimplicit-procedure\ + -Wno-unused-parameter +#-ffpe-summary=all only for newer gfortran +################################################################################################### +#COMPILE SWITCHES +#-shared +#-Wl,-undefined,dynamic_lookup:ensure to link against dynamic libraries +#-xf95-cpp-input: preprocessor +#-ffree-line-length-132: restrict line length to the standard 132 characters +#-ffpe-summary: print summary of floating point exeptions (‘invalid’, ‘zero’, ‘overflow’, ‘underflow’, ‘inexact’ and ‘denormal’) +#-fimplicit-none: assume "implicit-none" even if not present in source +#-fmodule-private: assume "private" even if not present in source +#-Wcharacter-truncation: warn if character expressions (strings) are truncated +#-Wunderflow: produce a warning when numerical constant expressions are encountered, which yield an UNDERFLOW during compilation +#-Wsuggest-attribute=pure: +#-Wsuggest-attribute=noreturn: +#-Wconversion-extra +#-Wimplicit-procedure +#-Wall: sets the following Fortran options: +# -Waliasing: warn about possible aliasing of dummy arguments. Specifically, it warns if the same actual argument is associated with a dummy argument with "INTENT(IN)" and a dummy argument with "INTENT(OUT)" in a call with an explicit interface. +# -Wampersand: checks if a character expression is continued proberly by an ampersand at the end of the line and at the beginning of the new line +# -Warray-bounds: checks if array reference is out of bounds at compile time. use -fcheck-bounds to also check during runtime +# -Wconversion: warn about implicit conversions between different type +# -Wsurprising: warn when "suspicious" code constructs are encountered. While technically legal these usually indicate that an error has been made. +# -Wc-binding-type: +# -Wintrinsics-std: only standard intrisics are available, e.g. "call flush(6)" will cause an error +# -Wno-tabs: do not allow tabs in source +# -Wintrinsic-shadow: warn if a user-defined procedure or module procedure has the same name as an intrinsic +# -Wline-truncation: +# -Wtarget-lifetime: +# -Wreal-q-constant: warn about real-literal-constants with 'q' exponent-letter +# -Wunused: a number of unused-xxx warnings +# these are general (non -Fortran options) implied by -Wall +# -Waddress +# -Warray-bounds (only with -O2) +# -Wc++11-compat +# -Wchar-subscripts +# -Wcomment +# -Wformat +# -Wmaybe-uninitialized +# -Wnonnull +# -Wparentheses +# -Wpointer-sign +# -Wreorder +# -Wreturn-type +# -Wsequence-point +# -Wstrict-aliasing +# -Wstrict-overflow=1 +# -Wswitch +# -Wtrigraphs +# -Wuninitialized +# -Wunknown-pragmas +# -Wunused-function +# -Wunused-label +# -Wunused-value +# -Wunused-variable +# -Wvolatile-register-var +#-Wextra: sets the following Fortran options: +# -Wunuses-parameter: +# -Wcompare-reals: +# these are general (non -Fortran options) implied by -Wextra +# -Wclobbered +# -Wempty-body +# -Wignored-qualifiers +# -Wmissing-field-initializers +# -Woverride-init +# -Wsign-compare +# -Wtype-limits +# -Wuninitialized +# -Wunused-but-set-parameter (only with -Wunused or -Wall) +# -Wno-globals + +################################################################################################### +#MORE OPTIONS FOR DEBUGGING DURING COMPILATION +#-Warray-temporarieswarnings: because we have many temporary arrays (performance issue?): +#-Wimplicit-interface: no interfaces for lapack routines +#-Wunsafe-loop-optimizations: warn if the loop cannot be optimized due to nontrivial assumptions. +#-Wstrict-overflow: + +DEBUG_OPTIONS_gfortran :=-g\ + -fbacktrace\ + -fdump-core\ + -fcheck=all\ + -ffpe-trap=invalid,zero,overflow + +################################################################################################### +#COMPILE SWITCHES FOR RUNTIME DEBUGGING +#-ffpe-trap=invalid,\ stop execution if floating point exception is detected (NaN is silent) +# zero,\ +# overflow +#-fcheck=all: sets the following Fortran options: +#array-temps +#bounds +#do +#mem +#pointer +#recursion +################################################################################################### +#MORE OPTIONS FOR RUNTIME DEBUGGING +#-ffpe-trap=precision,\ +# denormal, \ +# underflow + +ifeq "$(DEBUG)" "ON" +COMPILE_OPTIONS_$(F90) +=$(DEBUG_OPTIONS_$(F90)) +LINK_OPTIONS_$(F90) +=$(DEBUG_OPTIONS_$(F90)) +endif +LINK_OPTIONS_$(F90) += $(OPTIMIZATION_$(MAXOPTI)_$(F90)) + +PRECISION_ifort :=-real-size 64 -integer-size 32 -DFLOAT=8 -DINT=4 +#-real-size 32: set precision to one of those 32/64/128 (= 4/8/16 bytes) for standard real (=8 for pReal) +#-integer-size 16: set precision to one of those 16/32/64 (= 2/4/8 bytes) for standard integer (=4 for pInt) +PRECISION_gfortran :=-fdefault-real-8 -fdefault-double-8 -DFLOAT=8 -DINT=4 +#-fdefault-real-8: set precision to 8 bytes for standard real (=8 for pReal). Will set size of double to 16 bytes as long as -fdefault-double-8 is not set +#-fdefault-double-8: set precision to 8 bytes for double real, would be 16 bytes because -fdefault-real-8 is used +#-fdefault-integer-8: Use it to set precision to 8 bytes for integer, don't use it for the standard case of pInt=4 (there is no -fdefault-integer-4) + +################################################################################################### +COMPILE =$(OPENMP_FLAG_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) +COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) +################################################################################################### +SOURCE_FILES = \ + source_thermal_dissipation.o source_thermal_externalheat.o \ + source_damage_isoBrittle.o source_damage_isoDuctile.o source_damage_anisoBrittle.o source_damage_anisoDuctile.o \ + source_vacancy_phenoplasticity.o source_vacancy_irradiation.o source_vacancy_thermalfluc.o + +KINEMATICS_FILES = \ + kinematics_cleavage_opening.o kinematics_slipplane_opening.o \ + kinematics_thermal_expansion.o \ + kinematics_vacancy_strain.o kinematics_hydrogen_strain.o + +PLASTIC_FILES = \ + plastic_dislotwin.o plastic_disloUCLA.o plastic_isotropic.o plastic_j2.o \ + plastic_phenopowerlaw.o plastic_titanmod.o plastic_nonlocal.o plastic_none.o \ + plastic_phenoplus.o + +THERMAL_FILES = \ + thermal_isothermal.o thermal_adiabatic.o thermal_conduction.o + +DAMAGE_FILES = \ + damage_none.o damage_local.o damage_nonlocal.o + +VACANCYFLUX_FILES = \ + vacancyflux_isoconc.o vacancyflux_isochempot.o vacancyflux_cahnhilliard.o + +POROSITY_FILES = \ + porosity_none.o porosity_phasefield.o + +HYDROGENFLUX_FILES = \ + hydrogenflux_isoconc.o hydrogenflux_cahnhilliard.o + +HOMOGENIZATION_FILES = \ + homogenization_RGC.o homogenization_isostrain.o homogenization_none.o + +##################### +# Spectral Solver +##################### +DAMASK_spectral.exe: IGNORE := \# +DAMASK_spectral.exe: COMPILE += -DSpectral +DAMASK_spectral.exe: COMPILE_MAXOPTI += -DSpectral +DAMASK_spectral.exe: MESHNAME := mesh.f90 +DAMASK_spectral.exe: INTERFACENAME := spectral_interface.f90 + +DAMASK_spectral.o: IGNORE := \# +DAMASK_spectral.o: COMPILE += -DSpectral +DAMASK_spectral.o: COMPILE_MAXOPTI += -DSpectral +DAMASK_spectral.o: MESHNAME := mesh.f90 +DAMASK_spectral.o: INTERFACENAME := spectral_interface.f90 + + +SPECTRAL_SOLVER_FILES = spectral_mech_AL.o spectral_mech_Basic.o spectral_mech_Polarisation.o \ + spectral_thermal.o spectral_damage.o + +SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o damask_hdf5.o \ + FEsolving.o mesh.o material.o lattice.o \ + $(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \ + crystallite.o \ + $(THERMAL_FILES) $(DAMAGE_FILES) $(VACANCYFLUX_FILES) $(HYDROGENFLUX_FILES) $(POROSITY_FILES) \ + $(HOMOGENIZATION_FILES) homogenization.o \ + CPFEM2.o \ + spectral_utilities.o \ + $(SPECTRAL_SOLVER_FILES) + +DAMASK_spectral.exe: DAMASK_spectral.o \ + $(SPECTRAL_FILES) + $(PREFIX) $(LINKERNAME) $(OPENMP_FLAG_$(F90)) $(LINK_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) \ + -o DAMASK_spectral.exe DAMASK_spectral.o \ + $(SPECTRAL_FILES) $(LIBRARIES) $(HDFLIBS) $(HDFLIBZ) $(SUFFIX) + + +DAMASK_spectral.o: DAMASK_spectral.f90 \ + $(SPECTRAL_SOLVER_FILES) + $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral.f90 $(SUFFIX) + +spectral_mech_AL.o: spectral_mech_AL.f90 \ + spectral_utilities.o + +spectral_mech_Polarisation.o: spectral_mech_Polarisation.f90 \ + spectral_utilities.o + +spectral_mech_Basic.o: spectral_mech_Basic.f90 \ + spectral_utilities.o + +spectral_thermal.o: spectral_thermal.f90 \ + spectral_utilities.o + +spectral_damage.o: spectral_damage.f90 \ + spectral_utilities.o + +spectral_utilities.o: spectral_utilities.f90 \ + CPFEM2.o + +##################### +# FEM Solver +##################### +VPATH := ../private/FEM/code +DAMASK_FEM.exe: COMPILE += -DFEM +DAMASK_FEM.exe: COMPILE_MAXOPTI += -DFEM +DAMASK_FEM.exe: MESHNAME := ../private/FEM/code/meshFEM.f90 +DAMASK_FEM.exe: INTERFACENAME := ../private/FEM/code/DAMASK_FEM_interface.f90 +DAMASK_FEM.exe: INCLUDE_DIRS += -I./ + +FEM_SOLVER_FILES = FEM_mech.o FEM_thermal.o FEM_damage.o FEM_vacancyflux.o FEM_porosity.o FEM_hydrogenflux.o + +FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o libs.o numerics.o debug.o math.o \ + FEsolving.o mesh.o material.o lattice.o \ + $(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \ + crystallite.o \ + $(THERMAL_FILES) $(DAMAGE_FILES) $(VACANCYFLUX_FILES) $(HYDROGENFLUX_FILES) $(POROSITY_FILES) \ + $(HOMOGENIZATION_FILES) homogenization.o \ + CPFEM.o \ + FEM_utilities.o $(FEM_SOLVER_FILES) + +DAMASK_FEM.exe: DAMASK_FEM_driver.o + $(PREFIX) $(LINKERNAME) $(OPENMP_FLAG_$(F90)) $(LINK_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) \ + -o DAMASK_FEM.exe DAMASK_FEM_driver.o \ + $(FEM_FILES) $(LIBRARIES) $(HDFLIBS) $(HDFLIBZ) $(SUFFIX) + +DAMASK_FEM_driver.o: DAMASK_FEM_driver.f90 $(FEM_SOLVER_FILES) + $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c ../private/FEM/code/DAMASK_FEM_driver.f90 $(SUFFIX) + +FEM_mech.o: FEM_mech.f90 \ + FEM_utilities.o + +FEM_thermal.o: FEM_thermal.f90 \ + FEM_utilities.o + +FEM_damage.o: FEM_damage.f90 \ + FEM_utilities.o + +FEM_vacancyflux.o: FEM_vacancyflux.f90 \ + FEM_utilities.o + +FEM_porosity.o: FEM_porosity.f90 \ + FEM_utilities.o + +FEM_hydrogenflux.o: FEM_hydrogenflux.f90 \ + FEM_utilities.o + +FEM_utilities.o: FEM_utilities.f90 \ + CPFEM.o + +FEZoo.o: $(wildcard FEZoo.f90) \ + IO.o + $(IGNORE) $(PREFIX) $(COMPILERNAME) $(COMPILE) -c ../private/FEM/code/FEZoo.f90 $(SUFFIX) + touch FEZoo.o + +CPFEM.o: CPFEM.f90 \ + homogenization.o + +CPFEM2.o: CPFEM2.f90 \ + homogenization.o + +homogenization.o: homogenization.f90 \ + $(THERMAL_FILES) \ + $(DAMAGE_FILES) \ + $(VACANCYFLUX_FILES) \ + $(POROSITY_FILES) \ + $(HYDROGENFLUX_FILES) \ + $(HOMOGENIZATION_FILES) + +thermal_isothermal.o: thermal_isothermal.f90 \ + crystallite.o + +thermal_adiabatic.o: thermal_adiabatic.f90 \ + crystallite.o + +thermal_conduction.o: thermal_conduction.f90 \ + crystallite.o + +damage_none.o: damage_none.f90 \ + crystallite.o + +damage_local.o: damage_local.f90 \ + crystallite.o + +damage_nonlocal.o: damage_nonlocal.f90 \ + crystallite.o + +thermal_conduction.o: thermal_conduction.f90 \ + crystallite.o + +vacancyflux_isoconc.o: vacancyflux_isoconc.f90 \ + crystallite.o + +vacancyflux_isochempot.o: vacancyflux_isochempot.f90 \ + crystallite.o + +vacancyflux_cahnhilliard.o: vacancyflux_cahnhilliard.f90 \ + crystallite.o + +porosity_none.o: porosity_none.f90 \ + crystallite.o + +porosity_phasefield.o: porosity_phasefield.f90 \ + crystallite.o + +hydrogenflux_isoconc.o: hydrogenflux_isoconc.f90 \ + crystallite.o + +hydrogenflux_cahnhilliard.o: hydrogenflux_cahnhilliard.f90 \ + crystallite.o + +homogenization_RGC.o: homogenization_RGC.f90 \ + crystallite.o + +homogenization_isostrain.o: homogenization_isostrain.f90 \ + crystallite.o + +homogenization_none.o: homogenization_none.f90 \ + crystallite.o + +crystallite.o: crystallite.f90 \ + constitutive.o + +constitutive.o: constitutive.f90 \ + $(SOURCE_FILES) \ + $(KINEMATICS_FILES) \ + $(PLASTIC_FILES) + +source_thermal_dissipation.o: source_thermal_dissipation.f90 \ + lattice.o + +source_thermal_externalheat.o: source_thermal_externalheat.f90 \ + lattice.o + +source_damage_isoBrittle.o: source_damage_isoBrittle.f90 \ + lattice.o + +source_damage_isoDuctile.o: source_damage_isoDuctile.f90 \ + lattice.o + +source_damage_anisoBrittle.o: source_damage_anisoBrittle.f90 \ + lattice.o + +source_damage_anisoDuctile.o: source_damage_anisoDuctile.f90 \ + lattice.o + +source_vacancy_phenoplasticity.o: source_vacancy_phenoplasticity.f90 \ + lattice.o + +source_vacancy_irradiation.o: source_vacancy_irradiation.f90 \ + lattice.o + +source_vacancy_thermalfluc.o: source_vacancy_thermalfluc.f90 \ + lattice.o + +kinematics_cleavage_opening.o: kinematics_cleavage_opening.f90 \ + lattice.o + +kinematics_slipplane_opening.o: kinematics_slipplane_opening.f90 \ + lattice.o + +kinematics_thermal_expansion.o: kinematics_thermal_expansion.f90 \ + lattice.o + +kinematics_vacancy_strain.o: kinematics_vacancy_strain.f90 \ + lattice.o + +kinematics_hydrogen_strain.o: kinematics_hydrogen_strain.f90 \ + lattice.o + +plastic_nonlocal.o: plastic_nonlocal.f90 \ + lattice.o + +plastic_titanmod.o: plastic_titanmod.f90 \ + lattice.o + +plastic_disloUCLA.o: plastic_disloUCLA.f90 \ + lattice.o + +plastic_dislotwin.o: plastic_dislotwin.f90 \ + lattice.o + +plastic_phenopowerlaw.o: plastic_phenopowerlaw.f90 \ + lattice.o + +plastic_phenoplus.o: plastic_phenoplus.f90 \ + lattice.o + +plastic_isotropic.o: plastic_isotropic.f90 \ + lattice.o + +plastic_j2.o: plastic_j2.f90 \ + lattice.o + +plastic_none.o: plastic_none.f90 \ + lattice.o +ifeq "$(F90)" "gfortran" +lattice.o: lattice.f90 \ + material.o + $(PREFIX) $(COMPILERNAME) $(COMPILE) -ffree-line-length-240 -c lattice.f90 $(SUFFIX) +# long lines for interaction matrix +else +lattice.o: lattice.f90 \ + material.o +endif + +material.o: material.f90 \ + mesh.o + +mesh.o: mesh.f90 \ + $(wildcard meshFEM.f90) \ + FEsolving.o \ + math.o \ + FEZoo.o + $(PREFIX) $(COMPILERNAME) $(COMPILE) -c $(MESHNAME) -o mesh.o $(SUFFIX) + +FEsolving.o: FEsolving.f90 \ + debug.o + +math.o: math.f90 \ + debug.o + +debug.o: debug.f90 \ + numerics.o + +numerics.o: numerics.f90 \ + libs.o + +libs.o: libs.f90 \ + IO.o + +IO.o: IO.f90 \ + DAMASK_interface.o + +ifeq "$(F90)" "gfortran" +DAMASK_interface.o: spectral_interface.f90 \ + $(wildcard DAMASK_FEM_interface.f90) \ + prec.o + $(PREFIX) $(COMPILERNAME) $(COMPILE) -c $(INTERFACENAME) -fall-intrinsics -o DAMASK_interface.o $(SUFFIX) +#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored +# and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external +# --> allows the use of 'getcwd' +prec.o: prec.f90 + $(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 -fno-range-check -fall-intrinsics -fno-fast-math $(SUFFIX) +# fno-range-check: Disable range checking on results of simplification of constant expressions during compilation +# --> allows the definition of DAMASK_NaN +#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored +# and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external +# --> allows the use of 'isnan' +#-fno-fast-math: +# --> otherwise, when setting -ffast-math, isnan always evaluates to false (I would call it a bug) +else +DAMASK_interface.o: spectral_interface.f90 \ + $(wildcard DAMASK_FEM_interface.f90) \ + prec.o + $(PREFIX) $(COMPILERNAME) $(COMPILE) -c $(INTERFACENAME) -diag-remark 7410 -stand none -warn nostderrors -o DAMASK_interface.o $(SUFFIX) +# -diag-disable 7410 should disable warning about directory statement in inquire function, but does not work. hence the other 2 statements +prec.o: prec.f90 + $(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 $(SUFFIX) + +damask_hdf5.o: damask_hdf5.f90 \ + prec.o \ + IO.o + $(PREFIX) $(COMPILERNAME) $(HDFLIBS) $(HDFLIBZ) -c damask_hdf5.f90 $(SUFFIX) -lm +endif + +%.o : %.f90 + $(PREFIX) $(COMPILERNAME) $(COMPILE) -c $< $(SUFFIX) + +.PHONY: tidy +tidy: + @rm -rf *.o + @rm -rf *.mod + @rm -rf *.inst.f90 # for instrumentation + @rm -rf *.pomp.f90 # for instrumentation + @rm -rf *.pp.f90 # for instrumentation + @rm -rf *.pdb # for instrumnentation + @rm -rf *.opari.inc # for instrumnentation + +.PHONY: cleanDAMASK +cleanDAMASK: + @rm -rf *.exe + @rm -rf *.marc + @rm -rf *.o + @rm -rf *.mod + @rm -rf *.inst.f90 # for instrumentation + @rm -rf *.pomp.f90 # for instrumentation + @rm -rf *.pp.f90 # for instrumentation + @rm -rf *.pdb # for instrumentation + @rm -rf *.opari.inc # for instrumentation + +.PHONY: help +help: + F90="$(F90)" + COMPILERNAME="$(COMPILERNAME)" + COMPILEROUT="$(COMPILEROUT)" + diff --git a/code/quit__genmod.f90 b/code/quit__genmod.f90 new file mode 100644 index 000000000..c6b282470 --- /dev/null +++ b/code/quit__genmod.f90 @@ -0,0 +1,8 @@ + !COMPILER-GENERATED INTERFACE MODULE: Tue Feb 23 16:12:31 2016 + MODULE QUIT__genmod + INTERFACE + SUBROUTINE QUIT(STOP_ID) + INTEGER(KIND=4), INTENT(IN) :: STOP_ID + END SUBROUTINE QUIT + END INTERFACE + END MODULE QUIT__genmod diff --git a/code/spectral/CMakeLists.txt b/code/spectral/CMakeLists.txt new file mode 100644 index 000000000..a11ff1303 --- /dev/null +++ b/code/spectral/CMakeLists.txt @@ -0,0 +1,14 @@ +# group source for sepctral solver driver +set (SPECTRAL "spectral_damage" + "spectral_interface" + "spectral_mech_AL" + "spectral_mech_Basic" + "spectral_mech_Polarisation" + "spectral_thermal" + "spectral_utilities" + ) + +# compile spectral solver driver module +foreach (p ${SPECTRAL}) + add_library (${p} MODULE "${p}.f90") +endforeach (p) diff --git a/code/spectral/spectral_damage.f90 b/code/spectral/spectral_damage.f90 new file mode 100644 index 000000000..0b79d5e5d --- /dev/null +++ b/code/spectral/spectral_damage.f90 @@ -0,0 +1,414 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: spectral_damage.f90 4082 2015-04-11 20:28:07Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Spectral solver for nonlocal damage +!-------------------------------------------------------------------------------------------------- +module spectral_damage + use prec, only: & + pInt, & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + use numerics, only: & + worldrank, & + worldsize + + implicit none + private +#include + + character (len=*), parameter, public :: & + spectral_damage_label = 'spectraldamage' + +!-------------------------------------------------------------------------------------------------- +! derived types + type(tSolutionParams), private :: params + +!-------------------------------------------------------------------------------------------------- +! PETSc data + SNES, private :: damage_snes + Vec, private :: solution + PetscInt, private :: xstart, xend, ystart, yend, zstart, zend + real(pReal), private, dimension(:,:,:), allocatable :: & + damage_current, & !< field of current damage + damage_lastInc, & !< field of previous damage + damage_stagInc !< field of staggered damage + +!-------------------------------------------------------------------------------------------------- +! reference diffusion tensor, mobility etc. + integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment + real(pReal), dimension(3,3), private :: D_ref + real(pReal), private :: mobility_ref + character(len=1024), private :: incInfo + + public :: & + spectral_damage_init, & + spectral_damage_solution, & + spectral_damage_forward, & + spectral_damage_destroy + external :: & + VecDestroy, & + DMDestroy, & + DMDACreate3D, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + PETScFinalize, & + SNESDestroy, & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber, & + SNESSolve, & + SNESSetDM, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions, & + SNESCreate, & + MPI_Abort, & + MPI_Bcast, & + MPI_Allreduce + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!-------------------------------------------------------------------------------------------------- +subroutine spectral_damage_init() + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & + IO_intOut, & + IO_read_realFile, & + IO_timeStamp + use spectral_utilities, only: & + wgt + use mesh, only: & + grid, & + grid3 + use damage_nonlocal, only: & + damage_nonlocal_getDiffusion33, & + damage_nonlocal_getMobility + + implicit none + DM :: damage_grid + Vec :: uBound, lBound + PetscErrorCode :: ierr + PetscObject :: dummy + integer(pInt), dimension(:), allocatable :: localK + integer(pInt) :: proc + integer(pInt) :: i, j, k, cell + character(len=100) :: snes_type + + mainProcess: if (worldrank == 0_pInt) then + write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + +!-------------------------------------------------------------------------------------------------- +! initialize solver specific parts of PETSc + call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 + do proc = 1, worldsize + call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) + enddo + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & !< cut off stencil at boundary + DMDA_STENCIL_BOX, & !< Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & !< global grid + 1, 1, worldsize, & + 1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap) + grid(1),grid(2),localK, & !< local grid + damage_grid,ierr) !< handle, error + CHKERRQ(ierr) + call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da + call DMCreateGlobalVector(damage_grid,solution,ierr); CHKERRQ(ierr) !< global solution vector (grid x 1, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,spectral_damage_formResidual,dummy,ierr) !< residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional cli arguments + call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) + if (trim(snes_type) == 'vinewtonrsls' .or. & + trim(snes_type) == 'vinewtonssls') then + call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) + call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) + call VecSet(lBound,0.0,ierr); CHKERRQ(ierr) + call VecSet(uBound,1.0,ierr); CHKERRQ(ierr) + call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc. + call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) + call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) + endif + +!-------------------------------------------------------------------------------------------------- +! init fields + call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr) + CHKERRQ(ierr) + xend = xstart + xend - 1 + yend = ystart + yend - 1 + zend = zstart + zend - 1 + call VecSet(solution,1.0,ierr); CHKERRQ(ierr) + allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal) + allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) + allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! damage reference diffusion update + cell = 0_pInt + D_ref = 0.0_pReal + mobility_ref = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) + mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) + enddo; enddo; enddo + D_ref = D_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + mobility_ref = mobility_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + +end subroutine spectral_damage_init + +!-------------------------------------------------------------------------------------------------- +!> @brief solution for the spectral damage scheme with internal iterations +!-------------------------------------------------------------------------------------------------- +type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old,loadCaseTime) + use numerics, only: & + itmax, & + err_damage_tolAbs, & + err_damage_tolRel + use spectral_utilities, only: & + tBoundaryCondition, & + Utilities_maskedCompliance, & + Utilities_updateGamma + use mesh, only: & + grid, & + grid3 + use damage_nonlocal, only: & + damage_nonlocal_putNonLocalDamage + + implicit none + +!-------------------------------------------------------------------------------------------------- +! input data for solution + 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 + logical, intent(in) :: guess + integer(pInt) :: i, j, k, cell + PetscInt ::position + PetscReal :: minDamage, maxDamage, stagNorm, solnNorm + +!-------------------------------------------------------------------------------------------------- +! PETSc Data + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + spectral_damage_solution%converged =.false. + +!-------------------------------------------------------------------------------------------------- +! set module wide availabe data + params%timeinc = timeinc + params%timeincOld = timeinc_old + + call SNESSolve(damage_snes,PETSC_NULL_OBJECT,solution,ierr); CHKERRQ(ierr) + call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) + + if (reason < 1) then + spectral_damage_solution%converged = .false. + spectral_damage_solution%iterationsNeeded = itmax + else + spectral_damage_solution%converged = .true. + spectral_damage_solution%iterationsNeeded = totalIter + endif + stagNorm = maxval(abs(damage_current - damage_stagInc)) + solnNorm = maxval(abs(damage_current)) + call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) + call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) + damage_stagInc = damage_current + spectral_damage_solution%stagConverged = stagNorm < err_damage_tolAbs & + .or. stagNorm < err_damage_tolRel*solnNorm + +!-------------------------------------------------------------------------------------------------- +! updating damage state + cell = 0_pInt !< material point = 0 + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt !< material point increase + call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) + enddo; enddo; enddo + + call VecMin(solution,position,minDamage,ierr); CHKERRQ(ierr) + call VecMax(solution,position,maxDamage,ierr); CHKERRQ(ierr) + if (worldrank == 0) then + if (spectral_damage_solution%converged) & + write(6,'(/,a)') ' ... nonlocal damage converged .....................................' + write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& + minDamage, maxDamage, stagNorm + write(6,'(/,a)') ' ===========================================================================' + flush(6) + endif + +end function spectral_damage_solution + + +!-------------------------------------------------------------------------------------------------- +!> @brief forms the spectral damage residual vector +!-------------------------------------------------------------------------------------------------- +subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr) + use numerics, only: & + residualStiffness + use mesh, only: & + grid, & + grid3 + use math, only: & + math_mul33x3 + use spectral_utilities, only: & + scalarField_real, & + vectorField_real, & + utilities_FFTvectorForward, & + utilities_FFTvectorBackward, & + utilities_FFTscalarForward, & + utilities_FFTscalarBackward, & + utilities_fourierGreenConvolution, & + utilities_fourierScalarGradient, & + utilities_fourierVectorDivergence + use damage_nonlocal, only: & + damage_nonlocal_getSourceAndItsTangent,& + damage_nonlocal_getDiffusion33, & + damage_nonlocal_getMobility + + implicit none + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & + in + PetscScalar, dimension( & + XG_RANGE,YG_RANGE,ZG_RANGE) :: & + x_scal + PetscScalar, dimension( & + X_RANGE,Y_RANGE,Z_RANGE) :: & + f_scal + PetscObject :: dummy + PetscErrorCode :: ierr + integer(pInt) :: i, j, k, cell + real(pReal) :: phiDot, dPhiDot_dPhi, mobility + + damage_current = x_scal +!-------------------------------------------------------------------------------------------------- +! evaluate polarization field + scalarField_real = 0.0_pReal + scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current + call utilities_FFTscalarForward() + call utilities_fourierScalarGradient() !< calculate gradient of damage field + call utilities_FFTvectorBackward() + cell = 0_pInt + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, & + vectorField_real(1:3,i,j,k)) + enddo; enddo; enddo + call utilities_FFTvectorForward() + call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field + call utilities_FFTscalarBackward() + cell = 0_pInt + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell) + mobility = damage_nonlocal_getMobility(1,cell) + scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & + params%timeinc*phiDot + & + mobility*damage_lastInc(i,j,k) - & + mobility*damage_current(i,j,k) + & + mobility_ref*damage_current(i,j,k) + enddo; enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! convolution of damage field with green operator + call utilities_FFTscalarForward() + call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) + call utilities_FFTscalarBackward() + where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) & + scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc + where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) & + scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness + +!-------------------------------------------------------------------------------------------------- +! constructing residual + f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current + +end subroutine spectral_damage_formResidual + +!-------------------------------------------------------------------------------------------------- +!> @brief spectral damage forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime) + use mesh, only: & + grid, & + grid3 + use spectral_utilities, only: & + cutBack, & + wgt + use damage_nonlocal, only: & + damage_nonlocal_putNonLocalDamage, & + damage_nonlocal_getDiffusion33, & + damage_nonlocal_getMobility + + implicit none + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + logical, intent(in) :: guess + PetscErrorCode :: ierr + integer(pInt) :: i, j, k, cell + DM :: dm_local + PetscScalar, dimension(:,:,:), pointer :: x_scal + + if (cutBack) then + damage_current = damage_lastInc + damage_stagInc = damage_lastInc +!-------------------------------------------------------------------------------------------------- +! reverting damage field state + cell = 0_pInt + call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with + x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current + call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) + enddo; enddo; enddo + else +!-------------------------------------------------------------------------------------------------- +! update rate and forward last inc + damage_lastInc = damage_current + cell = 0_pInt + D_ref = 0.0_pReal + mobility_ref = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) + mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) + enddo; enddo; enddo + D_ref = D_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + mobility_ref = mobility_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + endif + + end subroutine spectral_damage_forward + +!-------------------------------------------------------------------------------------------------- +!> @brief destroy routine +!-------------------------------------------------------------------------------------------------- +subroutine spectral_damage_destroy() + + implicit none + PetscErrorCode :: ierr + + call VecDestroy(solution,ierr); CHKERRQ(ierr) + call SNESDestroy(damage_snes,ierr); CHKERRQ(ierr) + +end subroutine spectral_damage_destroy + +end module spectral_damage diff --git a/code/spectral/spectral_interface.f90 b/code/spectral/spectral_interface.f90 new file mode 100644 index 000000000..b24c5f747 --- /dev/null +++ b/code/spectral/spectral_interface.f90 @@ -0,0 +1,568 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Interfacing between the spectral solver and the material subroutines provided +!! by DAMASK +!> @details Interfacing between the spectral solver and the material subroutines provided +!> by DAMASK. Interpretating the command line arguments or, in case of called from f2py, +!> the arguments parsed to the init routine to get load case, geometry file, working +!> directory, etc. +!-------------------------------------------------------------------------------------------------- +module DAMASK_interface + use prec, only: & + pInt + + implicit none + private +#ifdef PETSc +#include +#endif + logical, public, protected :: appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding) + integer(pInt), public, protected :: spectralRestartInc = 1_pInt !< Increment at which calculation starts + character(len=1024), public, protected :: & + geometryFile = '', & !< parameter given for geometry file + loadCaseFile = '' !< parameter given for load case file + character(len=1024), private :: workingDirectory !< accessed by getSolverWorkingDirectoryName for compatibility reasons + + public :: & + getSolverWorkingDirectoryName, & + getSolverJobName, & + DAMASK_interface_init + private :: & + storeWorkingDirectory, & + getGeometryFile, & + getLoadCaseFile, & + rectifyPath, & + makeRelativePath, & + getPathSep, & + IIO_stringValue, & + IIO_intValue, & + IIO_lc, & + IIO_stringPos + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief initializes the solver by interpreting the command line arguments. Also writes +!! information on computation to screen +!-------------------------------------------------------------------------------------------------- +subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + + implicit none + character(len=1024), optional, intent(in) :: & + loadCaseParameterIn, & !< if using the f2py variant, the -l argument of DAMASK_spectral.exe + geometryParameterIn !< if using the f2py variant, the -g argument of DAMASK_spectral.exe + character(len=1024) :: & + commandLine, & !< command line call as string + loadCaseArg ='', & !< -l argument given to DAMASK_spectral.exe + geometryArg ='', & !< -g argument given to DAMASK_spectral.exe + workingDirArg ='', & !< -w argument given to DAMASK_spectral.exe + hostName, & !< name of machine on which DAMASK_spectral.exe is execute (might require export HOSTNAME) + userName, & !< name of user calling DAMASK_spectral.exe + tag + integer :: & + i, & + worldrank = 0 + integer, allocatable, dimension(:) :: & + chunkPos + integer, dimension(8) :: & + dateAndTime ! type default integer +#ifdef PETSc + PetscErrorCode :: ierr +#endif + external :: & + quit,& + MPI_Comm_rank,& + PETScInitialize, & + MPI_abort + +!-------------------------------------------------------------------------------------------------- +! PETSc Init +#ifdef PETSc + call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code + CHKERRQ(ierr) ! this is a macro definition, it is case sensitive + + open(6, encoding='UTF-8') ! modern fortran compilers (gfortran >4.4, ifort >11 support it) + call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr) +#endif + mainProcess: if (worldrank == 0) then + call date_and_time(values = dateAndTime) + write(6,'(/,a)') ' <<<+- DAMASK_spectral -+>>>' + write(6,'(/,a)') ' Version: '//DAMASKVERSION + write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',& + dateAndTime(2),'/',& + dateAndTime(1) + write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',& + dateAndTime(6),':',& + dateAndTime(7) + write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>' +#include "compilation_info.f90" + endif mainProcess + + if ( present(loadcaseParameterIn) .and. present(geometryParameterIn)) then ! both mandatory parameters given in function call + geometryArg = geometryParameterIn + loadcaseArg = loadcaseParameterIn + commandLine = 'n/a' + else if ( .not.( present(loadcaseParameterIn) .and. present(geometryParameterIn))) then ! none parameters given in function call, trying to get them from command line + call get_command(commandLine) + chunkPos = IIO_stringPos(commandLine) + do i = 1, chunkPos(1) + tag = IIO_lc(IIO_stringValue(commandLine,chunkPos,i)) ! extract key + select case(tag) + case ('-h','--help') + mainProcess2: if (worldrank == 0) then + write(6,'(a)') ' #######################################################################' + write(6,'(a)') ' DAMASK_spectral:' + write(6,'(a)') ' The spectral method boundary value problem solver for' + write(6,'(a)') ' the Düsseldorf Advanced Material Simulation Kit' + write(6,'(a,/)')' #######################################################################' + write(6,'(a,/)')' Valid command line switches:' + write(6,'(a)') ' --geom (-g, --geometry)' + write(6,'(a)') ' --load (-l, --loadcase)' + write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)' + write(6,'(a)') ' --restart (-r, --rs)' + write(6,'(a)') ' --regrid (--rg)' + write(6,'(a)') ' --help (-h)' + write(6,'(/,a)')' -----------------------------------------------------------------------' + write(6,'(a)') ' Mandatory arguments:' + write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom.geom' + write(6,'(a)') ' Specifies the location of the geometry definition file,' + write(6,'(a)') ' if no extension is given, .geom will be appended.' + write(6,'(a)') ' "PathToGeomFile" will be the working directory if not specified' + write(6,'(a)') ' via --workingdir.' + write(6,'(a)') ' Make sure the file "material.config" exists in the working' + write(6,'(a)') ' directory.' + write(6,'(a)') ' For further configuration place "numerics.config"' + write(6,'(a)')' and "numerics.config" in that directory.' + write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile.load' + write(6,'(a)') ' Specifies the location of the load case definition file,' + write(6,'(a)') ' if no extension is given, .load will be appended.' + write(6,'(/,a)')' -----------------------------------------------------------------------' + write(6,'(a)') ' Optional arguments:' + write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory' + write(6,'(a)') ' Specifies the working directory and overwrites the default' + write(6,'(a)') ' "PathToGeomFile".' + write(6,'(a)') ' Make sure the file "material.config" exists in the working' + write(6,'(a)') ' directory.' + write(6,'(a)') ' For further configuration place "numerics.config"' + write(6,'(a)')' and "numerics.config" in that directory.' + write(6,'(/,a)')' --restart XX' + write(6,'(a)') ' Reads in total increment No. XX-1 and continues to' + write(6,'(a)') ' calculate total increment No. XX.' + write(6,'(a)') ' Appends to existing results file ' + write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' + write(6,'(a)') ' Works only if the restart information for total increment' + write(6,'(a)') ' No. XX-1 is available in the working directory.' + write(6,'(/,a)')' --regrid XX' + write(6,'(a)') ' Reads in total increment No. XX-1 and continues to' + write(6,'(a)') ' calculate total increment No. XX.' + write(6,'(a)') ' Attention: Overwrites existing results file ' + write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' + write(6,'(a)') ' Works only if the restart information for total increment' + write(6,'(a)') ' No. XX-1 is available in the working directory.' + write(6,'(/,a)')' -----------------------------------------------------------------------' + write(6,'(a)') ' Help:' + write(6,'(/,a)')' --help' + write(6,'(a,/)')' Prints this message and exits' + call quit(0_pInt) ! normal Termination + endif mainProcess2 + case ('-l', '--load', '--loadcase') + loadcaseArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt) + case ('-g', '--geom', '--geometry') + geometryArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt) + case ('-w', '-d', '--wd', '--directory', '--workingdir', '--workingdirectory') + workingDirArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt) + case ('-r', '--rs', '--restart') + spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt) + appendToOutFile = .true. + case ('--rg', '--regrid') + spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt) + appendToOutFile = .false. + end select + enddo + endif + + if (len(trim(loadcaseArg)) == 0 .or. len(trim(geometryArg)) == 0) then + write(6,'(a)') ' Please specify geometry AND load case (-h for help)' + call quit(1_pInt) + endif + + workingDirectory = storeWorkingDirectory(trim(workingDirArg),trim(geometryArg)) + geometryFile = getGeometryFile(geometryArg) + loadCaseFile = getLoadCaseFile(loadCaseArg) + + call get_environment_variable('HOSTNAME',hostName) + call get_environment_variable('USER',userName) + mainProcess3: if (worldrank == 0) then + write(6,'(a,a)') ' Host name: ', trim(hostName) + write(6,'(a,a)') ' User name: ', trim(userName) + write(6,'(a,a)') ' Path separator: ', getPathSep() + write(6,'(a,a)') ' Command line call: ', trim(commandLine) + if (len(trim(workingDirArg))>0) & + write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) + write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) + write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg) + write(6,'(a,a)') ' Working directory: ', trim(getSolverWorkingDirectoryName()) + write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) + write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) + write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) + if (SpectralRestartInc > 1_pInt) & + write(6,'(a,i6.6)') ' Restart at increment: ', spectralRestartInc + write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile + endif mainProcess3 + +end subroutine DAMASK_interface_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief extract working directory from given argument or from location of geometry file, +!! possibly converting relative arguments to absolut path +!> @todo change working directory with call chdir(storeWorkingDirectory)? +!-------------------------------------------------------------------------------------------------- +character(len=1024) function storeWorkingDirectory(workingDirectoryArg,geometryArg) +#ifdef __INTEL_COMPILER + use IFPORT +#endif + + implicit none + character(len=*), intent(in) :: workingDirectoryArg !< working directory argument + character(len=*), intent(in) :: geometryArg !< geometry argument + character(len=1024) :: cwd + character :: pathSep + logical :: dirExists + external :: quit + integer :: error + + pathSep = getPathSep() + if (len(workingDirectoryArg)>0) then ! got working directory as input + if (workingDirectoryArg(1:1) == pathSep) then ! absolute path given as command line argument + storeWorkingDirectory = workingDirectoryArg + else + error = getcwd(cwd) ! relative path given as command line argument + storeWorkingDirectory = trim(cwd)//pathSep//workingDirectoryArg + endif + if (storeWorkingDirectory(len(trim(storeWorkingDirectory)):len(trim(storeWorkingDirectory))) & ! if path seperator is not given, append it + /= pathSep) storeWorkingDirectory = trim(storeWorkingDirectory)//pathSep +#ifdef __INTEL_COMPILER + inquire(directory = trim(storeWorkingDirectory)//'.', exist=dirExists) +#else + inquire(file = trim(storeWorkingDirectory), exist=dirExists) +#endif + if(.not. dirExists) then ! check if the directory exists + write(6,'(a20,a,a16)') ' working directory "',trim(storeWorkingDirectory),'" does not exist' + call quit(1_pInt) + endif + else ! using path to geometry file as working dir + if (geometryArg(1:1) == pathSep) then ! absolute path given as command line argument + storeWorkingDirectory = geometryArg(1:scan(geometryArg,pathSep,back=.true.)) + else + error = getcwd(cwd) ! relative path given as command line argument + storeWorkingDirectory = trim(cwd)//pathSep//& + geometryArg(1:scan(geometryArg,pathSep,back=.true.)) + endif + endif + storeWorkingDirectory = rectifyPath(storeWorkingDirectory) + +end function storeWorkingDirectory + + +!-------------------------------------------------------------------------------------------------- +!> @brief simply returns the private string workingDir +!-------------------------------------------------------------------------------------------------- +character(len=1024) function getSolverWorkingDirectoryName() + + implicit none + getSolverWorkingDirectoryName = workingDirectory + +end function getSolverWorkingDirectoryName + + +!-------------------------------------------------------------------------------------------------- +!> @brief solver job name (no extension) as combination of geometry and load case name +!-------------------------------------------------------------------------------------------------- +character(len=1024) function getSolverJobName() + + implicit none + integer :: posExt,posSep + character :: pathSep + character(len=1024) :: tempString + + pathSep = getPathSep() + + tempString = geometryFile + posExt = scan(tempString,'.',back=.true.) + posSep = scan(tempString,pathSep,back=.true.) + + getSolverJobName = tempString(posSep+1:posExt-1) + + tempString = loadCaseFile + posExt = scan(tempString,'.',back=.true.) + posSep = scan(tempString,pathSep,back=.true.) + + getSolverJobName = trim(getSolverJobName)//'_'//tempString(posSep+1:posExt-1) + +end function getSolverJobName + + +!-------------------------------------------------------------------------------------------------- +!> @brief basename of geometry file with extension from command line arguments +!-------------------------------------------------------------------------------------------------- +character(len=1024) function getGeometryFile(geometryParameter) +#ifdef __INTEL_COMPILER + use IFPORT +#endif + + implicit none + character(len=1024), intent(in) :: & + geometryParameter + character(len=1024) :: & + cwd + integer :: posExt, posSep + character :: pathSep + integer :: error + + getGeometryFile = geometryParameter + pathSep = getPathSep() + posExt = scan(getGeometryFile,'.',back=.true.) + posSep = scan(getGeometryFile,pathSep,back=.true.) + + if (posExt <= posSep) getGeometryFile = trim(getGeometryFile)//('.geom') ! no extension present + if (scan(getGeometryFile,pathSep) /= 1) then ! relative path given as command line argument + error = getcwd(cwd) + getGeometryFile = rectifyPath(trim(cwd)//pathSep//getGeometryFile) + else + getGeometryFile = rectifyPath(getGeometryFile) + endif + + getGeometryFile = makeRelativePath(getSolverWorkingDirectoryName(), getGeometryFile) + +end function getGeometryFile + + +!-------------------------------------------------------------------------------------------------- +!> @brief relative path of loadcase from command line arguments +!-------------------------------------------------------------------------------------------------- +character(len=1024) function getLoadCaseFile(loadCaseParameter) +#ifdef __INTEL_COMPILER + use IFPORT +#endif + + implicit none + character(len=1024), intent(in) :: & + loadCaseParameter + character(len=1024) :: & + cwd + integer :: posExt, posSep, error + character :: pathSep + + getLoadCaseFile = loadcaseParameter + pathSep = getPathSep() + posExt = scan(getLoadCaseFile,'.',back=.true.) + posSep = scan(getLoadCaseFile,pathSep,back=.true.) + + if (posExt <= posSep) getLoadCaseFile = trim(getLoadCaseFile)//('.load') ! no extension present + if (scan(getLoadCaseFile,pathSep) /= 1) then ! relative path given as command line argument + error = getcwd(cwd) + getLoadCaseFile = rectifyPath(trim(cwd)//pathSep//getLoadCaseFile) + else + getLoadCaseFile = rectifyPath(getLoadCaseFile) + endif + + getLoadCaseFile = makeRelativePath(getSolverWorkingDirectoryName(), getLoadCaseFile) + +end function getLoadCaseFile + + +!-------------------------------------------------------------------------------------------------- +!> @brief remove ../ and /./ from path +!-------------------------------------------------------------------------------------------------- +function rectifyPath(path) + + implicit none + character(len=*) :: path + character(len=len_trim(path)) :: rectifyPath + character :: pathSep + integer :: i,j,k,l ! no pInt + + pathSep = getPathSep() + +!-------------------------------------------------------------------------------------------------- +! remove /./ from path + l = len_trim(path) + rectifyPath = path + do i = l,3,-1 + if (rectifyPath(i-2:i) == pathSep//'.'//pathSep) & + rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' + enddo + +!-------------------------------------------------------------------------------------------------- +! remove ../ and corresponding directory from rectifyPath + l = len_trim(rectifyPath) + i = index(rectifyPath(i:l),'..'//pathSep) + j = 0 + do while (i > j) + j = scan(rectifyPath(1:i-2),pathSep,back=.true.) + rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j) + if (rectifyPath(j+1:j+1) == pathSep) then !search for '//' that appear in case of XXX/../../XXX + k = len_trim(rectifyPath) + rectifyPath(j+1:k-1) = rectifyPath(j+2:k) + rectifyPath(k:k) = ' ' + endif + i = j+index(rectifyPath(j+1:l),'..'//pathSep) + enddo + if(len_trim(rectifyPath) == 0) rectifyPath = pathSep + +end function rectifyPath + + +!-------------------------------------------------------------------------------------------------- +!> @brief relative path from absolute a to absolute b +!-------------------------------------------------------------------------------------------------- +character(len=1024) function makeRelativePath(a,b) + + implicit none + character (len=*) :: a,b + character :: pathSep + integer :: i,posLastCommonSlash,remainingSlashes !no pInt + + pathSep = getPathSep() + posLastCommonSlash = 0 + remainingSlashes = 0 + + do i = 1, min(1024,len_trim(a),len_trim(b)) + if (a(i:i) /= b(i:i)) exit + if (a(i:i) == pathSep) posLastCommonSlash = i + enddo + do i = posLastCommonSlash+1,len_trim(a) + if (a(i:i) == pathSep) remainingSlashes = remainingSlashes + 1 + enddo + makeRelativePath = repeat('..'//pathSep,remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) + +end function makeRelativePath + + +!-------------------------------------------------------------------------------------------------- +!> @brief counting / and \ in $PATH System variable the character occuring more often is assumed +! to be the path separator +!-------------------------------------------------------------------------------------------------- +character function getPathSep() + + implicit none + character(len=2048) :: & + path + integer(pInt) :: & + backslash = 0_pInt, & + slash = 0_pInt + integer :: i + + call get_environment_variable('PATH',path) + do i=1, len(trim(path)) + if (path(i:i)=='/') slash = slash + 1_pInt + if (path(i:i)=='\') backslash = backslash + 1_pInt + enddo + + if (backslash>slash) then + getPathSep = '\' + else + getPathSep = '/' + endif + +end function getPathSep + + +!-------------------------------------------------------------------------------------------------- +!> @brief taken from IO, check IO_stringValue for documentation +!-------------------------------------------------------------------------------------------------- +pure function IIO_stringValue(string,chunkPos,myChunk) + + implicit none + integer(pInt), dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string + integer(pInt), intent(in) :: myChunk !< position number of desired chunk + character(len=1+chunkPos(myChunk*2+1)-chunkPos(myChunk*2)) :: IIO_stringValue + character(len=*), intent(in) :: string !< raw input with known start and end of each chunk + + + valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then + IIO_stringValue = '' + else valuePresent + IIO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)) + endif valuePresent + +end function IIO_stringValue + + +!-------------------------------------------------------------------------------------------------- +!> @brief taken from IO, check IO_intValue for documentation +!-------------------------------------------------------------------------------------------------- +integer(pInt) pure function IIO_intValue(string,chunkPos,myChunk) + + implicit none + character(len=*), intent(in) :: string !< raw input with known start and end of each chunk + integer(pInt), intent(in) :: myChunk !< position number of desired sub string + integer(pInt), dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string + + + valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then + IIO_intValue = 0_pInt + else valuePresent + read(UNIT=string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)),ERR=100,FMT=*) IIO_intValue + endif valuePresent + return +100 IIO_intValue = huge(1_pInt) + +end function IIO_intValue + + +!-------------------------------------------------------------------------------------------------- +!> @brief taken from IO, check IO_lc for documentation +!-------------------------------------------------------------------------------------------------- +pure function IIO_lc(string) + + implicit none + character(len=*), intent(in) :: string !< string to convert + character(len=len(string)) :: IIO_lc + + character(26), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz' + character(26), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' + + integer :: i,n ! no pInt (len returns default integer) + + IIO_lc = string + do i=1,len(string) + n = index(UPPER,IIO_lc(i:i)) + if (n/=0) IIO_lc(i:i) = LOWER(n:n) + enddo + +end function IIO_lc + + +!-------------------------------------------------------------------------------------------------- +!> @brief taken from IO, check IO_stringPos for documentation +!-------------------------------------------------------------------------------------------------- +pure function IIO_stringPos(string) + + implicit none + integer(pInt), dimension(:), allocatable :: IIO_stringPos + character(len=*), intent(in) :: string !< string in which chunks are searched for + + character(len=*), parameter :: SEP=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces + integer :: left, right ! no pInt (verify and scan return default integer) + + allocate(IIO_stringPos(1), source=0_pInt) + right = 0 + + do while (verify(string(right+1:),SEP)>0) + left = right + verify(string(right+1:),SEP) + right = left + scan(string(left:),SEP) - 2 + if ( string(left:left) == '#' ) exit + IIO_stringPos = [IIO_stringPos,int(left, pInt), int(right, pInt)] + IIO_stringPos(1) = IIO_stringPos(1)+1_pInt + enddo + +end function IIO_stringPos + + +end module diff --git a/code/spectral/spectral_mech_AL.f90 b/code/spectral/spectral_mech_AL.f90 new file mode 100644 index 000000000..a937dcc86 --- /dev/null +++ b/code/spectral/spectral_mech_AL.f90 @@ -0,0 +1,715 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @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 +!-------------------------------------------------------------------------------------------------- +module spectral_mech_AL + use prec, only: & + pInt, & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private +#include + + character (len=*), parameter, public :: & + DAMASK_spectral_solverAL_label = 'al' + +!-------------------------------------------------------------------------------------------------- +! derived types + type(tSolutionParams), private :: params + real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal + +!-------------------------------------------------------------------------------------------------- +! PETSc data + DM, private :: da + SNES, private :: snes + Vec, private :: solution_vec + +!-------------------------------------------------------------------------------------------------- +! common pointwise data + real(pReal), private, dimension(:,:,:,:,:), allocatable :: & + F_lastInc, & !< field of previous compatible deformation gradients + F_lambda_lastInc, & !< field of previous incompatible deformation gradient + Fdot, & !< field of assumed rate of compatible deformation gradient + F_lambdaDot !< field of assumed rate of incopatible deformation gradient + +!-------------------------------------------------------------------------------------------------- +! stress, stiffness and compliance average etc. + real(pReal), private, dimension(3,3) :: & + 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 + F_av = 0.0_pReal, & !< average incompatible def grad field + 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 + character(len=1024), private :: incInfo !< time and increment information + real(pReal), private, dimension(3,3,3,3) :: & + 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 + S = 0.0_pReal, & !< current compliance (filled up with zeros) + C_scale = 0.0_pReal, & + S_scale = 0.0_pReal + + real(pReal), private :: & + err_BC, & !< deviation from stress BC + err_curl, & !< RMS of curl of F + err_div !< RMS of div of P + logical, private :: ForwardData + integer(pInt), private :: & + totalIter = 0_pInt !< total iteration in current increment + + public :: & + AL_init, & + AL_solution, & + AL_forward, & + AL_destroy + external :: & + VecDestroy, & + DMDestroy, & + DMDACreate3D, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + PETScFinalize, & + SNESDestroy, & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber, & + SNESSolve, & + SNESSetDM, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions, & + SNESCreate, & + MPI_Abort, & + MPI_Bcast, & + MPI_Allreduce + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc) +!-------------------------------------------------------------------------------------------------- +subroutine AL_init + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & + IO_intOut, & + IO_read_realFile, & + IO_timeStamp + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRestart + use FEsolving, only: & + restartInc + use numerics, only: & + worldrank, & + worldsize + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + Utilities_constitutiveResponse, & + Utilities_updateGamma, & + Utilities_updateIPcoords + use mesh, only: & + grid, & + grid3 + use math, only: & + math_invSym3333 + + implicit none + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal + + PetscErrorCode :: ierr + PetscObject :: dummy + PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda + integer(pInt), dimension(:), allocatable :: localK + integer(pInt) :: proc + character(len=1024) :: rankStr + + if (worldrank == 0_pInt) then + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif + +!-------------------------------------------------------------------------------------------------- +! allocate global fields + allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_lambda_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_lambdaDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! PETSc Init + call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 + do proc = 1, worldsize + call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) + enddo + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary + DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & ! global grid + 1 , 1, worldsize, & + 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) + grid(1),grid(2),localK, & ! local grid + da,ierr) ! handle, error + CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,dummy,ierr) + CHKERRQ(ierr) + call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr) + CHKERRQ(ierr) + call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! init fields + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data + F => xx_psc(0:8,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) + restart: if (restartInc > 1_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading values of increment ', restartInc - 1_pInt, ' from file' + flush(6) + write(rankStr,'(a1,i0)')'_',worldrank + call IO_read_realFile(777,'F'//trim(rankStr), trim(getSolverJobName()),size(F)) + read (777,rec=1) F + close (777) + call IO_read_realFile(777,'F_lastInc'//trim(rankStr), trim(getSolverJobName()),size(F_lastInc)) + read (777,rec=1) F_lastInc + close (777) + call IO_read_realFile(777,'F_lambda'//trim(rankStr),trim(getSolverJobName()),size(F_lambda)) + read (777,rec=1) F_lambda + close (777) + call IO_read_realFile(777,'F_lambda_lastInc'//trim(rankStr),& + trim(getSolverJobName()),size(F_lambda_lastInc)) + read (777,rec=1) F_lambda_lastInc + close (777) + call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim)) + read (777,rec=1) F_aim + close (777) + call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc)) + read (777,rec=1) F_aim_lastInc + close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) + read (777,rec=1) f_aimDot + close (777) + elseif (restartInc == 1_pInt) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_lambda = F + F_lambda_lastInc = F_lastInc + endif restart + + + call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) + call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & + 0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) + nullify(F) + nullify(F_lambda) + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc + + readRestart: if (restartInc > 1_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading more values of increment', restartInc - 1_pInt, 'from file' + flush(6) + call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) + read (777,rec=1) C_volAvg + close (777) + call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) + read (777,rec=1) C_volAvgLastInc + close (777) + call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) + read (777,rec=1) C_minMaxAvg + close (777) + endif readRestart + + call Utilities_updateGamma(C_minMaxAvg,.True.) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) + +end subroutine AL_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief solution for the AL scheme with internal iterations +!-------------------------------------------------------------------------------------------------- +type(tSolutionState) function & + AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,rotation_BC) + use IO, only: & + IO_error + use numerics, only: & + update_gamma + use math, only: & + math_invSym3333 + use spectral_utilities, only: & + tBoundaryCondition, & + Utilities_maskedCompliance, & + Utilities_updateGamma + use FEsolving, only: & + restartWrite, & + terminallyIll + + implicit none + +!-------------------------------------------------------------------------------------------------- +! input data for solution + 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 + logical, intent(in) :: & + guess + type(tBoundaryCondition), intent(in) :: & + P_BC, & + F_BC + character(len=*), intent(in) :: & + incInfoIn + real(pReal), dimension(3,3), intent(in) :: rotation_BC + +!-------------------------------------------------------------------------------------------------- +! PETSc Data + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn + +!-------------------------------------------------------------------------------------------------- +! update stiffness (and gamma operator) + S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) + if (update_gamma) then + call Utilities_updateGamma(C_minMaxAvg,restartWrite) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) + endif + +!-------------------------------------------------------------------------------------------------- +! set module wide availabe data + mask_stress = P_BC%maskFloat + params%P_BC = P_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + +!-------------------------------------------------------------------------------------------------- +! solve BVP + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) + CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! check convergence + call SNESGetConvergedReason(snes,reason,ierr) + CHKERRQ(ierr) + AL_solution%termIll = terminallyIll + terminallyIll = .false. + if (reason == -4) call IO_error(893_pInt) + if (reason < 1) AL_solution%converged = .false. + AL_solution%iterationsNeeded = totalIter + +end function AL_solution + + +!-------------------------------------------------------------------------------------------------- +!> @brief forms the AL residual vector +!-------------------------------------------------------------------------------------------------- +subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + polarAlpha, & + polarBeta, & + worldrank + use mesh, only: & + grid3, & + grid + use IO, only: & + IO_intOut + use math, only: & + math_rotate_backward33, & + math_transpose33, & + math_mul3333xx33, & + math_invSym3333, & + math_mul33x33 + use spectral_utilities, only: & + wgt, & + tensorField_real, & + utilities_FFTtensorForward, & + utilities_fourierGammaConvolution, & + utilities_FFTtensorBackward, & + Utilities_constitutiveResponse, & + Utilities_divergenceRMS, & + Utilities_curlRMS + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use homogenization, only: & + materialpoint_dPdF + use FEsolving, only: & + terminallyIll + + implicit none +!-------------------------------------------------------------------------------------------------- +! strange syntax in the next line because otherwise macros expand beyond 132 character limit + DMDALocalInfo, dimension(& + DMDA_LOCAL_INFO_SIZE) :: & + in + PetscScalar, target, dimension(3,3,2, & + XG_RANGE,YG_RANGE,ZG_RANGE) :: & + x_scal + PetscScalar, target, dimension(3,3,2, & + X_RANGE,Y_RANGE,Z_RANGE) :: & + f_scal + PetscScalar, pointer, dimension(:,:,:,:,:) :: & + F, & + F_lambda, & + residual_F, & + residual_F_lambda + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + integer(pInt) :: & + i, j, k, e + + F => x_scal(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_lambda => x_scal(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => f_scal(1:3,1:3,1,& + X_RANGE,Y_RANGE,Z_RANGE) + residual_F_lambda => f_scal(1:3,1:3,2,& + X_RANGE,Y_RANGE,Z_RANGE) + + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + + F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt + call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + + if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment + newIteration: if(totalIter <= PETScIter) then +!-------------------------------------------------------------------------------------------------- +! report begin of new iteration + totalIter = totalIter + 1_pInt + if (worldrank == 0_pInt) then + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & + math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & + math_transpose33(F_aim) + flush(6) + endif + endif newIteration + +!-------------------------------------------------------------------------------------------------- +! + tensorField_real = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + tensorField_real(1:3,1:3,i,j,k) = & + polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3)) + + enddo; enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! doing convolution in Fourier space + call utilities_FFTtensorForward() + call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call utilities_FFTtensorBackward() + +!-------------------------------------------------------------------------------------------------- +! constructing residual + residual_F_lambda = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + +!-------------------------------------------------------------------------------------------------- +! evaluate constitutive response + P_avLastEval = P_av + call Utilities_constitutiveResponse(F_lastInc,F - residual_F_lambda/polarBeta,params%timeinc, & + residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + ForwardData = .False. + +!-------------------------------------------------------------------------------------------------- +! calculate divergence + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F + call utilities_FFTtensorForward() + err_div = Utilities_divergenceRMS() + call utilities_FFTtensorBackward() + +!-------------------------------------------------------------------------------------------------- +! constructing residual + e = 0_pInt + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + e = e + 1_pInt + residual_F(1:3,1:3,i,j,k) = math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & + residual_F(1:3,1:3,i,j,k) - & + math_mul33x33(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3))) & + + residual_F_lambda(1:3,1:3,i,j,k) + enddo; enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! calculating curl + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + call utilities_FFTtensorForward() + err_curl = Utilities_curlRMS() + call utilities_FFTtensorBackward() + +end subroutine AL_formResidual + + +!-------------------------------------------------------------------------------------------------- +!> @brief convergence check +!-------------------------------------------------------------------------------------------------- +subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_curl_tolRel, & + err_curl_tolAbs, & + err_stress_tolAbs, & + err_stress_tolRel, & + worldrank + use math, only: & + math_mul3333xx33 + use FEsolving, only: & + terminallyIll + + implicit none + SNES :: snes_local + PetscInt :: PETScIter + PetscReal :: & + xnorm, & + snorm, & + fnorm + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode ::ierr + real(pReal) :: & + curlTol, & + divTol, & + BC_tol + +!-------------------------------------------------------------------------------------------------- +! stress BC handling + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc + err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & + mask_stress *(P_av - params%P_BC))) ! mask = 0.0 for no bc + +!-------------------------------------------------------------------------------------------------- +! error calculation + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel,err_div_tolAbs) + BC_tol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) + + converged: if ((totalIter >= itmin .and. & + all([ err_div/divTol, & + err_curl/curlTol, & + err_BC/BC_tol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then converged + reason = -1 + else converged + reason = 0 + endif converged + +!-------------------------------------------------------------------------------------------------- +! report + if (worldrank == 0_pInt) then + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' + write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' + write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) + endif + +end subroutine AL_converged + +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) + use math, only: & + math_mul33x33, & + math_mul3333xx33, & + math_transpose33, & + math_rotate_backward33 + use numerics, only: & + worldrank + use mesh, only: & + grid3, & + grid + use spectral_utilities, only: & + Utilities_calculateRate, & + Utilities_forwardField, & + Utilities_updateIPcoords, & + tBoundaryCondition, & + cutBack + use IO, only: & + IO_write_JobRealFile + use FEsolving, only: & + restartWrite + + implicit none + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + P_BC, & + F_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + logical, intent(in) :: & + guess + PetscErrorCode :: ierr + PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda + integer(pInt) :: i, j, k + real(pReal), dimension(3,3) :: F_lambda33 + character(len=1024) :: rankStr + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) + F => xx_psc(0:8,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) + if (restartWrite) then + if (worldrank == 0_pInt) then + write(6,'(/,a)') ' writing converged results for restart' + flush(6) + endif + write(rankStr,'(a1,i0)')'_',worldrank + call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file + write (777,rec=1) F + close (777) + call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_lastInc + close (777) + call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file + write (777,rec=1) F_lambda + close (777) + call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_lambda_lastInc + close (777) + if (worldrank == 0_pInt) then + call IO_write_jobRealFile(777,'F_aim',size(F_aim)) + write (777,rec=1) F_aim + close(777) + call IO_write_jobRealFile(777,'F_aim_lastInc',size(F_aim_lastInc)) + write (777,rec=1) F_aim_lastInc + close(777) + call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) + write (777,rec=1) F_aimDot + close(777) + call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg)) + write (777,rec=1) C_volAvg + close(777) + call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) + write (777,rec=1) C_volAvgLastInc + close(777) + endif + endif + + call utilities_updateIPcoords(F) + + if (cutBack) then + F_aim = F_aim_lastInc + F_lambda = reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid3]) + F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) + C_volAvg = C_volAvgLastInc + else + ForwardData = .True. + C_volAvgLastInc = C_volAvg +!-------------------------------------------------------------------------------------------------- +! 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) + elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed + f_aimDot = F_BC%maskFloat * F_BC%values + 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 + endif + if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + F_aim_lastInc = F_aim + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + call utilities_updateIPcoords(F) + Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3])) + 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),grid3])) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) + F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid3]) + endif + + F_aim = F_aim + f_aimDot * timeinc + +!-------------------------------------------------------------------------------------------------- +! update local deformation gradient + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim + math_rotate_backward33(F_aim,rotation_BC)), & + [9,grid(1),grid(2),grid3]) + F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & + [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + if (.not. guess) then ! large strain forwarding + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + F_lambda33 = reshape(F_lambda(1:9,i,j,k),[3,3]) + F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & + math_mul3333xx33(C_scale,& + math_mul33x33(math_transpose33(F_lambda33),& + F_lambda33) -math_I3))*0.5_pReal)& + + math_I3 + F_lambda(1:9,i,j,k) = reshape(F_lambda33,[9]) + enddo; enddo; enddo + endif + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + +end subroutine AL_forward + +!-------------------------------------------------------------------------------------------------- +!> @brief destroy routine +!-------------------------------------------------------------------------------------------------- +subroutine AL_destroy() + use spectral_utilities, only: & + Utilities_destroy + + implicit none + PetscErrorCode :: ierr + + call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) + call SNESDestroy(snes,ierr); CHKERRQ(ierr) + call DMDestroy(da,ierr); CHKERRQ(ierr) + +end subroutine AL_destroy + +end module spectral_mech_AL diff --git a/code/spectral/spectral_mech_Basic.f90 b/code/spectral/spectral_mech_Basic.f90 new file mode 100644 index 000000000..a8344fabe --- /dev/null +++ b/code/spectral/spectral_mech_Basic.f90 @@ -0,0 +1,569 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @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 Basic scheme PETSc solver +!-------------------------------------------------------------------------------------------------- +module spectral_mech_basic + use prec, only: & + pInt, & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private +#include + + character (len=*), parameter, public :: & + DAMASK_spectral_SolverBasicPETSC_label = 'basicpetsc' + +!-------------------------------------------------------------------------------------------------- +! derived types + type(tSolutionParams), private :: params + +!-------------------------------------------------------------------------------------------------- +! PETSc data + DM, private :: da + SNES, private :: snes + Vec, private :: solution_vec + +!-------------------------------------------------------------------------------------------------- +! common pointwise data + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot + +!-------------------------------------------------------------------------------------------------- +! stress, stiffness and compliance average etc. + real(pReal), private, dimension(3,3) :: & + F_aim = math_I3, & + F_aim_lastIter = math_I3, & + F_aim_lastInc = math_I3, & + P_av = 0.0_pReal, & + F_aimDot=0.0_pReal + character(len=1024), private :: incInfo + real(pReal), private, dimension(3,3,3,3) :: & + 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 + S = 0.0_pReal !< current compliance (filled up with zeros) + real(pReal), private :: err_stress, err_div + logical, private :: ForwardData + integer(pInt), private :: & + totalIter = 0_pInt !< total iteration in current increment + real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal + + public :: & + basicPETSc_init, & + basicPETSc_solution, & + BasicPETSc_forward, & + basicPETSc_destroy + external :: & + VecDestroy, & + DMDestroy, & + DMDACreate3D, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + PETScFinalize, & + SNESDestroy, & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber, & + SNESSolve, & + SNESSetDM, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions, & + SNESCreate, & + MPI_Abort, & + MPI_Bcast, & + MPI_Allreduce + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!-------------------------------------------------------------------------------------------------- +subroutine basicPETSc_init + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & + IO_intOut, & + IO_read_realFile, & + IO_timeStamp + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRestart + use FEsolving, only: & + restartInc + use numerics, only: & + worldrank, & + worldsize + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + Utilities_constitutiveResponse, & + Utilities_updateGamma, & + utilities_updateIPcoords, & + wgt + use mesh, only: & + grid, & + grid3 + use math, only: & + math_invSym3333 + + implicit none + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + PetscScalar, dimension(:,:,:,:), pointer :: F + PetscErrorCode :: ierr + PetscObject :: dummy + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal + integer(pInt), dimension(:), allocatable :: localK + integer(pInt) :: proc + character(len=1024) :: rankStr + + mainProcess: if (worldrank == 0_pInt) then + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + +!-------------------------------------------------------------------------------------------------- +! allocate global fields + allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! initialize solver specific parts of PETSc + call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 + do proc = 1, worldsize + call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) + enddo + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary + DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & ! global grid + 1, 1, worldsize, & + 9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) + grid (1),grid (2),localK, & ! local grid + da,ierr) ! handle, error + CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,dummy,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da + call SNESSetConvergenceTest(snes,BasicPETSC_converged,dummy,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" + CHKERRQ(ierr) + call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments + +!-------------------------------------------------------------------------------------------------- +! init fields + call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with + + restart: if (restartInc > 1_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading values of increment ', restartInc - 1_pInt, ' from file' + flush(6) + write(rankStr,'(a1,i0)')'_',worldrank + call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) + read (777,rec=1) F + close (777) + call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) + read (777,rec=1) F_lastInc + close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) + read (777,rec=1) f_aimDot + close (777) + 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 + elseif (restartInc == 1_pInt) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + endif restart + + call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) + call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & + 0.0_pReal, & + P, & + C_volAvg,C_minMaxAvg, & ! global average of stiffness and (min+max)/2 + temp33_Real, & + .false., & + math_I3) + + call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc + + restartRead: if (restartInc > 1_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading more values of increment', restartInc - 1_pInt, 'from file' + flush(6) + call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) + read (777,rec=1) C_volAvg + close (777) + call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) + read (777,rec=1) C_volAvgLastInc + close (777) + call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) + read (777,rec=1) C_minMaxAvg + close (777) + endif restartRead + + call Utilities_updateGamma(C_minmaxAvg,.True.) + +end subroutine basicPETSc_init + +!-------------------------------------------------------------------------------------------------- +!> @brief solution for the Basic PETSC scheme with internal iterations +!-------------------------------------------------------------------------------------------------- +type(tSolutionState) function & + basicPETSc_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,rotation_BC) + use IO, only: & + IO_error + use numerics, only: & + update_gamma + use spectral_utilities, only: & + tBoundaryCondition, & + Utilities_maskedCompliance, & + Utilities_updateGamma + use FEsolving, only: & + restartWrite, & + terminallyIll + + implicit none + +!-------------------------------------------------------------------------------------------------- +! input data for solution + 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 + type(tBoundaryCondition), intent(in) :: & + P_BC, & + F_BC + character(len=*), intent(in) :: & + incInfoIn + real(pReal), dimension(3,3), intent(in) :: rotation_BC + logical, intent(in) :: & + guess + +!-------------------------------------------------------------------------------------------------- +! PETSc Data + PetscErrorCode :: ierr + SNESConvergedReason :: reason + incInfo = incInfoIn + +!-------------------------------------------------------------------------------------------------- +! update stiffness (and gamma operator) + S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) + if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite) + + +!-------------------------------------------------------------------------------------------------- +! set module wide availabe data + mask_stress = P_BC%maskFloat + params%P_BC = P_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + +!-------------------------------------------------------------------------------------------------- +! solve BVP + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) + CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! check convergence + call SNESGetConvergedReason(snes,reason,ierr) + CHKERRQ(ierr) + basicPETSc_solution%termIll = terminallyIll + terminallyIll = .false. + BasicPETSc_solution%converged =.true. + if (reason == -4) call IO_error(893_pInt) + if (reason < 1) basicPETSC_solution%converged = .false. + basicPETSC_solution%iterationsNeeded = totalIter + +end function BasicPETSc_solution + + +!-------------------------------------------------------------------------------------------------- +!> @brief forms the AL residual vector +!-------------------------------------------------------------------------------------------------- +subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) + use numerics, only: & + itmax, & + itmin + use numerics, only: & + worldrank + use mesh, only: & + grid, & + grid3 + use math, only: & + math_rotate_backward33, & + math_transpose33, & + math_mul3333xx33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use spectral_utilities, only: & + tensorField_real, & + utilities_FFTtensorForward, & + utilities_FFTtensorBackward, & + utilities_fourierGammaConvolution, & + Utilities_constitutiveResponse, & + Utilities_divergenceRMS + use IO, only: & + IO_intOut + use FEsolving, only: & + terminallyIll + + implicit none + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & + in + PetscScalar, dimension(3,3, & + XG_RANGE,YG_RANGE,ZG_RANGE) :: & + x_scal + PetscScalar, dimension(3,3, & + X_RANGE,Y_RANGE,Z_RANGE) :: & + f_scal + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + + if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment + newIteration: if (totalIter <= PETScIter) then +!-------------------------------------------------------------------------------------------------- +! report begin of new iteration + totalIter = totalIter + 1_pInt + if (worldrank == 0_pInt) then + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & + math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & + math_transpose33(F_aim) + flush(6) + endif + endif newIteration + +!-------------------------------------------------------------------------------------------------- +! evaluate constitutive response + call Utilities_constitutiveResponse(F_lastInc,x_scal,params%timeinc, & + f_scal,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + ForwardData = .false. + +!-------------------------------------------------------------------------------------------------- +! stress BC handling + F_aim_lastIter = F_aim + F_aim = F_aim - 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 + +!-------------------------------------------------------------------------------------------------- +! updated deformation gradient using fix point algorithm of basic scheme + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = f_scal + call utilities_FFTtensorForward() + err_div = Utilities_divergenceRMS() + call utilities_fourierGammaConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC)) + call utilities_FFTtensorBackward() + +!-------------------------------------------------------------------------------------------------- +! constructing residual + f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + +end subroutine BasicPETSc_formResidual + + +!-------------------------------------------------------------------------------------------------- +!> @brief convergence check +!-------------------------------------------------------------------------------------------------- +subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_stress_tolRel, & + err_stress_tolAbs, & + worldrank + use FEsolving, only: & + terminallyIll + + implicit none + SNES :: snes_local + PetscInt :: PETScIter + PetscReal :: & + xnorm, & + snorm, & + fnorm + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal) :: & + divTol, & + stressTol + + divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs) + stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) + + converged: if ((totalIter >= itmin .and. & + all([ err_div/divTol, & + err_stress/stressTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then converged + reason = -1 + else converged + reason = 0 + endif converged + +!-------------------------------------------------------------------------------------------------- +! report + if (worldrank == 0_pInt) then + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) + endif + +end subroutine BasicPETSc_converged + +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) + use math, only: & + math_mul33x33 ,& + math_rotate_backward33 + use mesh, only: & + grid, & + grid3 + use spectral_utilities, only: & + Utilities_calculateRate, & + Utilities_forwardField, & + utilities_updateIPcoords, & + tBoundaryCondition, & + cutBack + use IO, only: & + IO_write_JobRealFile + use FEsolving, only: & + restartWrite + use numerics, only: & + worldrank + + implicit none + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + P_BC, & + F_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + logical, intent(in) :: & + guess + PetscScalar, pointer :: F(:,:,:,:) + PetscErrorCode :: ierr + character(len=1024) :: rankStr + + call DMDAVecGetArrayF90(da,solution_vec,F,ierr) +!-------------------------------------------------------------------------------------------------- +! restart information for spectral solver + if (restartWrite) then + if (worldrank == 0_pInt) then + write(6,'(/,a)') ' writing converged results for restart' + flush(6) + endif + write(rankStr,'(a1,i0)')'_',worldrank + call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file + write (777,rec=1) F + close (777) + call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_lastInc + close (777) + if (worldrank == 0_pInt) then + call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) + write (777,rec=1) F_aimDot + close(777) + call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg)) + write (777,rec=1) C_volAvg + close(777) + call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) + write (777,rec=1) C_volAvgLastInc + close(777) + endif + endif + + call utilities_updateIPcoords(F) + + if (cutBack) then + F_aim = F_aim_lastInc + F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) + C_volAvg = C_volAvgLastInc + else + ForwardData = .True. + C_volAvgLastInc = C_volAvg +!-------------------------------------------------------------------------------------------------- +! 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) + elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed + f_aimDot = F_BC%maskFloat * F_BC%values + 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 + endif + if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + F_aim_lastInc = F_aim + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + call utilities_updateIPcoords(F) + Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3])) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) + endif + F_aim = F_aim + f_aimDot * timeinc + +!-------------------------------------------------------------------------------------------------- +! update local deformation gradient + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim + math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) + call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) + +end subroutine BasicPETSc_forward + +!-------------------------------------------------------------------------------------------------- +!> @brief destroy routine +!-------------------------------------------------------------------------------------------------- +subroutine BasicPETSc_destroy() + use spectral_utilities, only: & + Utilities_destroy + + implicit none + PetscErrorCode :: ierr + + call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) + call SNESDestroy(snes,ierr); CHKERRQ(ierr) + call DMDestroy(da,ierr); CHKERRQ(ierr) + +end subroutine BasicPETSc_destroy + +end module spectral_mech_basic diff --git a/code/spectral/spectral_mech_Polarisation.f90 b/code/spectral/spectral_mech_Polarisation.f90 new file mode 100644 index 000000000..a28eb5adb --- /dev/null +++ b/code/spectral/spectral_mech_Polarisation.f90 @@ -0,0 +1,712 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @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 Polarisation scheme solver +!-------------------------------------------------------------------------------------------------- +module spectral_mech_Polarisation + use prec, only: & + pInt, & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private +#include + + character (len=*), parameter, public :: & + DAMASK_spectral_solverPolarisation_label = 'polarisation' + +!-------------------------------------------------------------------------------------------------- +! derived types + type(tSolutionParams), private :: params + real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal + +!-------------------------------------------------------------------------------------------------- +! PETSc data + DM, private :: da + SNES, private :: snes + Vec, private :: solution_vec + +!-------------------------------------------------------------------------------------------------- +! common pointwise data + real(pReal), private, dimension(:,:,:,:,:), allocatable :: & + F_lastInc, & !< field of previous compatible deformation gradients + F_tau_lastInc, & !< field of previous incompatible deformation gradient + Fdot, & !< field of assumed rate of compatible deformation gradient + F_tauDot !< field of assumed rate of incopatible deformation gradient + +!-------------------------------------------------------------------------------------------------- +! stress, stiffness and compliance average etc. + real(pReal), private, dimension(3,3) :: & + 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 + F_av = 0.0_pReal, & !< average incompatible def grad field + 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 + character(len=1024), private :: incInfo !< time and increment information + real(pReal), private, dimension(3,3,3,3) :: & + 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 + S = 0.0_pReal, & !< current compliance (filled up with zeros) + C_scale = 0.0_pReal, & + S_scale = 0.0_pReal + + real(pReal), private :: & + err_BC, & !< deviation from stress BC + err_curl, & !< RMS of curl of F + err_div !< RMS of div of P + logical, private :: ForwardData + integer(pInt), private :: & + totalIter = 0_pInt !< total iteration in current increment + + public :: & + Polarisation_init, & + Polarisation_solution, & + Polarisation_forward, & + Polarisation_destroy + external :: & + VecDestroy, & + DMDestroy, & + DMDACreate3D, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + PETScFinalize, & + SNESDestroy, & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber, & + SNESSolve, & + SNESSetDM, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions, & + SNESCreate, & + MPI_Abort, & + MPI_Bcast, & + MPI_Allreduce + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc) +!-------------------------------------------------------------------------------------------------- +subroutine Polarisation_init + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & + IO_intOut, & + IO_read_realFile, & + IO_timeStamp + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRestart + use FEsolving, only: & + restartInc + use numerics, only: & + worldrank, & + worldsize + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + Utilities_constitutiveResponse, & + Utilities_updateGamma, & + Utilities_updateIPcoords + use mesh, only: & + grid, & + grid3 + use math, only: & + math_invSym3333 + + implicit none + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal + + PetscErrorCode :: ierr + PetscObject :: dummy + PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_tau + integer(pInt), dimension(:), allocatable :: localK + integer(pInt) :: proc + character(len=1024) :: rankStr + + if (worldrank == 0_pInt) then + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif + +!-------------------------------------------------------------------------------------------------- +! allocate global fields + allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! PETSc Init + call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 + do proc = 1, worldsize + call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) + enddo + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary + DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & ! global grid + 1 , 1, worldsize, & + 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) + grid (1),grid (2),localK, & ! local grid + da,ierr) ! handle, error + CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,dummy,ierr) + CHKERRQ(ierr) + call SNESSetConvergenceTest(snes,Polarisation_converged,dummy,PETSC_NULL_FUNCTION,ierr) + CHKERRQ(ierr) + call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! init fields + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data + F => xx_psc(0:8,:,:,:) + F_tau => xx_psc(9:17,:,:,:) + restart: if (restartInc > 1_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading values of increment', restartInc - 1_pInt, 'from file' + flush(6) + write(rankStr,'(a1,i0)')'_',worldrank + call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) + read (777,rec=1) F + close (777) + call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) + read (777,rec=1) F_lastInc + close (777) + call IO_read_realFile(777,'F_tau'//trim(rankStr),trim(getSolverJobName()),size(F_tau)) + read (777,rec=1) F_tau + close (777) + call IO_read_realFile(777,'F_tau_lastInc'//trim(rankStr),& + trim(getSolverJobName()),size(F_tau_lastInc)) + read (777,rec=1) F_tau_lastInc + close (777) + call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim)) + read (777,rec=1) F_aim + close (777) + call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc)) + read (777,rec=1) F_aim_lastInc + close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) + read (777,rec=1) f_aimDot + close (777) + elseif (restartInc == 1_pInt) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_tau = 2.0_pReal* F + F_tau_lastInc = 2.0_pReal*F_lastInc + endif restart + + call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) + call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & + 0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) + nullify(F) + nullify(F_tau) + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc + + readRestart: if (restartInc > 1_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & + write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & + 'reading more values of increment', restartInc - 1_pInt, 'from file' + flush(6) + call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) + read (777,rec=1) C_volAvg + close (777) + call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) + read (777,rec=1) C_volAvgLastInc + close (777) + call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) + read (777,rec=1) C_minMaxAvg + close (777) + endif readRestart + + call Utilities_updateGamma(C_minMaxAvg,.True.) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) + +end subroutine Polarisation_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief solution for the Polarisation scheme with internal iterations +!-------------------------------------------------------------------------------------------------- +type(tSolutionState) function & + Polarisation_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,rotation_BC) + use IO, only: & + IO_error + use numerics, only: & + update_gamma + use math, only: & + math_invSym3333 + use spectral_utilities, only: & + tBoundaryCondition, & + Utilities_maskedCompliance, & + Utilities_updateGamma + use FEsolving, only: & + restartWrite, & + terminallyIll + + implicit none + +!-------------------------------------------------------------------------------------------------- +! input data for solution + 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 + logical, intent(in) :: & + guess + type(tBoundaryCondition), intent(in) :: & + P_BC, & + F_BC + character(len=*), intent(in) :: & + incInfoIn + real(pReal), dimension(3,3), intent(in) :: rotation_BC + +!-------------------------------------------------------------------------------------------------- +! PETSc Data + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn + +!-------------------------------------------------------------------------------------------------- +! update stiffness (and gamma operator) + S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) + if (update_gamma) then + call Utilities_updateGamma(C_minMaxAvg,restartWrite) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) + endif + +!-------------------------------------------------------------------------------------------------- +! set module wide availabe data + mask_stress = P_BC%maskFloat + params%P_BC = P_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + +!-------------------------------------------------------------------------------------------------- +! solve BVP + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) + CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! check convergence + call SNESGetConvergedReason(snes,reason,ierr) + CHKERRQ(ierr) + Polarisation_solution%termIll = terminallyIll + terminallyIll = .false. + if (reason == -4) call IO_error(893_pInt) + if (reason < 1) Polarisation_solution%converged = .false. + Polarisation_solution%iterationsNeeded = totalIter + +end function Polarisation_solution + + +!-------------------------------------------------------------------------------------------------- +!> @brief forms the Polarisation residual vector +!-------------------------------------------------------------------------------------------------- +subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + polarAlpha, & + polarBeta, & + worldrank + use mesh, only: & + grid3, & + grid + use IO, only: & + IO_intOut + use math, only: & + math_rotate_backward33, & + math_transpose33, & + math_mul3333xx33, & + math_invSym3333, & + math_mul33x33 + use spectral_utilities, only: & + wgt, & + tensorField_real, & + utilities_FFTtensorForward, & + utilities_fourierGammaConvolution, & + utilities_FFTtensorBackward, & + Utilities_constitutiveResponse, & + Utilities_divergenceRMS, & + Utilities_curlRMS + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use homogenization, only: & + materialpoint_dPdF + use FEsolving, only: & + terminallyIll + + implicit none +!-------------------------------------------------------------------------------------------------- +! strange syntax in the next line because otherwise macros expand beyond 132 character limit + DMDALocalInfo, dimension(& + DMDA_LOCAL_INFO_SIZE) :: & + in + PetscScalar, target, dimension(3,3,2, & + XG_RANGE,YG_RANGE,ZG_RANGE) :: & + x_scal + PetscScalar, target, dimension(3,3,2, & + X_RANGE,Y_RANGE,Z_RANGE) :: & + f_scal + PetscScalar, pointer, dimension(:,:,:,:,:) :: & + F, & + F_tau, & + residual_F, & + residual_F_tau + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + integer(pInt) :: & + i, j, k, e + + F => x_scal(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_tau => x_scal(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => f_scal(1:3,1:3,1,& + X_RANGE,Y_RANGE,Z_RANGE) + residual_F_tau => f_scal(1:3,1:3,2,& + X_RANGE,Y_RANGE,Z_RANGE) + + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + + F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt + call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + + if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment + newIteration: if(totalIter <= PETScIter) then +!-------------------------------------------------------------------------------------------------- +! report begin of new iteration + totalIter = totalIter + 1_pInt + if (worldrank == 0_pInt) then + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & + math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & + math_transpose33(F_aim) + flush(6) + endif + endif newIteration + +!-------------------------------------------------------------------------------------------------- +! + tensorField_real = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + tensorField_real(1:3,1:3,i,j,k) = & + polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) + enddo; enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! doing convolution in Fourier space + call utilities_FFTtensorForward() + call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call utilities_FFTtensorBackward() + +!-------------------------------------------------------------------------------------------------- +! constructing residual + residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + +!-------------------------------------------------------------------------------------------------- +! evaluate constitutive response + P_avLastEval = P_av + call Utilities_constitutiveResponse(F_lastInc,F - residual_F_tau/polarBeta,params%timeinc, & + residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + ForwardData = .False. + +!-------------------------------------------------------------------------------------------------- +! calculate divergence + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F + call utilities_FFTtensorForward() + err_div = Utilities_divergenceRMS() + call utilities_FFTtensorBackward() + +!-------------------------------------------------------------------------------------------------- +! constructing residual + e = 0_pInt + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + e = e + 1_pInt + residual_F(1:3,1:3,i,j,k) = & + math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & + residual_F(1:3,1:3,i,j,k) - math_mul33x33(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + + residual_F_tau(1:3,1:3,i,j,k) + enddo; enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! calculating curl + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + call utilities_FFTtensorForward() + err_curl = Utilities_curlRMS() + call utilities_FFTtensorBackward() + +end subroutine Polarisation_formResidual + + +!-------------------------------------------------------------------------------------------------- +!> @brief convergence check +!-------------------------------------------------------------------------------------------------- +subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_curl_tolRel, & + err_curl_tolAbs, & + err_stress_tolAbs, & + err_stress_tolRel, & + worldrank + use math, only: & + math_mul3333xx33 + use FEsolving, only: & + terminallyIll + + implicit none + SNES :: snes_local + PetscInt :: PETScIter + PetscReal :: & + xnorm, & + snorm, & + fnorm + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode ::ierr + real(pReal) :: & + curlTol, & + divTol, & + BC_tol + +!-------------------------------------------------------------------------------------------------- +! stress BC handling + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc + err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & + mask_stress *(P_av - params%P_BC))) ! mask = 0.0 for no bc + +!-------------------------------------------------------------------------------------------------- +! error calculation + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel,err_div_tolAbs) + BC_tol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) + + converged: if ((totalIter >= itmin .and. & + all([ err_div/divTol, & + err_curl/curlTol, & + err_BC/BC_tol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then converged + reason = -1 + else converged + reason = 0 + endif converged + +!-------------------------------------------------------------------------------------------------- +! report +if (worldrank == 0_pInt) then + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' + write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' + write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) +endif + +end subroutine Polarisation_converged + +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) + use math, only: & + math_mul33x33, & + math_mul3333xx33, & + math_transpose33, & + math_rotate_backward33 + use numerics, only: & + worldrank + use mesh, only: & + grid3, & + grid + use spectral_utilities, only: & + Utilities_calculateRate, & + Utilities_forwardField, & + Utilities_updateIPcoords, & + tBoundaryCondition, & + cutBack + use IO, only: & + IO_write_JobRealFile + use FEsolving, only: & + restartWrite + + implicit none + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + P_BC, & + F_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + logical, intent(in) :: & + guess + PetscErrorCode :: ierr + PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_tau + integer(pInt) :: i, j, k + real(pReal), dimension(3,3) :: F_lambda33 + character(len=1024) :: rankStr + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) + F => xx_psc(0:8,:,:,:) + F_tau => xx_psc(9:17,:,:,:) + if (restartWrite) then + if (worldrank == 0_pInt) write(6,'(/,a)') ' writing converged results for restart' + flush(6) + write(rankStr,'(a1,i0)')'_',worldrank + call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file + write (777,rec=1) F + close (777) + call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_lastInc + close (777) + call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file + write (777,rec=1) F_tau + close (777) + call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_tau_lastInc + close (777) + if (worldrank == 0_pInt) then + call IO_write_jobRealFile(777,'F_aim',size(F_aim)) + write (777,rec=1) F_aim + close(777) + call IO_write_jobRealFile(777,'F_aim_lastInc',size(F_aim_lastInc)) + write (777,rec=1) F_aim_lastInc + close (777) + call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) + write (777,rec=1) F_aimDot + close(777) + call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg)) + write (777,rec=1) C_volAvg + close(777) + call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) + write (777,rec=1) C_volAvgLastInc + close(777) + endif + endif + call utilities_updateIPcoords(F) + + if (cutBack) then + F_aim = F_aim_lastInc + F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid3]) + F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) + C_volAvg = C_volAvgLastInc + else + ForwardData = .True. + C_volAvgLastInc = C_volAvg +!-------------------------------------------------------------------------------------------------- +! 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) + elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed + f_aimDot = F_BC%maskFloat * F_BC%values + 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 + endif + if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + F_aim_lastInc = F_aim + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + call utilities_updateIPcoords(F) + Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc_old,guess,F_lastInc, & + reshape(F,[3,3,grid(1),grid(2),grid3])) + F_tauDot = Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), & + timeinc_old,guess,F_tau_lastInc, & + reshape(F_tau,[3,3,grid(1),grid(2),grid3])) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) + F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) + endif + + F_aim = F_aim + f_aimDot * timeinc + +!-------------------------------------------------------------------------------------------------- +! update local deformation gradient + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim + math_rotate_backward33(F_aim,rotation_BC)), & + [9,grid(1),grid(2),grid3]) + F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & ! does not have any average value as boundary condition + [9,grid(1),grid(2),grid3]) + if (.not. guess) then ! large strain forwarding + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) + F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & + math_mul3333xx33(C_scale,& + math_mul33x33(math_transpose33(F_lambda33),& + F_lambda33) -math_I3))*0.5_pReal)& + + math_I3 + F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) + enddo; enddo; enddo + endif + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + +end subroutine Polarisation_forward + +!-------------------------------------------------------------------------------------------------- +!> @brief destroy routine +!-------------------------------------------------------------------------------------------------- +subroutine Polarisation_destroy() + use spectral_utilities, only: & + Utilities_destroy + + implicit none + PetscErrorCode :: ierr + + call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) + call SNESDestroy(snes,ierr); CHKERRQ(ierr) + call DMDestroy(da,ierr); CHKERRQ(ierr) + +end subroutine Polarisation_destroy + +end module spectral_mech_Polarisation diff --git a/code/spectral/spectral_thermal.f90 b/code/spectral/spectral_thermal.f90 new file mode 100644 index 000000000..843642394 --- /dev/null +++ b/code/spectral/spectral_thermal.f90 @@ -0,0 +1,419 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: spectral_thermal.f90 4082 2015-04-11 20:28:07Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Spectral solver for thermal conduction +!-------------------------------------------------------------------------------------------------- +module spectral_thermal + use prec, only: & + pInt, & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + use numerics, only: & + worldrank, & + worldsize + + implicit none + private +#include + + character (len=*), parameter, public :: & + spectral_thermal_label = 'spectralthermal' + +!-------------------------------------------------------------------------------------------------- +! derived types + type(tSolutionParams), private :: params + +!-------------------------------------------------------------------------------------------------- +! PETSc data + SNES, private :: thermal_snes + Vec, private :: solution + PetscInt, private :: xstart, xend, ystart, yend, zstart, zend + real(pReal), private, dimension(:,:,:), allocatable :: & + temperature_current, & !< field of current temperature + temperature_lastInc, & !< field of previous temperature + temperature_stagInc !< field of staggered temperature + +!-------------------------------------------------------------------------------------------------- +! reference diffusion tensor, mobility etc. + integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment + real(pReal), dimension(3,3), private :: D_ref + real(pReal), private :: mobility_ref + character(len=1024), private :: incInfo + + public :: & + spectral_thermal_init, & + spectral_thermal_solution, & + spectral_thermal_forward, & + spectral_thermal_destroy + external :: & + VecDestroy, & + DMDestroy, & + DMDACreate3D, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + PETScFinalize, & + SNESDestroy, & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber, & + SNESSolve, & + SNESSetDM, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions, & + SNESCreate, & + MPI_Abort, & + MPI_Bcast, & + MPI_Allreduce + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!-------------------------------------------------------------------------------------------------- +subroutine spectral_thermal_init + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & + IO_intOut, & + IO_read_realFile, & + IO_timeStamp + use spectral_utilities, only: & + wgt + use mesh, only: & + grid, & + grid3 + use thermal_conduction, only: & + thermal_conduction_getConductivity33, & + thermal_conduction_getMassDensity, & + thermal_conduction_getSpecificHeat + use material, only: & + mappingHomogenization, & + temperature, & + thermalMapping + + implicit none + integer(pInt), dimension(:), allocatable :: localK + integer(pInt) :: proc + integer(pInt) :: i, j, k, cell + DM :: thermal_grid + PetscScalar, pointer :: x_scal(:,:,:) + PetscErrorCode :: ierr + PetscObject :: dummy + + mainProcess: if (worldrank == 0_pInt) then + write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + +!-------------------------------------------------------------------------------------------------- +! initialize solver specific parts of PETSc + call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 + do proc = 1, worldsize + call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) + enddo + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary + DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & ! global grid + 1, 1, worldsize, & + 1, 0, & ! #dof (temperature field), ghost boundary width (domain overlap) + grid (1),grid(2),localK, & ! local grid + thermal_grid,ierr) ! handle, error + CHKERRQ(ierr) + call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da + call DMCreateGlobalVector(thermal_grid,solution ,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,spectral_thermal_formResidual,dummy,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments + +!-------------------------------------------------------------------------------------------------- +! init fields + call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,ierr) + CHKERRQ(ierr) + xend = xstart + xend - 1 + yend = ystart + yend - 1 + zend = zstart + zend - 1 + allocate(temperature_current(grid(1),grid(2),grid3), source=0.0_pReal) + allocate(temperature_lastInc(grid(1),grid(2),grid3), source=0.0_pReal) + allocate(temperature_stagInc(grid(1),grid(2),grid3), source=0.0_pReal) + cell = 0_pInt + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + temperature_current(i,j,k) = temperature(mappingHomogenization(2,1,cell))% & + p(thermalMapping(mappingHomogenization(2,1,cell))%p(1,cell)) + temperature_lastInc(i,j,k) = temperature_current(i,j,k) + temperature_stagInc(i,j,k) = temperature_current(i,j,k) + enddo; enddo; enddo + call DMDAVecGetArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with + x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current + call DMDAVecRestoreArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) + + cell = 0_pInt + D_ref = 0.0_pReal + mobility_ref = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) + mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* & + thermal_conduction_getSpecificHeat(1,cell) + enddo; enddo; enddo + D_ref = D_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + mobility_ref = mobility_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + +end subroutine spectral_thermal_init + +!-------------------------------------------------------------------------------------------------- +!> @brief solution for the Basic PETSC scheme with internal iterations +!-------------------------------------------------------------------------------------------------- +type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_old,loadCaseTime) + use numerics, only: & + itmax, & + err_thermal_tolAbs, & + err_thermal_tolRel + use spectral_utilities, only: & + tBoundaryCondition, & + Utilities_maskedCompliance, & + Utilities_updateGamma + use mesh, only: & + grid, & + grid3 + use thermal_conduction, only: & + thermal_conduction_putTemperatureAndItsRate + + implicit none + +!-------------------------------------------------------------------------------------------------- +! input data for solution + 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 + logical, intent(in) :: guess + integer(pInt) :: i, j, k, cell + PetscInt :: position + PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm + +!-------------------------------------------------------------------------------------------------- +! PETSc Data + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + spectral_thermal_solution%converged =.false. + +!-------------------------------------------------------------------------------------------------- +! set module wide availabe data + params%timeinc = timeinc + params%timeincOld = timeinc_old + + call SNESSolve(thermal_snes,PETSC_NULL_OBJECT,solution,ierr); CHKERRQ(ierr) + call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) + + if (reason < 1) then + spectral_thermal_solution%converged = .false. + spectral_thermal_solution%iterationsNeeded = itmax + else + spectral_thermal_solution%converged = .true. + spectral_thermal_solution%iterationsNeeded = totalIter + endif + stagNorm = maxval(abs(temperature_current - temperature_stagInc)) + solnNorm = maxval(abs(temperature_current)) + call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) + call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) + temperature_stagInc = temperature_current + spectral_thermal_solution%stagConverged = stagNorm < err_thermal_tolAbs & + .or. stagNorm < err_thermal_tolRel*solnNorm + +!-------------------------------------------------------------------------------------------------- +! updating thermal state + cell = 0_pInt !< material point = 0 + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt !< material point increase + call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & + (temperature_current(i,j,k)-temperature_lastInc(i,j,k))/params%timeinc, & + 1,cell) + enddo; enddo; enddo + + call VecMin(solution,position,minTemperature,ierr); CHKERRQ(ierr) + call VecMax(solution,position,maxTemperature,ierr); CHKERRQ(ierr) + if (worldrank == 0) then + if (spectral_thermal_solution%converged) & + write(6,'(/,a)') ' ... thermal conduction converged ..................................' + write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature = ',& + minTemperature, maxTemperature, stagNorm + write(6,'(/,a)') ' ===========================================================================' + flush(6) + endif + +end function spectral_thermal_solution + + +!-------------------------------------------------------------------------------------------------- +!> @brief forms the spectral thermal residual vector +!-------------------------------------------------------------------------------------------------- +subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr) + use mesh, only: & + grid, & + grid3 + use math, only: & + math_mul33x3 + use spectral_utilities, only: & + scalarField_real, & + vectorField_real, & + utilities_FFTvectorForward, & + utilities_FFTvectorBackward, & + utilities_FFTscalarForward, & + utilities_FFTscalarBackward, & + utilities_fourierGreenConvolution, & + utilities_fourierScalarGradient, & + utilities_fourierVectorDivergence + use thermal_conduction, only: & + thermal_conduction_getSourceAndItsTangent, & + thermal_conduction_getConductivity33, & + thermal_conduction_getMassDensity, & + thermal_conduction_getSpecificHeat + + implicit none + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & + in + PetscScalar, dimension( & + XG_RANGE,YG_RANGE,ZG_RANGE) :: & + x_scal + PetscScalar, dimension( & + X_RANGE,Y_RANGE,Z_RANGE) :: & + f_scal + PetscObject :: dummy + PetscErrorCode :: ierr + integer(pInt) :: i, j, k, cell + real(pReal) :: Tdot, dTdot_dT + + temperature_current = x_scal +!-------------------------------------------------------------------------------------------------- +! evaluate polarization field + scalarField_real = 0.0_pReal + scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current + call utilities_FFTscalarForward() + call utilities_fourierScalarGradient() !< calculate gradient of damage field + call utilities_FFTvectorBackward() + cell = 0_pInt + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & + vectorField_real(1:3,i,j,k)) + enddo; enddo; enddo + call utilities_FFTvectorForward() + call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field + call utilities_FFTscalarBackward() + cell = 0_pInt + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell) + scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & + params%timeinc*Tdot + & + thermal_conduction_getMassDensity (1,cell)* & + thermal_conduction_getSpecificHeat(1,cell)*(temperature_lastInc(i,j,k) - & + temperature_current(i,j,k)) + & + mobility_ref*temperature_current(i,j,k) + enddo; enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! convolution of damage field with green operator + call utilities_FFTscalarForward() + call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) + call utilities_FFTscalarBackward() + +!-------------------------------------------------------------------------------------------------- +! constructing residual + f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) + +end subroutine spectral_thermal_formResidual + +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime) + use mesh, only: & + grid, & + grid3 + use spectral_utilities, only: & + cutBack, & + wgt + use thermal_conduction, only: & + thermal_conduction_putTemperatureAndItsRate, & + thermal_conduction_getConductivity33, & + thermal_conduction_getMassDensity, & + thermal_conduction_getSpecificHeat + + implicit none + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + logical, intent(in) :: guess + integer(pInt) :: i, j, k, cell + DM :: dm_local + PetscScalar, pointer :: x_scal(:,:,:) + PetscErrorCode :: ierr + + if (cutBack) then + temperature_current = temperature_lastInc + temperature_stagInc = temperature_lastInc + +!-------------------------------------------------------------------------------------------------- +! reverting thermal field state + cell = 0_pInt !< material point = 0 + call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with + x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current + call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt !< material point increase + call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & + (temperature_current(i,j,k) - & + temperature_lastInc(i,j,k))/params%timeinc, & + 1,cell) + enddo; enddo; enddo + else +!-------------------------------------------------------------------------------------------------- +! update rate and forward last inc + temperature_lastInc = temperature_current + cell = 0_pInt + D_ref = 0.0_pReal + mobility_ref = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) + cell = cell + 1_pInt + D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) + mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* & + thermal_conduction_getSpecificHeat(1,cell) + enddo; enddo; enddo + D_ref = D_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + mobility_ref = mobility_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + endif + + end subroutine spectral_thermal_forward + +!-------------------------------------------------------------------------------------------------- +!> @brief destroy routine +!-------------------------------------------------------------------------------------------------- +subroutine spectral_thermal_destroy() + + implicit none + PetscErrorCode :: ierr + + call VecDestroy(solution,ierr); CHKERRQ(ierr) + call SNESDestroy(thermal_snes,ierr); CHKERRQ(ierr) + +end subroutine spectral_thermal_destroy + +end module spectral_thermal diff --git a/code/spectral/spectral_utilities.f90 b/code/spectral/spectral_utilities.f90 new file mode 100644 index 000000000..bde088ccb --- /dev/null +++ b/code/spectral/spectral_utilities.f90 @@ -0,0 +1,1262 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @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 +!-------------------------------------------------------------------------------------------------- +module spectral_utilities + use, intrinsic :: iso_c_binding + use prec, only: & + pReal, & + pInt + use math, only: & + math_I3 + + implicit none + private +#include + include 'fftw3-mpi.f03' + + logical, public :: cutBack =.false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill + integer(pInt), public, parameter :: maxPhaseFields = 2_pInt + integer(pInt), public :: nActiveFields = 0_pInt + +!-------------------------------------------------------------------------------------------------- +! field labels information + enum, bind(c) + enumerator :: FIELD_UNDEFINED_ID, & + FIELD_MECH_ID, & + FIELD_THERMAL_ID, & + FIELD_DAMAGE_ID, & + FIELD_VACANCYDIFFUSION_ID + end enum + +!-------------------------------------------------------------------------------------------------- +! grid related information information + real(pReal), public :: wgt !< weighting factor 1/Nelems + +!-------------------------------------------------------------------------------------------------- +! variables storing information for spectral method and FFTW + integer(pInt), public :: grid1Red !< grid(1)/2 + 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 + 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 + 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) + +!-------------------------------------------------------------------------------------------------- +! plans for FFTW + type(C_PTR), private :: & + 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) + +!-------------------------------------------------------------------------------------------------- +! variables controlling debugging + logical, private :: & + debugGeneral, & !< general debugging of spectral solver + debugRotation, & !< also printing out results in lab frame + debugPETSc !< use some in debug defined options for more verbose PETSc solution + +!-------------------------------------------------------------------------------------------------- +! derived types + type, public :: tSolutionState !< return type of solution from spectral solver variants + logical :: converged = .true. + logical :: regrid = .false. + logical :: stagConverged = .true. + logical :: termIll = .false. + integer(pInt) :: iterationsNeeded = 0_pInt + end type tSolutionState + + type, public :: tBoundaryCondition !< set of parameters defining a boundary condition + real(pReal), dimension(3,3) :: values = 0.0_pReal + real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal + logical, dimension(3,3) :: maskLogical = .false. + character(len=64) :: myType = 'None' + end type tBoundaryCondition + + 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) + real(pReal) :: time = 0.0_pReal !< length of increment + 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 + logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase + integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) + end type tLoadCase + + 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 + + type(tSolutionParams), private :: params + + type, public :: phaseFieldDataBin !< set of parameters defining a phase field + real(pReal) :: diffusion = 0.0_pReal, & !< thermal conductivity + mobility = 0.0_pReal, & !< thermal mobility + phaseField0 = 0.0_pReal !< homogeneous damage field starting condition + logical :: active = .false. + character(len=64) :: label = '' + end type phaseFieldDataBin + + enum, bind(c) + enumerator :: DERIVATIVE_CONTINUOUS_ID, & + DERIVATIVE_CENTRAL_DIFF_ID, & + DERIVATIVE_FWBW_DIFF_ID + end enum + integer(kind(DERIVATIVE_CONTINUOUS_ID)) :: & + spectral_derivative_ID + + public :: & + utilities_init, & + utilities_updateGamma, & + utilities_FFTtensorForward, & + utilities_FFTtensorBackward, & + utilities_FFTvectorForward, & + utilities_FFTvectorBackward, & + utilities_FFTscalarForward, & + utilities_FFTscalarBackward, & + utilities_fourierGammaConvolution, & + utilities_fourierGreenConvolution, & + utilities_divergenceRMS, & + utilities_curlRMS, & + utilities_fourierScalarGradient, & + utilities_fourierVectorDivergence, & + utilities_fourierVectorGradient, & + utilities_fourierTensorDivergence, & + utilities_maskedCompliance, & + utilities_constitutiveResponse, & + utilities_calculateRate, & + utilities_forwardField, & + utilities_destroy, & + utilities_updateIPcoords, & + FIELD_UNDEFINED_ID, & + FIELD_MECH_ID, & + FIELD_THERMAL_ID, & + FIELD_DAMAGE_ID + private :: & + utilities_getFreqDerivative + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields, sets debug flags, create plans for FFTW +!> @details Sets the debug levels for general, divergence, restart and FFTW from the biwise coding +!> 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. +!-------------------------------------------------------------------------------------------------- +subroutine utilities_init() + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & + IO_error, & + IO_warning, & + IO_timeStamp, & + IO_open_file + use numerics, only: & + spectral_derivative, & + fftw_planner_flag, & + fftw_timelimit, & + memory_efficient, & + petsc_defaultOptions, & + petsc_options, & + divergence_correction, & + worldrank + use debug, only: & + debug_level, & + debug_SPECTRAL, & + debug_LEVELBASIC, & + debug_SPECTRALDIVERGENCE, & + debug_SPECTRALFFTW, & + debug_SPECTRALPETSC, & + debug_SPECTRALROTATION + use debug, only: & + PETSCDEBUG + use math + use mesh, only: & + grid, & + grid3, & + grid3Offset, & + geomSize + + implicit none + + external :: & + PETScOptionsClear, & + PETScOptionsInsertString, & + MPI_Abort + + PetscErrorCode :: ierr + integer(pInt) :: i, j, k + integer(pInt), dimension(3) :: k_s + type(C_PTR) :: & + 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) + integer(C_INTPTR_T), dimension(3) :: gridFFTW + integer(C_INTPTR_T) :: alloc_local, local_K, local_K_offset + integer(C_INTPTR_T), parameter :: & + scalarSize = 1_C_INTPTR_T, & + vecSize = 3_C_INTPTR_T, & + tensorSize = 9_C_INTPTR_T + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + +!-------------------------------------------------------------------------------------------------- +! set debugging parameters + 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 + + if(debugPETSc .and. worldrank == 0_pInt) write(6,'(3(/,a),/)') & + ' Initializing PETSc with debug options: ', & + trim(PETScDebug), & + ' add more using the PETSc_Options keyword in numerics.config ' + flush(6) + call PetscOptionsClear(ierr); CHKERRQ(ierr) + if(debugPETSc) call PetscOptionsInsertString(trim(PETSCDEBUG),ierr); CHKERRQ(ierr) + call PetscOptionsInsertString(trim(petsc_defaultOptions),ierr); CHKERRQ(ierr) + call PetscOptionsInsertString(trim(petsc_options),ierr); CHKERRQ(ierr) + + grid1Red = grid(1)/2_pInt + 1_pInt + wgt = 1.0/real(product(grid),pReal) + + if (worldrank == 0) then + write(6,'(a,3(i12 ))') ' grid a b c: ', grid + write(6,'(a,3(es12.5))') ' size x y z: ', geomSize + endif + + 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 + +!-------------------------------------------------------------------------------------------------- +! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and +! resolution-independent divergence + if (divergence_correction == 1_pInt) then + do j = 1_pInt, 3_pInt + if (j /= minloc(geomSize,1) .and. j /= maxloc(geomSize,1)) & + scaledGeomSize = geomSize/geomSize(j) + enddo + elseif (divergence_correction == 2_pInt) then + do j = 1_pInt, 3_pInt + if (j /= minloc(geomSize/grid,1) .and. j /= maxloc(geomSize/grid,1)) & + scaledGeomSize = geomSize/geomSize(j)*grid(j) + enddo + else + scaledGeomSize = geomSize + endif + + +!-------------------------------------------------------------------------------------------------- +! MPI allocation + gridFFTW = int(grid,C_INTPTR_T) + alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, & + MPI_COMM_WORLD, local_K, local_K_offset) + 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 + + tensorField = fftw_alloc_complex(tensorSize*alloc_local) + call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, & + 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 + call c_f_pointer(tensorField, tensorField_fourier, [3_C_INTPTR_T,3_C_INTPTR_T, & + gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T , gridFFTW(2),local_K]) ! place a pointer for a fourier tensor representation + + vectorField = fftw_alloc_complex(vecSize*alloc_local) + call c_f_pointer(vectorField, vectorField_real, [3_C_INTPTR_T,& + 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 + call c_f_pointer(vectorField, vectorField_fourier,[3_C_INTPTR_T,& + gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T, gridFFTW(2),local_K]) ! place a pointer for a fourier vector representation + + scalarField = fftw_alloc_complex(scalarSize*alloc_local) ! allocate data for real representation (no in place transform) + call c_f_pointer(scalarField, scalarField_real, & + [2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1),gridFFTW(2),local_K]) ! place a pointer for a real scalar representation + call c_f_pointer(scalarField, scalarField_fourier, & + [ gridFFTW(1)/2_C_INTPTR_T + 1 ,gridFFTW(2),local_K]) ! place a pointer for a fourier scarlar representation + +!-------------------------------------------------------------------------------------------------- +! tensor MPI fftw plans + planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order + tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock + tensorField_real, tensorField_fourier, & ! input data, output data + MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision + if (.not. C_ASSOCIATED(planTensorForth)) call IO_error(810, ext_msg='planTensorForth') + planTensorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order + tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock + tensorField_fourier,tensorField_real, & ! input data, output data + MPI_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision + if (.not. C_ASSOCIATED(planTensorBack)) call IO_error(810, ext_msg='planTensorBack') + +!-------------------------------------------------------------------------------------------------- +! vector MPI fftw plans + planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order + vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock + vectorField_real, vectorField_fourier, & ! input data, output data + MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision + if (.not. C_ASSOCIATED(planVectorForth)) call IO_error(810, ext_msg='planVectorForth') + planVectorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order + vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock + vectorField_fourier,vectorField_real, & ! input data, output data + MPI_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision + if (.not. C_ASSOCIATED(planVectorBack)) call IO_error(810, ext_msg='planVectorBack') + +!-------------------------------------------------------------------------------------------------- +! 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 + scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock + scalarField_real, scalarField_fourier, & ! input data, output data + MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision + 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 + scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock + scalarField_fourier,scalarField_real, & ! input data, output data + MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision + if (.not. C_ASSOCIATED(planScalarBack)) call IO_error(810, ext_msg='planScalarBack') + +!-------------------------------------------------------------------------------------------------- +! 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 + + if (debugGeneral .and. worldrank == 0_pInt) write(6,'(/,a)') ' FFTW initialized' + flush(6) + +!-------------------------------------------------------------------------------------------------- +! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) + do k = grid3Offset+1_pInt, grid3Offset+grid3 + k_s(3) = k - 1_pInt + 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) + k_s(2) = j - 1_pInt + 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 + do i = 1_pInt, grid1Red + k_s(1) = i - 1_pInt ! symmetry, junst running from 0,1,...,N/2,N/2+1 + 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) + elsewhere + xi1st(1:3,i,j,k-grid3Offset) = xi2nd(1:3,i,j,k-grid3Offset) + endwhere + enddo; enddo; enddo + + if(memory_efficient) then ! allocate just single fourth order tensor + allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal)) + else ! precalculation of gamma_hat field + allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal)) + endif + +end subroutine utilities_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief updates references stiffness and potentially precalculated gamma operator +!> @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. +!> Also writes out the current reference stiffness for restart. +!-------------------------------------------------------------------------------------------------- +subroutine utilities_updateGamma(C,saveReference) + use IO, only: & + IO_write_jobRealFile + use numerics, only: & + memory_efficient, & + worldrank + use mesh, only: & + grid3Offset, & + grid3,& + grid + use math, only: & + math_det33, & + math_invert + + implicit none + 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 + complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx + real(pReal), dimension(6,6) :: matA, matInvA + integer(pInt) :: & + i, j, k, & + l, m, n, o + logical :: ierr + + C_ref = C + if (saveReference) then + 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 + endif + + if(.not. memory_efficient) then + gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A + do k = grid3Offset+1_pInt, grid3Offset+grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red + if (any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset) + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + 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) + endif + endif + enddo; enddo; enddo + endif + +end subroutine utilities_updateGamma + +!-------------------------------------------------------------------------------------------------- +!> @brief forward FFT of data in field_real to field_fourier +!> @details Does an unweighted filtered FFT transform from real to complex +!-------------------------------------------------------------------------------------------------- +subroutine utilities_FFTtensorForward() + implicit none + +!-------------------------------------------------------------------------------------------------- +! doing the tensor FFT + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) + +end subroutine utilities_FFTtensorForward + + +!-------------------------------------------------------------------------------------------------- +!> @brief backward FFT of data in field_fourier to field_real +!> @details Does an weighted inverse FFT transform from complex to real +!-------------------------------------------------------------------------------------------------- +subroutine utilities_FFTtensorBackward() + implicit none + + call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) + tensorField_real = tensorField_real * wgt ! normalize the result by number of elements + +end subroutine utilities_FFTtensorBackward + +!-------------------------------------------------------------------------------------------------- +!> @brief forward FFT of data in scalarField_real to scalarField_fourier +!> @details Does an unweighted filtered FFT transform from real to complex +!-------------------------------------------------------------------------------------------------- +subroutine utilities_FFTscalarForward() + implicit none + +!-------------------------------------------------------------------------------------------------- +! doing the scalar FFT + call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) + +end subroutine utilities_FFTscalarForward + + +!-------------------------------------------------------------------------------------------------- +!> @brief backward FFT of data in scalarField_fourier to scalarField_real +!> @details Does an weighted inverse FFT transform from complex to real +!-------------------------------------------------------------------------------------------------- +subroutine utilities_FFTscalarBackward() + 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 + +end subroutine utilities_FFTscalarBackward + + +!-------------------------------------------------------------------------------------------------- +!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed +!> @details Does an unweighted filtered FFT transform from real to complex. +!-------------------------------------------------------------------------------------------------- +subroutine utilities_FFTvectorForward() + implicit none + +!-------------------------------------------------------------------------------------------------- +! doing the vector FFT + call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) + +end subroutine utilities_FFTvectorForward + + +!-------------------------------------------------------------------------------------------------- +!> @brief backward FFT of data in field_fourier to field_real +!> @details Does an weighted inverse FFT transform from complex to real +!-------------------------------------------------------------------------------------------------- +subroutine utilities_FFTvectorBackward() + 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 + +end subroutine utilities_FFTvectorBackward + + +!-------------------------------------------------------------------------------------------------- +!> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim +!-------------------------------------------------------------------------------------------------- +subroutine utilities_fourierGammaConvolution(fieldAim) + use numerics, only: & + memory_efficient + use math, only: & + math_det33, & + math_invert + use numerics, only: & + worldrank + use mesh, only: & + grid3, & + grid, & + grid3Offset + + implicit none + real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution + complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx + real(pReal) :: matA(6,6), matInvA(6,6) + + integer(pInt) :: & + i, j, k, & + l, m, n, o + logical :: ierr + + + if (worldrank == 0_pInt) then + write(6,'(/,a)') ' ... doing gamma convolution ...............................................' + flush(6) + endif + +!-------------------------------------------------------------------------------------------------- +! do the actual spectral method calculation (mechanical equilibrium) + memoryEfficient: if(memory_efficient) then + 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 + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + 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 + gamma_hat(1:3,1:3,1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal) + endif + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + 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)) + tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex + endif + enddo; enddo; enddo + else memoryEfficient + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k)) + tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex + enddo; enddo; enddo + endif memoryEfficient + + if (grid3Offset == 0_pInt) & + tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) + +end subroutine utilities_fourierGammaConvolution + + +!-------------------------------------------------------------------------------------------------- +!> @brief doing convolution DamageGreenOp_hat * field_real +!-------------------------------------------------------------------------------------------------- +subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT) + + use math, only: & + math_mul33x3, & + PI + use mesh, only: & + grid, & + grid3 + + implicit none + 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 + complex(pReal) :: GreenOp_hat + integer(pInt) :: i, j, k + +!-------------------------------------------------------------------------------------------------- +! do the actual spectral method calculation + do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red + 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)))) ! why not use dot_product + scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat + enddo; enddo; enddo + +end subroutine utilities_fourierGreenConvolution + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate root mean square of divergence of field_fourier +!-------------------------------------------------------------------------------------------------- +real(pReal) function utilities_divergenceRMS() + use numerics, only: & + worldrank + use mesh, only: & + geomSize, & + grid, & + grid3 + + implicit none + integer(pInt) :: i, j, k + PetscErrorCode :: ierr + + external :: & + MPI_Allreduce + + if (worldrank == 0_pInt) then + write(6,'(/,a)') ' ... calculating divergence ................................................' + flush(6) + endif + +!-------------------------------------------------------------------------------------------------- +! calculating RMS divergence criterion in Fourier space + utilities_divergenceRMS = 0.0_pReal + do k = 1_pInt, grid3; do j = 1_pInt, grid(2) + do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. + utilities_divergenceRMS = utilities_divergenceRMS & + + 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)) + enddo + utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1) + + 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) + enddo; enddo + if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space + + +end function utilities_divergenceRMS + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate max of curl of field_fourier +!-------------------------------------------------------------------------------------------------- +real(pReal) function utilities_curlRMS() + use numerics, only: & + worldrank + use mesh, only: & + geomSize, & + grid, & + grid3 + + implicit none + integer(pInt) :: i, j, k, l + complex(pReal), dimension(3,3) :: curl_fourier + PetscErrorCode :: ierr + + external :: & + MPI_Reduce, & + MPI_Allreduce + + if (worldrank == 0_pInt) then + write(6,'(/,a)') ' ... calculating curl ......................................................' + flush(6) + endif + + !-------------------------------------------------------------------------------------------------- +! calculating max curl criterion in Fourier space + utilities_curlRMS = 0.0_pReal + + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); + do i = 2_pInt, grid1Red - 1_pInt + do l = 1_pInt, 3_pInt + 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)) + 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. + enddo + do l = 1_pInt, 3_pInt + 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)) + enddo + utilities_curlRMS = utilities_curlRMS + & + 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) + do l = 1_pInt, 3_pInt + 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)) + enddo + utilities_curlRMS = utilities_curlRMS + & + 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) + enddo; enddo + + call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + utilities_curlRMS = sqrt(utilities_curlRMS) * wgt + if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + +end function utilities_curlRMS + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC +!-------------------------------------------------------------------------------------------------- +function utilities_maskedCompliance(rot_BC,mask_stress,C) + use prec, only: & + prec_isNaN + use IO, only: & + IO_error + use numerics, only: & + worldrank + use math, only: & + math_Plain3333to99, & + math_plain99to3333, & + math_rotate_forward3333, & + math_rotate_forward33, & + math_invert + + implicit none + 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 + integer(pInt) :: j, k, m, n + logical, dimension(9) :: mask_stressVector + real(pReal), dimension(9,9) :: temp99_Real + integer(pInt) :: size_reduced = 0_pInt + real(pReal), dimension(:,:), allocatable :: & + s_reduced, & !< reduced compliance matrix (depending on number of stress BC) + c_reduced, & !< reduced stiffness (depending on number of stress BC) + sTimesC !< temp variable to check inversion + logical :: errmatinv + character(len=1024):: formatString + + mask_stressVector = reshape(transpose(mask_stress), [9]) + size_reduced = int(count(mask_stressVector), pInt) + 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) + allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal) + temp99_Real = math_Plain3333to99(math_rotate_forward3333(C,rot_BC)) + + if(debugGeneral .and. worldrank == 0_pInt) then + write(6,'(/,a)') ' ... updating masked compliance ............................................' + write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& + transpose(temp99_Real)/1.e9_pReal + flush(6) + endif + k = 0_pInt ! calculate reduced stiffness + 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 + c_reduced(k,j) = temp99_Real(n,m) + endif; enddo; endif; enddo + + call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness + if (any(prec_isNaN(s_reduced))) errmatinv = .true. + if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') + temp99_Real = 0.0_pReal ! fill up compliance with zeros + 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 + +!-------------------------------------------------------------------------------------------------- +! check if inversion was successful + sTimesC = matmul(c_reduced,s_reduced) + do m=1_pInt, size_reduced + do n=1_pInt, size_reduced + if(m==n .and. abs(sTimesC(m,n)) > (1.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! diagonal elements of S*C should be 1 + if(m/=n .and. abs(sTimesC(m,n)) > (0.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! off diagonal elements of S*C should be 0 + enddo + enddo + if((debugGeneral .or. errmatinv) .and. (worldrank == 0_pInt)) then ! report + write(formatString, '(I16.16)') size_reduced + formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' + 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) + endif + if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') + deallocate(c_reduced) + deallocate(s_reduced) + deallocate(sTimesC) + else + temp99_real = 0.0_pReal + endif + if(debugGeneral .and. worldrank == 0_pInt) & ! report + write(6,'(/,a,/,9(9(2x,f12.7,1x)/),/)',advance='no') ' Masked Compliance (load) * GPa =', & + transpose(temp99_Real*1.e9_pReal) + flush(6) + utilities_maskedCompliance = math_Plain99to3333(temp99_Real) + +end function utilities_maskedCompliance + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate scalar gradient in fourier field +!-------------------------------------------------------------------------------------------------- +subroutine utilities_fourierScalarGradient() + use mesh, only: & + grid3, & + grid + + implicit none + integer(pInt) :: i, j, k + + vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) & + vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) + +end subroutine utilities_fourierScalarGradient + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate vector divergence in fourier field +!-------------------------------------------------------------------------------------------------- +subroutine utilities_fourierVectorDivergence() + use mesh, only: & + grid3, & + grid + + implicit none + integer(pInt) :: i, j, k + + scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) & + scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k) + & + sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k))) + +end subroutine utilities_fourierVectorDivergence + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate vector gradient in fourier field +!-------------------------------------------------------------------------------------------------- +subroutine utilities_fourierVectorGradient() + use mesh, only: & + grid3, & + grid + + 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 + + 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 + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates constitutive response +!-------------------------------------------------------------------------------------------------- +subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & + P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC) + use debug, only: & + debug_reset, & + debug_info + use numerics, only: & + worldrank + use math, only: & + math_transpose33, & + math_rotate_forward33, & + math_det33 + use mesh, only: & + grid,& + grid3 + use FEsolving, only: & + restartWrite + use CPFEM2, only: & + CPFEM_general + use homogenization, only: & + materialpoint_F0, & + materialpoint_F, & + materialpoint_P, & + materialpoint_dPdF + + implicit none + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & + F_lastInc, & !< target deformation gradient + F !< previous deformation gradient + real(pReal), intent(in) :: timeinc !< loading time + logical, intent(in) :: forwardData !< age results + real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame + + 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 + real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress + + logical :: & + age + + integer(pInt) :: & + j,k + real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF + real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet + PetscErrorCode :: ierr + + external :: & + MPI_Reduce, & + MPI_Allreduce + + if (worldrank == 0_pInt) then + write(6,'(/,a)') ' ... evaluating constitutive response ......................................' + flush(6) + endif + age = .False. + + if (forwardData) then ! aging results + age = .True. + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) + endif + if (cutBack) then ! restore saved variables + age = .False. + endif + + materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) + call debug_reset() + +!-------------------------------------------------------------------------------------------------- +! calculate bounds of det(F) and report + if(debugGeneral) then + defgradDetMax = -huge(1.0_pReal) + defgradDetMin = +huge(1.0_pReal) + do j = 1_pInt, product(grid(1:2))*grid3 + defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j)) + defgradDetMax = max(defgradDetMax,defgradDet) + defgradDetMin = min(defgradDetMin,defgradDet) + end do + 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) + if (worldrank == 0_pInt) then + 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 + endif + + call CPFEM_general(age,timeinc) + + max_dPdF = 0.0_pReal + max_dPdF_norm = 0.0_pReal + min_dPdF = huge(1.0_pReal) + min_dPdF_norm = huge(1.0_pReal) + do k = 1_pInt, product(grid(1:2))*grid3 + 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) + endif + 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) + endif + end do + + 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) + + C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF) + C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt + + call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + + call debug_info() + + restartWrite = .false. ! reset restartWrite status + cutBack = .false. ! reset cutBack status + + P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) + P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P + call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if (debugRotation .and. worldrank == 0_pInt) & + write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& + math_transpose33(P_av)*1.e-6_pReal + P_av = math_rotate_forward33(P_av,rotation_BC) + if (worldrank == 0_pInt) then + write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& + math_transpose33(P_av)*1.e-6_pReal + flush(6) + endif + +end subroutine utilities_constitutiveResponse + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates forward rate, either guessing or just add delta/timeinc +!-------------------------------------------------------------------------------------------------- +pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,field) + use mesh, only: & + grid3, & + grid + + implicit none + real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon + real(pReal), intent(in) :: & + timeinc_old !< timeinc of last step + logical, intent(in) :: & + guess !< guess along former trajectory + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & + field_lastInc, & !< data of previous step + field !< data of current step + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & + utilities_calculateRate + + if (guess) then + utilities_calculateRate = (field-field_lastInc) / timeinc_old + else + utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3) + endif + +end function utilities_calculateRate + + +!-------------------------------------------------------------------------------------------------- +!> @brief forwards a field with a pointwise given rate, if aim is given, +!> ensures that the average matches the aim +!-------------------------------------------------------------------------------------------------- +function utilities_forwardField(timeinc,field_lastInc,rate,aim) + use mesh, only: & + grid3, & + grid + + implicit none + real(pReal), intent(in) :: & + timeinc !< timeinc of current step + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & + field_lastInc, & !< initial field + rate !< rate by which to forward + real(pReal), intent(in), optional, dimension(3,3) :: & + aim !< average field value aim + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & + utilities_forwardField + real(pReal), dimension(3,3) :: fieldDiff !< - aim + PetscErrorCode :: ierr + + external :: & + MPI_Allreduce + + utilities_forwardField = field_lastInc + rate*timeinc + if (present(aim)) then !< correct to match average + 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 + utilities_forwardField = utilities_forwardField - & + spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid3) + endif + +end function utilities_forwardField + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates filter for fourier convolution depending on type given in numerics.config +!-------------------------------------------------------------------------------------------------- +pure function utilities_getFreqDerivative(k_s) + use math, only: & + PI + use mesh, only: & + geomSize, & + grid + + implicit none + 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) + + end select + +end function utilities_getFreqDerivative + + +!-------------------------------------------------------------------------------------------------- +!> @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) + use math, only: & + math_mul33x3 + use mesh, only: & + grid, & + grid3, & + grid3Offset, & + geomSize, & + mesh_ipCoordinates + implicit none + + real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F + integer(pInt) :: i, j, k, m + real(pReal), dimension(3) :: step, offset_coords + real(pReal), dimension(3,3) :: Favg + PetscErrorCode :: ierr + external & + MPI_Bcast + +!-------------------------------------------------------------------------------------------------- +! integration in Fourier space + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + call utilities_FFTtensorForward() + call utilities_fourierTensorDivergence() + + 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) + +!-------------------------------------------------------------------------------------------------- +! average F + if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt + call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) + +!-------------------------------------------------------------------------------------------------- +! add average to fluctuation and put (0,0,0) on (0,0,0) + step = geomSize/real(grid, pReal) + if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1) + call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) + offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords + m = 1_pInt + 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) & + + offset_coords & + + math_mul33x3(Favg,step*real([i,j,k+grid3Offset]-1_pInt,pReal)) + m = m+1_pInt + enddo; enddo; enddo + +end subroutine utilities_updateIPcoords + + +!-------------------------------------------------------------------------------------------------- +!> @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 + + +end module spectral_utilities diff --git a/code/thermal/CMakeLists.txt b/code/thermal/CMakeLists.txt new file mode 100644 index 000000000..b9d3f8ff9 --- /dev/null +++ b/code/thermal/CMakeLists.txt @@ -0,0 +1,10 @@ +# group source related themal module +set (THERMAL "thermal_isothermal" + "thermal_adiabatic" + "thermal_conduction" + ) + +# compiler theraml module +foreach (p ${THERMAL}) + add_library (${p} MODULE "${p}.f90") +endforeach (p) diff --git a/code/thermal/thermal_adiabatic.f90 b/code/thermal/thermal_adiabatic.f90 new file mode 100644 index 000000000..7bb8620e7 --- /dev/null +++ b/code/thermal/thermal_adiabatic.f90 @@ -0,0 +1,422 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for adiabatic temperature evolution +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module thermal_adiabatic + use prec, only: & + pReal, & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + thermal_adiabatic_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + thermal_adiabatic_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + thermal_adiabatic_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, target, public :: & + thermal_adiabatic_Noutput !< number of outputs per instance of this thermal model + + enum, bind(c) + enumerator :: undefined_ID, & + temperature_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + thermal_adiabatic_outputID !< ID of each post result output + + + public :: & + thermal_adiabatic_init, & + thermal_adiabatic_updateState, & + thermal_adiabatic_getSourceAndItsTangent, & + thermal_adiabatic_getSpecificHeat, & + thermal_adiabatic_getMassDensity, & + thermal_adiabatic_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine thermal_adiabatic_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_warning, & + IO_error, & + IO_timeStamp, & + IO_EOF + use material, only: & + thermal_type, & + thermal_typeInstance, & + homogenization_Noutput, & + THERMAL_ADIABATIC_label, & + THERMAL_adiabatic_ID, & + material_homog, & + mappingHomogenization, & + thermalState, & + thermalMapping, & + thermal_initialT, & + temperature, & + temperatureRate, & + material_partHomogenization + use numerics,only: & + worldrank + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), allocatable, dimension(:) :: chunkPos + integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o + integer(pInt) :: sizeState + integer(pInt) :: NofMyHomog + character(len=65536) :: & + tag = '', & + line = '' + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + maxNinstance = int(count(thermal_type == THERMAL_adiabatic_ID),pInt) + if (maxNinstance == 0_pInt) return + + allocate(thermal_adiabatic_sizePostResults(maxNinstance), source=0_pInt) + allocate(thermal_adiabatic_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) + allocate(thermal_adiabatic_output (maxval(homogenization_Noutput),maxNinstance)) + thermal_adiabatic_output = '' + allocate(thermal_adiabatic_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) + allocate(thermal_adiabatic_Noutput (maxNinstance), source=0_pInt) + + rewind(fileUnit) + section = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homog part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next homog section + section = section + 1_pInt ! advance homog section counter + cycle ! skip to next line + endif + + if (section > 0_pInt ) then; if (thermal_type(section) == THERMAL_adiabatic_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + + instance = thermal_typeInstance(section) ! which instance of my thermal is present homog + chunkPos = IO_stringPos(line) + tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key + select case(tag) + case ('(output)') + select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) + case ('temperature') + thermal_adiabatic_Noutput(instance) = thermal_adiabatic_Noutput(instance) + 1_pInt + thermal_adiabatic_outputID(thermal_adiabatic_Noutput(instance),instance) = temperature_ID + thermal_adiabatic_output(thermal_adiabatic_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + end select + + end select + endif; endif + enddo parsingFile + + initializeInstances: do section = 1_pInt, size(thermal_type) + if (thermal_type(section) == THERMAL_adiabatic_ID) then + NofMyHomog=count(material_homog==section) + instance = thermal_typeInstance(section) + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,thermal_adiabatic_Noutput(instance) + select case(thermal_adiabatic_outputID(o,instance)) + case(temperature_ID) + mySize = 1_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + thermal_adiabatic_sizePostResult(o,instance) = mySize + thermal_adiabatic_sizePostResults(instance) = thermal_adiabatic_sizePostResults(instance) + mySize + endif + enddo outputsLoop + +! allocate state arrays + sizeState = 1_pInt + thermalState(section)%sizeState = sizeState + thermalState(section)%sizePostResults = thermal_adiabatic_sizePostResults(instance) + allocate(thermalState(section)%state0 (sizeState,NofMyHomog), source=thermal_initialT(section)) + allocate(thermalState(section)%subState0(sizeState,NofMyHomog), source=thermal_initialT(section)) + allocate(thermalState(section)%state (sizeState,NofMyHomog), source=thermal_initialT(section)) + + nullify(thermalMapping(section)%p) + thermalMapping(section)%p => mappingHomogenization(1,:,:) + deallocate(temperature(section)%p) + temperature(section)%p => thermalState(section)%state(1,:) + deallocate(temperatureRate(section)%p) + allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal) + + endif + + enddo initializeInstances +end subroutine thermal_adiabatic_init + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates adiabatic change in temperature based on local heat generation model +!-------------------------------------------------------------------------------------------------- +function thermal_adiabatic_updateState(subdt, ip, el) + use numerics, only: & + err_thermal_tolAbs, & + err_thermal_tolRel + use material, only: & + mappingHomogenization, & + thermalState, & + temperature, & + temperatureRate, & + thermalMapping + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + subdt + logical, dimension(2) :: & + thermal_adiabatic_updateState + integer(pInt) :: & + homog, & + offset + real(pReal) :: & + T, Tdot, dTdot_dT + + homog = mappingHomogenization(2,ip,el) + offset = mappingHomogenization(1,ip,el) + + T = thermalState(homog)%subState0(1,offset) + call thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) + T = T + subdt*Tdot/(thermal_adiabatic_getSpecificHeat(ip,el)*thermal_adiabatic_getMassDensity(ip,el)) + + thermal_adiabatic_updateState = [ abs(T - thermalState(homog)%state(1,offset)) & + <= err_thermal_tolAbs & + .or. abs(T - thermalState(homog)%state(1,offset)) & + <= err_thermal_tolRel*abs(thermalState(homog)%state(1,offset)), & + .true.] + + temperature (homog)%p(thermalMapping(homog)%p(ip,el)) = T + temperatureRate(homog)%p(thermalMapping(homog)%p(ip,el)) = & + (thermalState(homog)%state(1,offset) - thermalState(homog)%subState0(1,offset))/(subdt+tiny(0.0_pReal)) + +end function thermal_adiabatic_updateState + +!-------------------------------------------------------------------------------------------------- +!> @brief returns heat generation rate +!-------------------------------------------------------------------------------------------------- +subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) + use math, only: & + math_Mandel6to33 + use material, only: & + homogenization_Ngrains, & + mappingHomogenization, & + phaseAt, phasememberAt, & + thermal_typeInstance, & + phase_Nsources, & + phase_source, & + SOURCE_thermal_dissipation_ID, & + SOURCE_thermal_externalheat_ID + use source_thermal_dissipation, only: & + source_thermal_dissipation_getRateAndItsTangent + use source_thermal_externalheat, only: & + source_thermal_externalheat_getRateAndItsTangent + use crystallite, only: & + crystallite_Tstar_v, & + crystallite_Lp + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(out) :: & + Tdot, dTdot_dT + real(pReal) :: & + my_Tdot, my_dTdot_dT + integer(pInt) :: & + phase, & + homog, & + offset, & + instance, & + grain, & + source + + homog = mappingHomogenization(2,ip,el) + offset = mappingHomogenization(1,ip,el) + instance = thermal_typeInstance(homog) + + Tdot = 0.0_pReal + dTdot_dT = 0.0_pReal + do grain = 1, homogenization_Ngrains(homog) + phase = phaseAt(grain,ip,el) + do source = 1, phase_Nsources(phase) + select case(phase_source(source,phase)) + case (SOURCE_thermal_dissipation_ID) + call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + crystallite_Tstar_v(1:6,grain,ip,el), & + crystallite_Lp(1:3,1:3,grain,ip,el), & + grain, ip, el) + + case (SOURCE_thermal_externalheat_ID) + call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + grain, ip, el) + + case default + my_Tdot = 0.0_pReal + my_dTdot_dT = 0.0_pReal + end select + Tdot = Tdot + my_Tdot + dTdot_dT = dTdot_dT + my_dTdot_dT + enddo + enddo + + Tdot = Tdot/homogenization_Ngrains(homog) + dTdot_dT = dTdot_dT/homogenization_Ngrains(homog) + +end subroutine thermal_adiabatic_getSourceAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized specific heat capacity +!-------------------------------------------------------------------------------------------------- +function thermal_adiabatic_getSpecificHeat(ip,el) + use lattice, only: & + lattice_specificHeat + use material, only: & + homogenization_Ngrains, & + mappingHomogenization, & + material_phase + use mesh, only: & + mesh_element + use crystallite, only: & + crystallite_push33ToRef + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal) :: & + thermal_adiabatic_getSpecificHeat + integer(pInt) :: & + homog, grain + + thermal_adiabatic_getSpecificHeat = 0.0_pReal + + homog = mappingHomogenization(2,ip,el) + + do grain = 1, homogenization_Ngrains(mesh_element(3,el)) + thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat + & + lattice_specificHeat(material_phase(grain,ip,el)) + enddo + + thermal_adiabatic_getSpecificHeat = & + thermal_adiabatic_getSpecificHeat/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function thermal_adiabatic_getSpecificHeat + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized mass density +!-------------------------------------------------------------------------------------------------- +function thermal_adiabatic_getMassDensity(ip,el) + use lattice, only: & + lattice_massDensity + use material, only: & + homogenization_Ngrains, & + mappingHomogenization, & + material_phase + use mesh, only: & + mesh_element + use crystallite, only: & + crystallite_push33ToRef + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal) :: & + thermal_adiabatic_getMassDensity + integer(pInt) :: & + homog, grain + + thermal_adiabatic_getMassDensity = 0.0_pReal + + homog = mappingHomogenization(2,ip,el) + + do grain = 1, homogenization_Ngrains(mesh_element(3,el)) + thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity + & + lattice_massDensity(material_phase(grain,ip,el)) + enddo + + thermal_adiabatic_getMassDensity = & + thermal_adiabatic_getMassDensity/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function thermal_adiabatic_getMassDensity + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of thermal results +!-------------------------------------------------------------------------------------------------- +function thermal_adiabatic_postResults(ip,el) + use material, only: & + mappingHomogenization, & + thermal_typeInstance, & + thermalMapping, & + temperature + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point + el !< element + real(pReal), dimension(thermal_adiabatic_sizePostResults(thermal_typeInstance(mappingHomogenization(2,ip,el)))) :: & + thermal_adiabatic_postResults + + integer(pInt) :: & + instance, homog, offset, o, c + + homog = mappingHomogenization(2,ip,el) + offset = thermalMapping(homog)%p(ip,el) + instance = thermal_typeInstance(homog) + + c = 0_pInt + thermal_adiabatic_postResults = 0.0_pReal + + do o = 1_pInt,thermal_adiabatic_Noutput(instance) + select case(thermal_adiabatic_outputID(o,instance)) + + case (temperature_ID) + thermal_adiabatic_postResults(c+1_pInt) = temperature(homog)%p(offset) + c = c + 1 + end select + enddo +end function thermal_adiabatic_postResults + +end module thermal_adiabatic diff --git a/code/thermal/thermal_conduction.f90 b/code/thermal/thermal_conduction.f90 new file mode 100644 index 000000000..2f9b766eb --- /dev/null +++ b/code/thermal/thermal_conduction.f90 @@ -0,0 +1,444 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for temperature evolution from heat conduction +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module thermal_conduction + use prec, only: & + pReal, & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + thermal_conduction_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + thermal_conduction_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + thermal_conduction_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, target, public :: & + thermal_conduction_Noutput !< number of outputs per instance of this damage + + enum, bind(c) + enumerator :: undefined_ID, & + temperature_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + thermal_conduction_outputID !< ID of each post result output + + + public :: & + thermal_conduction_init, & + thermal_conduction_getSourceAndItsTangent, & + thermal_conduction_getConductivity33, & + thermal_conduction_getSpecificHeat, & + thermal_conduction_getMassDensity, & + thermal_conduction_putTemperatureAndItsRate, & + thermal_conduction_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine thermal_conduction_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_warning, & + IO_error, & + IO_timeStamp, & + IO_EOF + use material, only: & + thermal_type, & + thermal_typeInstance, & + homogenization_Noutput, & + THERMAL_conduction_label, & + THERMAL_conduction_ID, & + material_homog, & + mappingHomogenization, & + thermalState, & + thermalMapping, & + thermal_initialT, & + temperature, & + temperatureRate, & + material_partHomogenization + use numerics,only: & + worldrank + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), allocatable, dimension(:) :: chunkPos + integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o + integer(pInt) :: sizeState + integer(pInt) :: NofMyHomog + character(len=65536) :: & + tag = '', & + line = '' + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + maxNinstance = int(count(thermal_type == THERMAL_conduction_ID),pInt) + if (maxNinstance == 0_pInt) return + + allocate(thermal_conduction_sizePostResults(maxNinstance), source=0_pInt) + allocate(thermal_conduction_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) + allocate(thermal_conduction_output (maxval(homogenization_Noutput),maxNinstance)) + thermal_conduction_output = '' + allocate(thermal_conduction_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) + allocate(thermal_conduction_Noutput (maxNinstance), source=0_pInt) + + rewind(fileUnit) + section = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homog part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next homog section + section = section + 1_pInt ! advance homog section counter + cycle ! skip to next line + endif + + if (section > 0_pInt ) then; if (thermal_type(section) == THERMAL_conduction_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + + instance = thermal_typeInstance(section) ! which instance of my thermal is present homog + chunkPos = IO_stringPos(line) + tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key + select case(tag) + case ('(output)') + select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) + case ('temperature') + thermal_conduction_Noutput(instance) = thermal_conduction_Noutput(instance) + 1_pInt + thermal_conduction_outputID(thermal_conduction_Noutput(instance),instance) = temperature_ID + thermal_conduction_output(thermal_conduction_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + end select + + end select + endif; endif + enddo parsingFile + + initializeInstances: do section = 1_pInt, size(thermal_type) + if (thermal_type(section) == THERMAL_conduction_ID) then + NofMyHomog=count(material_homog==section) + instance = thermal_typeInstance(section) + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,thermal_conduction_Noutput(instance) + select case(thermal_conduction_outputID(o,instance)) + case(temperature_ID) + mySize = 1_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + thermal_conduction_sizePostResult(o,instance) = mySize + thermal_conduction_sizePostResults(instance) = thermal_conduction_sizePostResults(instance) + mySize + endif + enddo outputsLoop + +! allocate state arrays + sizeState = 0_pInt + thermalState(section)%sizeState = sizeState + thermalState(section)%sizePostResults = thermal_conduction_sizePostResults(instance) + allocate(thermalState(section)%state0 (sizeState,NofMyHomog)) + allocate(thermalState(section)%subState0(sizeState,NofMyHomog)) + allocate(thermalState(section)%state (sizeState,NofMyHomog)) + + nullify(thermalMapping(section)%p) + thermalMapping(section)%p => mappingHomogenization(1,:,:) + deallocate(temperature (section)%p) + allocate (temperature (section)%p(NofMyHomog), source=thermal_initialT(section)) + deallocate(temperatureRate(section)%p) + allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal) + + endif + + enddo initializeInstances +end subroutine thermal_conduction_init + +!-------------------------------------------------------------------------------------------------- +!> @brief returns heat generation rate +!-------------------------------------------------------------------------------------------------- +subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) + use math, only: & + math_Mandel6to33 + use material, only: & + homogenization_Ngrains, & + mappingHomogenization, & + phaseAt, phasememberAt, & + thermal_typeInstance, & + phase_Nsources, & + phase_source, & + SOURCE_thermal_dissipation_ID, & + SOURCE_thermal_externalheat_ID + use source_thermal_dissipation, only: & + source_thermal_dissipation_getRateAndItsTangent + use source_thermal_externalheat, only: & + source_thermal_externalheat_getRateAndItsTangent + use crystallite, only: & + crystallite_Tstar_v, & + crystallite_Lp + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(out) :: & + Tdot, dTdot_dT + real(pReal) :: & + my_Tdot, my_dTdot_dT + integer(pInt) :: & + phase, & + homog, & + offset, & + instance, & + grain, & + source + + homog = mappingHomogenization(2,ip,el) + offset = mappingHomogenization(1,ip,el) + instance = thermal_typeInstance(homog) + + Tdot = 0.0_pReal + dTdot_dT = 0.0_pReal + do grain = 1, homogenization_Ngrains(homog) + phase = phaseAt(grain,ip,el) + do source = 1, phase_Nsources(phase) + select case(phase_source(source,phase)) + case (SOURCE_thermal_dissipation_ID) + call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + crystallite_Tstar_v(1:6,grain,ip,el), & + crystallite_Lp(1:3,1:3,grain,ip,el), & + grain, ip, el) + + case (SOURCE_thermal_externalheat_ID) + call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + grain, ip, el) + + case default + my_Tdot = 0.0_pReal + my_dTdot_dT = 0.0_pReal + + end select + Tdot = Tdot + my_Tdot + dTdot_dT = dTdot_dT + my_dTdot_dT + enddo + enddo + + Tdot = Tdot/homogenization_Ngrains(homog) + dTdot_dT = dTdot_dT/homogenization_Ngrains(homog) + +end subroutine thermal_conduction_getSourceAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized thermal conductivity in reference configuration +!-------------------------------------------------------------------------------------------------- +function thermal_conduction_getConductivity33(ip,el) + use lattice, only: & + lattice_thermalConductivity33 + use material, only: & + homogenization_Ngrains, & + mappingHomogenization, & + material_phase + use mesh, only: & + mesh_element + use crystallite, only: & + crystallite_push33ToRef + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + thermal_conduction_getConductivity33 + integer(pInt) :: & + homog, & + grain + + homog = mappingHomogenization(2,ip,el) + + thermal_conduction_getConductivity33 = 0.0_pReal + do grain = 1, homogenization_Ngrains(mesh_element(3,el)) + thermal_conduction_getConductivity33 = thermal_conduction_getConductivity33 + & + crystallite_push33ToRef(grain,ip,el,lattice_thermalConductivity33(:,:,material_phase(grain,ip,el))) + enddo + + thermal_conduction_getConductivity33 = & + thermal_conduction_getConductivity33/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function thermal_conduction_getConductivity33 + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized specific heat capacity +!-------------------------------------------------------------------------------------------------- +function thermal_conduction_getSpecificHeat(ip,el) + use lattice, only: & + lattice_specificHeat + use material, only: & + homogenization_Ngrains, & + mappingHomogenization, & + material_phase + use mesh, only: & + mesh_element + use crystallite, only: & + crystallite_push33ToRef + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal) :: & + thermal_conduction_getSpecificHeat + integer(pInt) :: & + homog, grain + + thermal_conduction_getSpecificHeat = 0.0_pReal + + homog = mappingHomogenization(2,ip,el) + + do grain = 1, homogenization_Ngrains(mesh_element(3,el)) + thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat + & + lattice_specificHeat(material_phase(grain,ip,el)) + enddo + + thermal_conduction_getSpecificHeat = & + thermal_conduction_getSpecificHeat/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function thermal_conduction_getSpecificHeat + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized mass density +!-------------------------------------------------------------------------------------------------- +function thermal_conduction_getMassDensity(ip,el) + use lattice, only: & + lattice_massDensity + use material, only: & + homogenization_Ngrains, & + mappingHomogenization, & + material_phase + use mesh, only: & + mesh_element + use crystallite, only: & + crystallite_push33ToRef + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal) :: & + thermal_conduction_getMassDensity + integer(pInt) :: & + homog, grain + + thermal_conduction_getMassDensity = 0.0_pReal + + homog = mappingHomogenization(2,ip,el) + + do grain = 1, homogenization_Ngrains(mesh_element(3,el)) + thermal_conduction_getMassDensity = thermal_conduction_getMassDensity + & + lattice_massDensity(material_phase(grain,ip,el)) + enddo + + thermal_conduction_getMassDensity = & + thermal_conduction_getMassDensity/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function thermal_conduction_getMassDensity + +!-------------------------------------------------------------------------------------------------- +!> @brief updates thermal state with solution from heat conduction PDE +!-------------------------------------------------------------------------------------------------- +subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el) + use material, only: & + mappingHomogenization, & + temperature, & + temperatureRate, & + thermalMapping + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T, & + Tdot + integer(pInt) :: & + homog, & + offset + + homog = mappingHomogenization(2,ip,el) + offset = thermalMapping(homog)%p(ip,el) + temperature (homog)%p(offset) = T + temperatureRate(homog)%p(offset) = Tdot + +end subroutine thermal_conduction_putTemperatureAndItsRate + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of thermal results +!-------------------------------------------------------------------------------------------------- +function thermal_conduction_postResults(ip,el) + use material, only: & + mappingHomogenization, & + thermal_typeInstance, & + temperature, & + thermalMapping + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point + el !< element + real(pReal), dimension(thermal_conduction_sizePostResults(thermal_typeInstance(mappingHomogenization(2,ip,el)))) :: & + thermal_conduction_postResults + + integer(pInt) :: & + instance, homog, offset, o, c + + homog = mappingHomogenization(2,ip,el) + offset = thermalMapping(homog)%p(ip,el) + instance = thermal_typeInstance(homog) + + c = 0_pInt + thermal_conduction_postResults = 0.0_pReal + + do o = 1_pInt,thermal_conduction_Noutput(instance) + select case(thermal_conduction_outputID(o,instance)) + + case (temperature_ID) + thermal_conduction_postResults(c+1_pInt) = temperature(homog)%p(offset) + c = c + 1 + end select + enddo +end function thermal_conduction_postResults + +end module thermal_conduction diff --git a/code/thermal/thermal_isothermal.f90 b/code/thermal/thermal_isothermal.f90 new file mode 100644 index 000000000..8c9d3a782 --- /dev/null +++ b/code/thermal/thermal_isothermal.f90 @@ -0,0 +1,65 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for isothermal temperature field +!-------------------------------------------------------------------------------------------------- +module thermal_isothermal + + implicit none + private + + public :: & + thermal_isothermal_init + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields, reads information from material configuration file +!-------------------------------------------------------------------------------------------------- +subroutine thermal_isothermal_init() + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use prec, only: & + pReal, & + pInt + use IO, only: & + IO_timeStamp + use material + use numerics, only: & + worldrank + + implicit none + integer(pInt) :: & + homog, & + NofMyHomog, & + sizeState + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + initializeInstances: do homog = 1_pInt, material_Nhomogenization + + myhomog: if (thermal_type(homog) == THERMAL_isothermal_ID) then + NofMyHomog = count(material_homog == homog) + sizeState = 0_pInt + thermalState(homog)%sizeState = sizeState + thermalState(homog)%sizePostResults = sizeState + allocate(thermalState(homog)%state0 (sizeState,NofMyHomog), source=0.0_pReal) + allocate(thermalState(homog)%subState0(sizeState,NofMyHomog), source=0.0_pReal) + allocate(thermalState(homog)%state (sizeState,NofMyHomog), source=0.0_pReal) + + deallocate(temperature (homog)%p) + allocate (temperature (homog)%p(1), source=thermal_initialT(homog)) + deallocate(temperatureRate(homog)%p) + allocate (temperatureRate(homog)%p(1), source=0.0_pReal) + + endif myhomog + enddo initializeInstances + + +end subroutine thermal_isothermal_init + +end module thermal_isothermal diff --git a/code/vacancyflux/CMakeLists.txt b/code/vacancyflux/CMakeLists.txt new file mode 100644 index 000000000..4a4c4f31c --- /dev/null +++ b/code/vacancyflux/CMakeLists.txt @@ -0,0 +1,10 @@ +# group source file +set (VACANCYFLUX "vacancyflux_isoconc" + "vacancyflux_isochempot" + "vacancyflux_cahnhilliard" + ) + +# compiler as module +foreach (p ${VACANCYFLUX}) + add_library (${p} MODULE "${p}.f90") +endforeach (p) diff --git a/code/vacancyflux/vacancyflux_cahnhilliard.f90 b/code/vacancyflux/vacancyflux_cahnhilliard.f90 new file mode 100644 index 000000000..16a380ffc --- /dev/null +++ b/code/vacancyflux/vacancyflux_cahnhilliard.f90 @@ -0,0 +1,606 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for conservative transport of vacancy concentration field +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module vacancyflux_cahnhilliard + use prec, only: & + pReal, & + pInt, & + p_vec + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + vacancyflux_cahnhilliard_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + vacancyflux_cahnhilliard_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + vacancyflux_cahnhilliard_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, target, public :: & + vacancyflux_cahnhilliard_Noutput !< number of outputs per instance of this damage + + real(pReal), dimension(:), allocatable, private :: & + vacancyflux_cahnhilliard_flucAmplitude + + type(p_vec), dimension(:), allocatable, private :: & + vacancyflux_cahnhilliard_thermalFluc + + real(pReal), parameter, private :: & + kB = 1.3806488e-23_pReal !< Boltzmann constant in J/Kelvin + + enum, bind(c) + enumerator :: undefined_ID, & + vacancyConc_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + vacancyflux_cahnhilliard_outputID !< ID of each post result output + + + public :: & + vacancyflux_cahnhilliard_init, & + vacancyflux_cahnhilliard_getSourceAndItsTangent, & + vacancyflux_cahnhilliard_getMobility33, & + vacancyflux_cahnhilliard_getDiffusion33, & + vacancyflux_cahnhilliard_getChemPotAndItsTangent, & + vacancyflux_cahnhilliard_putVacancyConcAndItsRate, & + vacancyflux_cahnhilliard_postResults + private :: & + vacancyflux_cahnhilliard_getFormationEnergy, & + vacancyflux_cahnhilliard_getEntropicCoeff, & + vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine vacancyflux_cahnhilliard_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_warning, & + IO_error, & + IO_timeStamp, & + IO_EOF + use material, only: & + vacancyflux_type, & + vacancyflux_typeInstance, & + homogenization_Noutput, & + VACANCYFLUX_cahnhilliard_label, & + VACANCYFLUX_cahnhilliard_ID, & + material_homog, & + mappingHomogenization, & + vacancyfluxState, & + vacancyfluxMapping, & + vacancyConc, & + vacancyConcRate, & + vacancyflux_initialCv, & + material_partHomogenization, & + material_partPhase + use numerics,only: & + worldrank + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), allocatable, dimension(:) :: chunkPos + integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o,offset + integer(pInt) :: sizeState + integer(pInt) :: NofMyHomog + character(len=65536) :: & + tag = '', & + line = '' + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_cahnhilliard_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_cahnhilliard_ID),pInt) + if (maxNinstance == 0_pInt) return + + allocate(vacancyflux_cahnhilliard_sizePostResults(maxNinstance), source=0_pInt) + allocate(vacancyflux_cahnhilliard_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) + allocate(vacancyflux_cahnhilliard_output (maxval(homogenization_Noutput),maxNinstance)) + vacancyflux_cahnhilliard_output = '' + allocate(vacancyflux_cahnhilliard_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) + allocate(vacancyflux_cahnhilliard_Noutput (maxNinstance), source=0_pInt) + + allocate(vacancyflux_cahnhilliard_flucAmplitude (maxNinstance)) + allocate(vacancyflux_cahnhilliard_thermalFluc (maxNinstance)) + + rewind(fileUnit) + section = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to + line = IO_read(fileUnit) + enddo + + parsingHomog: do while (trim(line) /= IO_EOF) ! read through sections of homog part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next homog section + section = section + 1_pInt ! advance homog section counter + cycle ! skip to next line + endif + + if (section > 0_pInt ) then; if (vacancyflux_type(section) == VACANCYFLUX_cahnhilliard_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + + instance = vacancyflux_typeInstance(section) ! which instance of my vacancyflux is present homog + chunkPos = IO_stringPos(line) + tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key + select case(tag) + case ('(output)') + select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) + case ('vacancyconc') + vacancyflux_cahnhilliard_Noutput(instance) = vacancyflux_cahnhilliard_Noutput(instance) + 1_pInt + vacancyflux_cahnhilliard_outputID(vacancyflux_cahnhilliard_Noutput(instance),instance) = vacancyConc_ID + vacancyflux_cahnhilliard_output(vacancyflux_cahnhilliard_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + end select + + case ('vacancyflux_flucamplitude') + vacancyflux_cahnhilliard_flucAmplitude(instance) = IO_floatValue(line,chunkPos,2_pInt) + + end select + endif; endif + enddo parsingHomog + + initializeInstances: do section = 1_pInt, size(vacancyflux_type) + if (vacancyflux_type(section) == VACANCYFLUX_cahnhilliard_ID) then + NofMyHomog=count(material_homog==section) + instance = vacancyflux_typeInstance(section) + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,vacancyflux_cahnhilliard_Noutput(instance) + select case(vacancyflux_cahnhilliard_outputID(o,instance)) + case(vacancyConc_ID) + mySize = 1_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + vacancyflux_cahnhilliard_sizePostResult(o,instance) = mySize + vacancyflux_cahnhilliard_sizePostResults(instance) = vacancyflux_cahnhilliard_sizePostResults(instance) + mySize + endif + enddo outputsLoop + +! allocate state arrays + sizeState = 0_pInt + vacancyfluxState(section)%sizeState = sizeState + vacancyfluxState(section)%sizePostResults = vacancyflux_cahnhilliard_sizePostResults(instance) + allocate(vacancyfluxState(section)%state0 (sizeState,NofMyHomog)) + allocate(vacancyfluxState(section)%subState0(sizeState,NofMyHomog)) + allocate(vacancyfluxState(section)%state (sizeState,NofMyHomog)) + + allocate(vacancyflux_cahnhilliard_thermalFluc(instance)%p(NofMyHomog)) + do offset = 1_pInt, NofMyHomog + call random_number(vacancyflux_cahnhilliard_thermalFluc(instance)%p(offset)) + vacancyflux_cahnhilliard_thermalFluc(instance)%p(offset) = & + 1.0_pReal - & + vacancyflux_cahnhilliard_flucAmplitude(instance)* & + (vacancyflux_cahnhilliard_thermalFluc(instance)%p(offset) - 0.5_pReal) + enddo + + nullify(vacancyfluxMapping(section)%p) + vacancyfluxMapping(section)%p => mappingHomogenization(1,:,:) + deallocate(vacancyConc (section)%p) + allocate (vacancyConc (section)%p(NofMyHomog), source=vacancyflux_initialCv(section)) + deallocate(vacancyConcRate(section)%p) + allocate (vacancyConcRate(section)%p(NofMyHomog), source=0.0_pReal) + + endif + + enddo initializeInstances + +end subroutine vacancyflux_cahnhilliard_init + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates homogenized vacancy driving forces +!-------------------------------------------------------------------------------------------------- +subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv, ip, el) + use material, only: & + homogenization_Ngrains, & + mappingHomogenization, & + phaseAt, phasememberAt, & + phase_source, & + phase_Nsources, & + SOURCE_vacancy_phenoplasticity_ID, & + SOURCE_vacancy_irradiation_ID, & + SOURCE_vacancy_thermalfluc_ID + use source_vacancy_phenoplasticity, only: & + source_vacancy_phenoplasticity_getRateAndItsTangent + use source_vacancy_irradiation, only: & + source_vacancy_irradiation_getRateAndItsTangent + use source_vacancy_thermalfluc, only: & + source_vacancy_thermalfluc_getRateAndItsTangent + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + Cv + integer(pInt) :: & + phase, & + grain, & + source + real(pReal) :: & + CvDot, dCvDot_dCv, localCvDot, dLocalCvDot_dCv + + CvDot = 0.0_pReal + dCvDot_dCv = 0.0_pReal + do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el)) + phase = phaseAt(grain,ip,el) + do source = 1_pInt, phase_Nsources(phase) + select case(phase_source(source,phase)) + case (SOURCE_vacancy_phenoplasticity_ID) + call source_vacancy_phenoplasticity_getRateAndItsTangent (localCvDot, dLocalCvDot_dCv, grain, ip, el) + + case (SOURCE_vacancy_irradiation_ID) + call source_vacancy_irradiation_getRateAndItsTangent (localCvDot, dLocalCvDot_dCv, grain, ip, el) + + case (SOURCE_vacancy_thermalfluc_ID) + call source_vacancy_thermalfluc_getRateAndItsTangent(localCvDot, dLocalCvDot_dCv, grain, ip, el) + + end select + CvDot = CvDot + localCvDot + dCvDot_dCv = dCvDot_dCv + dLocalCvDot_dCv + enddo + enddo + + CvDot = CvDot/homogenization_Ngrains(mappingHomogenization(2,ip,el)) + dCvDot_dCv = dCvDot_dCv/homogenization_Ngrains(mappingHomogenization(2,ip,el)) + +end subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized vacancy mobility tensor in reference configuration +!-------------------------------------------------------------------------------------------------- +function vacancyflux_cahnhilliard_getMobility33(ip,el) + use lattice, only: & + lattice_vacancyfluxMobility33 + use material, only: & + homogenization_Ngrains, & + material_phase + use mesh, only: & + mesh_element + use crystallite, only: & + crystallite_push33ToRef + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + vacancyflux_cahnhilliard_getMobility33 + integer(pInt) :: & + grain + + vacancyflux_cahnhilliard_getMobility33 = 0.0_pReal + do grain = 1, homogenization_Ngrains(mesh_element(3,el)) + vacancyflux_cahnhilliard_getMobility33 = vacancyflux_cahnhilliard_getMobility33 + & + crystallite_push33ToRef(grain,ip,el,lattice_vacancyfluxMobility33(:,:,material_phase(grain,ip,el))) + enddo + + vacancyflux_cahnhilliard_getMobility33 = & + vacancyflux_cahnhilliard_getMobility33/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function vacancyflux_cahnhilliard_getMobility33 + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized vacancy diffusion tensor in reference configuration +!-------------------------------------------------------------------------------------------------- +function vacancyflux_cahnhilliard_getDiffusion33(ip,el) + use lattice, only: & + lattice_vacancyfluxDiffusion33 + use material, only: & + homogenization_Ngrains, & + material_phase + use mesh, only: & + mesh_element + use crystallite, only: & + crystallite_push33ToRef + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + vacancyflux_cahnhilliard_getDiffusion33 + integer(pInt) :: & + grain + + vacancyflux_cahnhilliard_getDiffusion33 = 0.0_pReal + do grain = 1, homogenization_Ngrains(mesh_element(3,el)) + vacancyflux_cahnhilliard_getDiffusion33 = vacancyflux_cahnhilliard_getDiffusion33 + & + crystallite_push33ToRef(grain,ip,el,lattice_vacancyfluxDiffusion33(:,:,material_phase(grain,ip,el))) + enddo + + vacancyflux_cahnhilliard_getDiffusion33 = & + vacancyflux_cahnhilliard_getDiffusion33/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function vacancyflux_cahnhilliard_getDiffusion33 + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized vacancy formation energy +!-------------------------------------------------------------------------------------------------- +real(pReal) function vacancyflux_cahnhilliard_getFormationEnergy(ip,el) + use lattice, only: & + lattice_vacancyFormationEnergy, & + lattice_vacancyVol, & + lattice_vacancySurfaceEnergy + use material, only: & + homogenization_Ngrains, & + material_phase + use mesh, only: & + mesh_element + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + grain + + vacancyflux_cahnhilliard_getFormationEnergy = 0.0_pReal + do grain = 1, homogenization_Ngrains(mesh_element(3,el)) + vacancyflux_cahnhilliard_getFormationEnergy = vacancyflux_cahnhilliard_getFormationEnergy + & + lattice_vacancyFormationEnergy(material_phase(grain,ip,el))/ & + lattice_vacancyVol(material_phase(grain,ip,el))/ & + lattice_vacancySurfaceEnergy(material_phase(grain,ip,el)) + enddo + + vacancyflux_cahnhilliard_getFormationEnergy = & + vacancyflux_cahnhilliard_getFormationEnergy/ & + homogenization_Ngrains(mesh_element(3,el)) + +end function vacancyflux_cahnhilliard_getFormationEnergy + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized vacancy entropy coefficient +!-------------------------------------------------------------------------------------------------- +real(pReal) function vacancyflux_cahnhilliard_getEntropicCoeff(ip,el) + use lattice, only: & + lattice_vacancyVol, & + lattice_vacancySurfaceEnergy + use material, only: & + homogenization_Ngrains, & + material_homog, & + material_phase, & + temperature, & + thermalMapping + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + grain + + vacancyflux_cahnhilliard_getEntropicCoeff = 0.0_pReal + do grain = 1, homogenization_Ngrains(material_homog(ip,el)) + vacancyflux_cahnhilliard_getEntropicCoeff = vacancyflux_cahnhilliard_getEntropicCoeff + & + kB/ & + lattice_vacancyVol(material_phase(grain,ip,el))/ & + lattice_vacancySurfaceEnergy(material_phase(grain,ip,el)) + enddo + + vacancyflux_cahnhilliard_getEntropicCoeff = & + vacancyflux_cahnhilliard_getEntropicCoeff* & + temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))/ & + homogenization_Ngrains(material_homog(ip,el)) + +end function vacancyflux_cahnhilliard_getEntropicCoeff + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized kinematic contribution to chemical potential +!-------------------------------------------------------------------------------------------------- +subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCv, Cv, ip, el) + use lattice, only: & + lattice_vacancySurfaceEnergy + use material, only: & + homogenization_Ngrains, & + material_homog, & + phase_kinematics, & + phase_Nkinematics, & + material_phase, & + KINEMATICS_vacancy_strain_ID + use crystallite, only: & + crystallite_Tstar_v, & + crystallite_Fi0, & + crystallite_Fi + use kinematics_vacancy_strain, only: & + kinematics_vacancy_strain_ChemPotAndItsTangent + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + Cv + real(pReal), intent(out) :: & + KPot, dKPot_dCv + real(pReal) :: & + my_KPot, my_dKPot_dCv + integer(pInt) :: & + grain, kinematics + + KPot = 0.0_pReal + dKPot_dCv = 0.0_pReal + do grain = 1_pInt,homogenization_Ngrains(material_homog(ip,el)) + do kinematics = 1_pInt, phase_Nkinematics(material_phase(grain,ip,el)) + select case (phase_kinematics(kinematics,material_phase(grain,ip,el))) + case (KINEMATICS_vacancy_strain_ID) + call kinematics_vacancy_strain_ChemPotAndItsTangent(my_KPot, my_dKPot_dCv, & + crystallite_Tstar_v(1:6,grain,ip,el), & + crystallite_Fi0(1:3,1:3,grain,ip,el), & + crystallite_Fi (1:3,1:3,grain,ip,el), & + grain,ip, el) + + case default + my_KPot = 0.0_pReal + my_dKPot_dCv = 0.0_pReal + + end select + KPot = KPot + my_KPot/lattice_vacancySurfaceEnergy(material_phase(grain,ip,el)) + dKPot_dCv = dKPot_dCv + my_dKPot_dCv/lattice_vacancySurfaceEnergy(material_phase(grain,ip,el)) + enddo + enddo + + KPot = KPot/homogenization_Ngrains(material_homog(ip,el)) + dKPot_dCv = dKPot_dCv/homogenization_Ngrains(material_homog(ip,el)) + +end subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized chemical potential and its tangent +!-------------------------------------------------------------------------------------------------- +subroutine vacancyflux_cahnhilliard_getChemPotAndItsTangent(ChemPot,dChemPot_dCv,Cv,ip,el) + use numerics, only: & + vacancyBoundPenalty, & + vacancyPolyOrder + use material, only: & + mappingHomogenization, & + vacancyflux_typeInstance, & + porosity, & + porosityMapping + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + Cv + real(pReal), intent(out) :: & + ChemPot, & + dChemPot_dCv + real(pReal) :: & + VoidPhaseFrac, kBT, KPot, dKPot_dCv + integer(pInt) :: & + homog, o + + homog = mappingHomogenization(2,ip,el) + VoidPhaseFrac = porosity(homog)%p(porosityMapping(homog)%p(ip,el)) + kBT = vacancyflux_cahnhilliard_getEntropicCoeff(ip,el) + + ChemPot = vacancyflux_cahnhilliard_getFormationEnergy(ip,el) + dChemPot_dCv = 0.0_pReal + do o = 1_pInt, vacancyPolyOrder + ChemPot = ChemPot + kBT*((2.0_pReal*Cv - 1.0_pReal)**real(2_pInt*o-1_pInt,pReal))/ & + real(2_pInt*o-1_pInt,pReal) + dChemPot_dCv = dChemPot_dCv + 2.0_pReal*kBT*(2.0_pReal*Cv - 1.0_pReal)**real(2_pInt*o-2_pInt,pReal) + enddo + + ChemPot = VoidPhaseFrac*VoidPhaseFrac*ChemPot & + - 2.0_pReal*(1.0_pReal - Cv)*(1.0_pReal - VoidPhaseFrac)*(1.0_pReal - VoidPhaseFrac) + + dChemPot_dCv = VoidPhaseFrac*VoidPhaseFrac*dChemPot_dCv & + + 2.0_pReal*(1.0_pReal - VoidPhaseFrac)*(1.0_pReal - VoidPhaseFrac) + + call vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCv, Cv, ip, el) + ChemPot = ChemPot + KPot + dChemPot_dCv = dChemPot_dCv + dKPot_dCv + + if (Cv < 0.0_pReal) then + ChemPot = ChemPot - 3.0_pReal*vacancyBoundPenalty*Cv*Cv + dChemPot_dCv = dChemPot_dCv - 6.0_pReal*vacancyBoundPenalty*Cv + elseif (Cv > 1.0_pReal) then + ChemPot = ChemPot + 3.0_pReal*vacancyBoundPenalty*(1.0_pReal - Cv)*(1.0_pReal - Cv) + dChemPot_dCv = dChemPot_dCv - 6.0_pReal*vacancyBoundPenalty*(1.0_pReal - Cv) + endif + + ChemPot = ChemPot* & + vacancyflux_cahnhilliard_thermalFluc(vacancyflux_typeInstance(homog))%p(mappingHomogenization(1,ip,el)) + dChemPot_dCv = dChemPot_dCv* & + vacancyflux_cahnhilliard_thermalFluc(vacancyflux_typeInstance(homog))%p(mappingHomogenization(1,ip,el)) + +end subroutine vacancyflux_cahnhilliard_getChemPotAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief updated vacancy concentration and its rate with solution from transport PDE +!-------------------------------------------------------------------------------------------------- +subroutine vacancyflux_cahnhilliard_putVacancyConcAndItsRate(Cv,Cvdot,ip,el) + use material, only: & + mappingHomogenization, & + vacancyConc, & + vacancyConcRate, & + vacancyfluxMapping + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + Cv, & + Cvdot + integer(pInt) :: & + homog, & + offset + + homog = mappingHomogenization(2,ip,el) + offset = vacancyfluxMapping(homog)%p(ip,el) + vacancyConc (homog)%p(offset) = Cv + vacancyConcRate(homog)%p(offset) = Cvdot + +end subroutine vacancyflux_cahnhilliard_putVacancyConcAndItsRate + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of vacancy transport results +!-------------------------------------------------------------------------------------------------- +function vacancyflux_cahnhilliard_postResults(ip,el) + use material, only: & + mappingHomogenization, & + vacancyflux_typeInstance, & + vacancyConc, & + vacancyfluxMapping + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point + el !< element + real(pReal), dimension(vacancyflux_cahnhilliard_sizePostResults(vacancyflux_typeInstance(mappingHomogenization(2,ip,el)))) :: & + vacancyflux_cahnhilliard_postResults + + integer(pInt) :: & + instance, homog, offset, o, c + + homog = mappingHomogenization(2,ip,el) + offset = vacancyfluxMapping(homog)%p(ip,el) + instance = vacancyflux_typeInstance(homog) + + c = 0_pInt + vacancyflux_cahnhilliard_postResults = 0.0_pReal + + do o = 1_pInt,vacancyflux_cahnhilliard_Noutput(instance) + select case(vacancyflux_cahnhilliard_outputID(o,instance)) + + case (vacancyConc_ID) + vacancyflux_cahnhilliard_postResults(c+1_pInt) = vacancyConc(homog)%p(offset) + c = c + 1 + end select + enddo +end function vacancyflux_cahnhilliard_postResults + +end module vacancyflux_cahnhilliard diff --git a/code/vacancyflux/vacancyflux_isochempot.f90 b/code/vacancyflux/vacancyflux_isochempot.f90 new file mode 100644 index 000000000..35db8d159 --- /dev/null +++ b/code/vacancyflux/vacancyflux_isochempot.f90 @@ -0,0 +1,329 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for locally evolving vacancy concentration +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module vacancyflux_isochempot + use prec, only: & + pReal, & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + vacancyflux_isochempot_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + vacancyflux_isochempot_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + vacancyflux_isochempot_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, target, public :: & + vacancyflux_isochempot_Noutput !< number of outputs per instance of this damage + + enum, bind(c) + enumerator :: undefined_ID, & + vacancyconc_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + vacancyflux_isochempot_outputID !< ID of each post result output + + + public :: & + vacancyflux_isochempot_init, & + vacancyflux_isochempot_updateState, & + vacancyflux_isochempot_getSourceAndItsTangent, & + vacancyflux_isochempot_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine vacancyflux_isochempot_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_warning, & + IO_error, & + IO_timeStamp, & + IO_EOF + use material, only: & + vacancyflux_type, & + vacancyflux_typeInstance, & + homogenization_Noutput, & + VACANCYFLUX_isochempot_label, & + VACANCYFLUX_isochempot_ID, & + material_homog, & + mappingHomogenization, & + vacancyfluxState, & + vacancyfluxMapping, & + vacancyConc, & + vacancyConcRate, & + vacancyflux_initialCv, & + material_partHomogenization + use numerics,only: & + worldrank + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), allocatable, dimension(:) :: chunkPos + integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o + integer(pInt) :: sizeState + integer(pInt) :: NofMyHomog + character(len=65536) :: & + tag = '', & + line = '' + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isochempot_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_isochempot_ID),pInt) + if (maxNinstance == 0_pInt) return + + allocate(vacancyflux_isochempot_sizePostResults(maxNinstance), source=0_pInt) + allocate(vacancyflux_isochempot_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) + allocate(vacancyflux_isochempot_output (maxval(homogenization_Noutput),maxNinstance)) + vacancyflux_isochempot_output = '' + allocate(vacancyflux_isochempot_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) + allocate(vacancyflux_isochempot_Noutput (maxNinstance), source=0_pInt) + + rewind(fileUnit) + section = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homog part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next homog section + section = section + 1_pInt ! advance homog section counter + cycle ! skip to next line + endif + + if (section > 0_pInt ) then; if (vacancyflux_type(section) == VACANCYFLUX_isochempot_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + + instance = vacancyflux_typeInstance(section) ! which instance of my vacancyflux is present homog + chunkPos = IO_stringPos(line) + tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key + select case(tag) + case ('(output)') + select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) + case ('vacancyconc') + vacancyflux_isochempot_Noutput(instance) = vacancyflux_isochempot_Noutput(instance) + 1_pInt + vacancyflux_isochempot_outputID(vacancyflux_isochempot_Noutput(instance),instance) = vacancyconc_ID + vacancyflux_isochempot_output(vacancyflux_isochempot_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + end select + + end select + endif; endif + enddo parsingFile + + initializeInstances: do section = 1_pInt, size(vacancyflux_type) + if (vacancyflux_type(section) == VACANCYFLUX_isochempot_ID) then + NofMyHomog=count(material_homog==section) + instance = vacancyflux_typeInstance(section) + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,vacancyflux_isochempot_Noutput(instance) + select case(vacancyflux_isochempot_outputID(o,instance)) + case(vacancyconc_ID) + mySize = 1_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + vacancyflux_isochempot_sizePostResult(o,instance) = mySize + vacancyflux_isochempot_sizePostResults(instance) = vacancyflux_isochempot_sizePostResults(instance) + mySize + endif + enddo outputsLoop + +! allocate state arrays + sizeState = 1_pInt + vacancyfluxState(section)%sizeState = sizeState + vacancyfluxState(section)%sizePostResults = vacancyflux_isochempot_sizePostResults(instance) + allocate(vacancyfluxState(section)%state0 (sizeState,NofMyHomog), source=vacancyflux_initialCv(section)) + allocate(vacancyfluxState(section)%subState0(sizeState,NofMyHomog), source=vacancyflux_initialCv(section)) + allocate(vacancyfluxState(section)%state (sizeState,NofMyHomog), source=vacancyflux_initialCv(section)) + + nullify(vacancyfluxMapping(section)%p) + vacancyfluxMapping(section)%p => mappingHomogenization(1,:,:) + deallocate(vacancyConc(section)%p) + vacancyConc(section)%p => vacancyfluxState(section)%state(1,:) + deallocate(vacancyConcRate(section)%p) + allocate(vacancyConcRate(section)%p(NofMyHomog), source=0.0_pReal) + + endif + + enddo initializeInstances +end subroutine vacancyflux_isochempot_init + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates change in vacancy concentration based on local vacancy generation model +!-------------------------------------------------------------------------------------------------- +function vacancyflux_isochempot_updateState(subdt, ip, el) + use numerics, only: & + err_vacancyflux_tolAbs, & + err_vacancyflux_tolRel + use material, only: & + mappingHomogenization, & + vacancyflux_typeInstance, & + vacancyfluxState, & + vacancyConc, & + vacancyConcRate, & + vacancyfluxMapping + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + subdt + logical, dimension(2) :: & + vacancyflux_isochempot_updateState + integer(pInt) :: & + homog, & + offset, & + instance + real(pReal) :: & + Cv, Cvdot, dCvDot_dCv + + homog = mappingHomogenization(2,ip,el) + offset = mappingHomogenization(1,ip,el) + instance = vacancyflux_typeInstance(homog) + + Cv = vacancyfluxState(homog)%subState0(1,offset) + call vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv, ip, el) + Cv = Cv + subdt*Cvdot + + vacancyflux_isochempot_updateState = [ abs(Cv - vacancyfluxState(homog)%state(1,offset)) & + <= err_vacancyflux_tolAbs & + .or. abs(Cv - vacancyfluxState(homog)%state(1,offset)) & + <= err_vacancyflux_tolRel*abs(vacancyfluxState(homog)%state(1,offset)), & + .true.] + + vacancyConc (homog)%p(vacancyfluxMapping(homog)%p(ip,el)) = Cv + vacancyConcRate(homog)%p(vacancyfluxMapping(homog)%p(ip,el)) = & + (vacancyfluxState(homog)%state(1,offset) - vacancyfluxState(homog)%subState0(1,offset))/(subdt+tiny(0.0_pReal)) + +end function vacancyflux_isochempot_updateState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates homogenized vacancy driving forces +!-------------------------------------------------------------------------------------------------- +subroutine vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv, ip, el) + use material, only: & + homogenization_Ngrains, & + mappingHomogenization, & + phaseAt, phasememberAt, & + phase_source, & + phase_Nsources, & + SOURCE_vacancy_phenoplasticity_ID, & + SOURCE_vacancy_irradiation_ID, & + SOURCE_vacancy_thermalfluc_ID + use source_vacancy_phenoplasticity, only: & + source_vacancy_phenoplasticity_getRateAndItsTangent + use source_vacancy_irradiation, only: & + source_vacancy_irradiation_getRateAndItsTangent + use source_vacancy_thermalfluc, only: & + source_vacancy_thermalfluc_getRateAndItsTangent + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + Cv + integer(pInt) :: & + phase, & + grain, & + source + real(pReal) :: & + CvDot, dCvDot_dCv, localCvDot, dLocalCvDot_dCv + + CvDot = 0.0_pReal + dCvDot_dCv = 0.0_pReal + do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el)) + phase = phaseAt(grain,ip,el) + do source = 1_pInt, phase_Nsources(phase) + select case(phase_source(source,phase)) + case (SOURCE_vacancy_phenoplasticity_ID) + call source_vacancy_phenoplasticity_getRateAndItsTangent (localCvDot, dLocalCvDot_dCv, grain, ip, el) + + case (SOURCE_vacancy_irradiation_ID) + call source_vacancy_irradiation_getRateAndItsTangent (localCvDot, dLocalCvDot_dCv, grain, ip, el) + + case (SOURCE_vacancy_thermalfluc_ID) + call source_vacancy_thermalfluc_getRateAndItsTangent(localCvDot, dLocalCvDot_dCv, grain, ip, el) + + end select + CvDot = CvDot + localCvDot + dCvDot_dCv = dCvDot_dCv + dLocalCvDot_dCv + enddo + enddo + + CvDot = CvDot/homogenization_Ngrains(mappingHomogenization(2,ip,el)) + dCvDot_dCv = dCvDot_dCv/homogenization_Ngrains(mappingHomogenization(2,ip,el)) + +end subroutine vacancyflux_isochempot_getSourceAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of vacancy transport results +!-------------------------------------------------------------------------------------------------- +function vacancyflux_isochempot_postResults(ip,el) + use material, only: & + mappingHomogenization, & + vacancyflux_typeInstance, & + vacancyConc, & + vacancyfluxMapping + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point + el !< element + real(pReal), dimension(vacancyflux_isochempot_sizePostResults(vacancyflux_typeInstance(mappingHomogenization(2,ip,el)))) :: & + vacancyflux_isochempot_postResults + + integer(pInt) :: & + instance, homog, offset, o, c + + homog = mappingHomogenization(2,ip,el) + offset = vacancyfluxMapping(homog)%p(ip,el) + instance = vacancyflux_typeInstance(homog) + + c = 0_pInt + vacancyflux_isochempot_postResults = 0.0_pReal + + do o = 1_pInt,vacancyflux_isochempot_Noutput(instance) + select case(vacancyflux_isochempot_outputID(o,instance)) + + case (vacancyconc_ID) + vacancyflux_isochempot_postResults(c+1_pInt) = vacancyConc(homog)%p(offset) + c = c + 1 + end select + enddo +end function vacancyflux_isochempot_postResults + +end module vacancyflux_isochempot diff --git a/code/vacancyflux/vacancyflux_isoconc.f90 b/code/vacancyflux/vacancyflux_isoconc.f90 new file mode 100644 index 000000000..63cfb1b62 --- /dev/null +++ b/code/vacancyflux/vacancyflux_isoconc.f90 @@ -0,0 +1,63 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for constant vacancy concentration +!-------------------------------------------------------------------------------------------------- +module vacancyflux_isoconc + + implicit none + private + + public :: & + vacancyflux_isoconc_init + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields, reads information from material configuration file +!-------------------------------------------------------------------------------------------------- +subroutine vacancyflux_isoconc_init() + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use prec, only: & + pReal, & + pInt + use IO, only: & + IO_timeStamp + use material + use numerics, only: & + worldrank + + implicit none + integer(pInt) :: & + homog, & + NofMyHomog + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isoconc_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + initializeInstances: do homog = 1_pInt, material_Nhomogenization + + myhomog: if (vacancyflux_type(homog) == VACANCYFLUX_isoconc_ID) then + NofMyHomog = count(material_homog == homog) + vacancyfluxState(homog)%sizeState = 0_pInt + vacancyfluxState(homog)%sizePostResults = 0_pInt + allocate(vacancyfluxState(homog)%state0 (0_pInt,NofMyHomog)) + allocate(vacancyfluxState(homog)%subState0(0_pInt,NofMyHomog)) + allocate(vacancyfluxState(homog)%state (0_pInt,NofMyHomog)) + + deallocate(vacancyConc (homog)%p) + allocate (vacancyConc (homog)%p(1), source=vacancyflux_initialCv(homog)) + deallocate(vacancyConcRate(homog)%p) + allocate (vacancyConcRate(homog)%p(1), source=0.0_pReal) + + endif myhomog + enddo initializeInstances + + +end subroutine vacancyflux_isoconc_init + +end module vacancyflux_isoconc diff --git a/lib/damask/core.so b/lib/damask/core.so new file mode 100755 index 000000000..7f744dafa Binary files /dev/null and b/lib/damask/core.so differ diff --git a/lib/damask/corientation.so b/lib/damask/corientation.so new file mode 100755 index 000000000..76fcf94f9 Binary files /dev/null and b/lib/damask/corientation.so differ