set up build system with CMake

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}
+

code/DAMASK_marc2011.f90

@@ -0,0 +1 @@
+DAMASK_marc.f90

code/DAMASK_marc2012.f90

@@ -0,0 +1 @@
+DAMASK_marc.f90

code/DAMASK_marc2013.1.f90

@@ -0,0 +1 @@
+DAMASK_marc.f90

code/DAMASK_marc2013.f90

@@ -0,0 +1 @@
+DAMASK_marc.f90

code/DAMASK_marc2014.2.f90

@@ -0,0 +1 @@
+DAMASK_marc.f90

code/DAMASK_marc2014.f90

@@ -0,0 +1 @@
+DAMASK_marc.f90

code/DAMASK_marc2015.f90

@@ -0,0 +1 @@
+DAMASK_marc.f90

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
+LIBZ    = /mnt/research/CMM/opt/hdf5/lib/libz.a
-LIBSZ   = $(HDF5)/lib/libsz.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
+HDFLIBS = -I$(HDF5)/include -L$(HDF5)/lib
-HDFLIBZ = $(LIBZ) $(LIBSZ) -lm
+HDFLIBZ = -L$(LIBZ) -L$(LIBSZ)
 
 # MPI compiler wrappers will tell if they are pointing to ifort or gfortran
 COMPILEROUT :=$(shell $(FC) -show)
@@ -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
+damask_hdf5.o:                 damask_hdf5.f90 \
-	$(PREFIX) $(COMPILERNAME) $(COMPILE) -c damask_hdf5.f90 $(HDFLIBS) $(HDFLIBZ)
+                               prec.o \
+                               IO.o
+	$(PREFIX) $(COMPILERNAME) $(HDFLIBS) $(HDFLIBZ) -c damask_hdf5.f90 $(SUFFIX) -lm
+endif
 
 %.o : %.f90
 	$(PREFIX) $(COMPILERNAME) $(COMPILE) -c $< $(SUFFIX)

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)"
+

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

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)

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 <petsc/finclude/petsc.h90>
+
+ 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

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 <petsc/finclude/petscsys.h>
+#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

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 <petsc/finclude/petsc.h90>
+
+ 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

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 <petsc/finclude/petsc.h90>
+
+ 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

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 <petsc/finclude/petsc.h90>
+
+ 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

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 <petsc/finclude/petsc.h90>
+
+ 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

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)

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 <homogenization>
+   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

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 <homogenization>
+   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

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

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)

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 <homogenization>
+   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

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 <homogenization>
+   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

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