diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index f1e67a115..17728632d 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -99,6 +99,7 @@ checkout: - git submodule update --init - source env/DAMASK.sh - make processing + - mkdir /tmp/$CI_PIPELINE_ID except: - master - release @@ -108,7 +109,7 @@ Pytest_python: stage: python script: - cd $DAMASKROOT/python - - pytest + - pytest --basetemp=/tmp/${CI_PIPELINE_ID}/python except: - master - release @@ -194,7 +195,7 @@ grid_mech_compile_Intel: - cp -r grid_mech_compile grid_mech_compile_Intel - grid_mech_compile_Intel/test.py - cd pytest - - pytest -k 'compile and grid' + - pytest -k 'compile and grid' --basetemp=/tmp/${CI_PIPELINE_ID}/compile_grid_Intel except: - master - release @@ -206,7 +207,7 @@ Compile_FEM_Intel: - cp -r FEM_compile FEM_compile_Intel - FEM_compile_Intel/test.py - cd pytest - - pytest -k 'compile and mesh' + - pytest -k 'compile and mesh' --basetemp=/tmp/${CI_PIPELINE_ID}/compile_mesh_Intel except: - master - release @@ -218,7 +219,7 @@ grid_mech_compile_GNU: - cp -r grid_mech_compile grid_mech_compile_GNU - grid_mech_compile_GNU/test.py - cd pytest - - pytest -k 'compile and grid' + - pytest -k 'compile and grid' --basetemp=/tmp/${CI_PIPELINE_ID}/compile_grid_GNU except: - master - release @@ -230,7 +231,7 @@ Compile_FEM_GNU: - cp -r FEM_compile FEM_compile_GNU - FEM_compile_GNU/test.py - cd pytest - - pytest -k 'compile and mesh' + - pytest -k 'compile and mesh' --basetemp=/tmp/${CI_PIPELINE_ID}/compile_mesh_GNU except: - master - release @@ -252,7 +253,7 @@ Pytest_grid: script: - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel - cd pytest - - pytest -m 'not compile' + - pytest -m 'not compile' --basetemp=/tmp/${CI_PIPELINE_ID}/fortran except: - master - release diff --git a/PRIVATE b/PRIVATE index dde6a5836..1e8c66897 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit dde6a5836089b14bd39d4393a47b82b351d2bd85 +Subproject commit 1e8c66897820468ab46958d995005e2b69204d0e diff --git a/VERSION b/VERSION index c7934861f..2e3222238 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha-523-gf4b5d5f24 +v3.0.0-alpha-556-gf379aecd8 diff --git a/examples/SpectralMethod/Polycrystal/shearXY.load b/examples/SpectralMethod/Polycrystal/shearXY.load deleted file mode 100644 index 9f89ef886..000000000 --- a/examples/SpectralMethod/Polycrystal/shearXY.load +++ /dev/null @@ -1 +0,0 @@ -fdot 0 0 0 1e-3 0 0 0 0 0 stress * * * * * * * * * time 60 incs 120 freq 20 diff --git a/examples/SpectralMethod/Polycrystal/shearXY.yaml b/examples/SpectralMethod/Polycrystal/shearXY.yaml new file mode 100644 index 000000000..39966d181 --- /dev/null +++ b/examples/SpectralMethod/Polycrystal/shearXY.yaml @@ -0,0 +1,9 @@ +step: + - mech: + dot_F: [0, 0, 0, + 1e-3, 0, 0, + 0, 0, 0] + discretization: + t: 60 + N: 120 + f_out: 20 diff --git a/examples/SpectralMethod/Polycrystal/shearZX.load b/examples/SpectralMethod/Polycrystal/shearZX.load deleted file mode 100644 index 6a3606340..000000000 --- a/examples/SpectralMethod/Polycrystal/shearZX.load +++ /dev/null @@ -1 +0,0 @@ -fdot 0 0 1.0e-3 0 0 0 0 0 0 stress * * * * * * * * * time 60 incs 120 freq 20 diff --git a/examples/SpectralMethod/Polycrystal/shearZX.yaml b/examples/SpectralMethod/Polycrystal/shearZX.yaml new file mode 100644 index 000000000..4395ecf17 --- /dev/null +++ b/examples/SpectralMethod/Polycrystal/shearZX.yaml @@ -0,0 +1,10 @@ +--- +step: + - mech: + dot_F: [0, 0, 1e-3, + 0, 0, 0, + 0, 0, 0] + discretization: + t: 60 + N: 120 + f_out: 20 diff --git a/examples/SpectralMethod/Polycrystal/tensionX.load b/examples/SpectralMethod/Polycrystal/tensionX.load deleted file mode 100644 index b0af80ea8..000000000 --- a/examples/SpectralMethod/Polycrystal/tensionX.load +++ /dev/null @@ -1,2 +0,0 @@ -fdot 1.0e-3 0 0 0 * 0 0 0 * stress * * * * 0 * * * 0 time 10 incs 40 freq 4 -fdot 1.0e-3 0 0 0 * 0 0 0 * stress * * * * 0 * * * 0 time 60 incs 60 diff --git a/examples/SpectralMethod/Polycrystal/tensionX.yaml b/examples/SpectralMethod/Polycrystal/tensionX.yaml new file mode 100644 index 000000000..6e86fdcf2 --- /dev/null +++ b/examples/SpectralMethod/Polycrystal/tensionX.yaml @@ -0,0 +1,25 @@ +--- + +step: + - mech: + dot_F: [1.0e-3, 0, 0, + 0, x, 0, + 0, 0, x] + P: [x, x, x, + x, 0, x, + x, x, 0] + discretization: + t: 10 + N: 40 + f_out: 4 + - mech: + dot_F: [1.0e-3, 0, 0, + 0, x, 0, + 0, 0, x] + P: [x, x, x, + x, 0, x, + x, x, 0] + discretization: + t: 60 + N: 60 + f_out: 4 diff --git a/src/IO.f90 b/src/IO.f90 index a95938a7f..260bd94b8 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -491,7 +491,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (705) msg = 'Unsupported feature' case (706) - msg = 'Access by incorrect node type' + msg = 'Type mismatch in YAML data node' case (707) msg = 'Abrupt end of file' case (708) diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index 09a2d0592..f11faf54b 100644 --- a/src/YAML_types.f90 +++ b/src/YAML_types.f90 @@ -311,7 +311,7 @@ function tNode_asScalar(self) result(scalar) class is(tScalar) scalar => self class default - call IO_error(706,ext_msg='tNode_asScalar') + call IO_error(706,ext_msg='Expected "scalar"') end select end function tNode_asScalar @@ -329,7 +329,7 @@ function tNode_asList(self) result(list) class is(tList) list => self class default - call IO_error(706,ext_msg='tNode_asList') + call IO_error(706,ext_msg='Expected "list"') end select end function tNode_asList @@ -347,7 +347,7 @@ function tNode_asDict(self) result(dict) class is(tDict) dict => self class default - call IO_error(706,ext_msg='tNode_asDict') + call IO_error(706,ext_msg='Expected "dict"') end select end function tNode_asDict @@ -641,7 +641,7 @@ function tNode_contains(self,k) result(exists) endif enddo else - call IO_error(706,ext_msg='tNode_contains') + call IO_error(706,ext_msg='Expected "list" or "dict"') endif end function tNode_contains diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index a1bb23375..4b8302def 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -24,20 +24,13 @@ program DAMASK_grid use grid_damage_spectral use grid_thermal_spectral use results - + implicit none !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0) logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors - integer, allocatable, dimension(:) :: chunkPos - integer :: & - N_t = 0, & !< # of time indicators found in load case file - N_n = 0, & !< # of increment specifiers found in load case file - N_def = 0 !< # of rate of deformation specifiers found in load case file - character(len=:), allocatable :: & - line !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. @@ -57,26 +50,23 @@ program DAMASK_grid stagIterate, & cutBack = .false. integer :: & - i, j, k, l, field, & + i, j, m, field, & errorID = 0, & cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ - stepFraction = 0 !< fraction of current time interval - integer :: & - currentLoadcase = 0, & !< current load case + stepFraction = 0, & !< fraction of current time interval + l = 0, & !< current load case inc, & !< current increment in current load case totalIncsCounter = 0, & !< total # of increments statUnit = 0, & !< file unit for statistics output stagIter, & - nActiveFields = 0 - character(len=pStringLen), dimension(:), allocatable :: fileContent + nActiveFields = 0, & + maxCutBack, & !< max number of cut backs + stagItMax !< max number of field level staggered iterations character(len=pStringLen) :: & incInfo, & loadcase_string - integer :: & - maxCutBack, & !< max number of cut backs - stagItMax !< max number of field level staggered iterations + type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases - type(tLoadCase) :: newLoadCase type(tSolutionState), allocatable, dimension(:) :: solres procedure(grid_mech_spectral_basic_init), pointer :: & mech_init @@ -88,19 +78,26 @@ program DAMASK_grid mech_updateCoords procedure(grid_mech_spectral_basic_restartWrite), pointer :: & mech_restartWrite - + external :: & quit class (tNode), pointer :: & num_grid, & - debug_grid ! pointer to grid debug options + debug_grid, & ! pointer to grid debug options + config_load, & + load_steps, & + load_step, & + step_mech, & + step_discretization, & + step_deformation, & + step_stress !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) - - call CPFEM_initAll + + call CPFEM_initAll print'(/,a)', ' <<<+- DAMASK_spectral init -+>>>'; flush(IO_STDOUT) - + print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' @@ -110,7 +107,6 @@ program DAMASK_grid if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1 if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 allocate(solres(nActiveFields)) - allocate(newLoadCase%ID(nActiveFields)) !------------------------------------------------------------------------------------------------- ! reading field paramters from numerics file and do sanity checks @@ -123,16 +119,16 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! assign mechanics solver depending on selected type - + debug_grid => config_debug%get('grid',defaultVal=emptyList) - select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) + select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) case ('Basic') mech_init => grid_mech_spectral_basic_init mech_forward => grid_mech_spectral_basic_forward mech_solution => grid_mech_spectral_basic_solution mech_updateCoords => grid_mech_spectral_basic_updateCoords mech_restartWrite => grid_mech_spectral_basic_restartWrite - + case ('Polarisation') if(debug_grid%contains('basic')) & call IO_warning(42, ext_msg='debug Divergence') @@ -141,7 +137,7 @@ program DAMASK_grid mech_solution => grid_mech_spectral_polarisation_solution mech_updateCoords => grid_mech_spectral_polarisation_updateCoords mech_restartWrite => grid_mech_spectral_polarisation_restartWrite - + case ('FEM') if(debug_grid%contains('basic')) & call IO_warning(42, ext_msg='debug Divergence') @@ -150,169 +146,128 @@ program DAMASK_grid mech_solution => grid_mech_FEM_solution mech_updateCoords => grid_mech_FEM_updateCoords mech_restartWrite => grid_mech_FEM_restartWrite - + case default call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver'))) - + end select - + !-------------------------------------------------------------------------------------------------- -! reading information from load case file and to sanity checks - fileContent = IO_readlines(trim(interface_loadFile)) - if(size(fileContent) == 0) call IO_error(307,ext_msg='No load case specified') - - allocate (loadCases(0)) ! array of load cases - do currentLoadCase = 1, size(fileContent) - line = fileContent(currentLoadCase) - if (IO_isBlank(line)) cycle - chunkPos = IO_stringPos(line) - - do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase - select case (IO_lc(IO_stringValue(line,chunkPos,i))) - case('l','fdot','dotf','f') - N_def = N_def + 1 - case('t','time','delta') - N_t = N_t + 1 - case('n','incs','increments','logincs','logincrements') - N_n = N_n + 1 - end select - enddo - if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check - call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase - - newLoadCase%stress%myType='p' +! reading information from load case file and to sanity checks + config_load => YAML_parse_file(trim(interface_loadFile)) + + load_steps => config_load%get('step') + allocate(loadCases(load_steps%length)) ! array of load cases + + do l = 1, load_steps%length + + allocate(loadCases(l)%ID(nActiveFields)) field = 1 - newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default + loadCases(l)%ID(field) = FIELD_MECH_ID ! mechanical active by default thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then field = field + 1 - newLoadCase%ID(field) = FIELD_THERMAL_ID + loadCases(l)%ID(field) = FIELD_THERMAL_ID endif thermalActive damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then field = field + 1 - newLoadCase%ID(field) = FIELD_DAMAGE_ID + loadCases(l)%ID(field) = FIELD_DAMAGE_ID endif damageActive - - call newLoadCase%rot%fromEulers(real([0.0,0.0,0.0],pReal)) - readIn: do i = 1, chunkPos(1) - select case (IO_lc(IO_stringValue(line,chunkPos,i))) - case('fdot','dotf','l','f') ! assign values for the deformation BC matrix + + load_step => load_steps%get(l) + + step_mech => load_step%get('mech') + loadCases(l)%stress%myType='P' + readMech: do m = 1, step_mech%length + select case (step_mech%getKey(m)) + case('dot_F','L','F') ! assign values for the deformation BC matrix + loadCases(l)%deformation%myType = step_mech%getKey(m) temp_valueVector = 0.0_pReal - if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'fdot'.or. & ! in case of Fdot, set type to fdot - IO_lc(IO_stringValue(line,chunkPos,i)) == 'dotf') then - newLoadCase%deformation%myType = 'fdot' - else if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'f') then - newLoadCase%deformation%myType = 'f' - else - newLoadCase%deformation%myType = 'l' - endif + + step_deformation => step_mech%get(m) do j = 1, 9 - temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a * - if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable + temp_maskVector(j) = step_deformation%get_asString(j) /= 'x' ! true if not a 'x' + if (temp_maskVector(j)) temp_valueVector(j) = step_deformation%get_asFloat(j) ! read value where applicable enddo - newLoadCase%deformation%mask = transpose(reshape(temp_maskVector,[ 3,3])) ! mask in 3x3 notation - newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation - case('p','stress', 's') + loadCases(l)%deformation%mask = transpose(reshape(temp_maskVector,[ 3,3])) ! mask in 3x3 notation + loadCases(l)%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation + case('P') temp_valueVector = 0.0_pReal + step_stress => step_mech%get(m) do j = 1, 9 - temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk - if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable + temp_maskVector(j) = step_stress%get_asString(j) /= 'x' ! true if not a 'x' + if (temp_maskVector(j)) temp_valueVector(j) = step_stress%get_asFloat(j) ! read value where applicable enddo - newLoadCase%stress%mask = transpose(reshape(temp_maskVector,[ 3,3])) - newLoadCase%stress%values = math_9to33(temp_valueVector) - case('t','time','delta') ! increment time - newLoadCase%time = IO_floatValue(line,chunkPos,i+1) - case('n','incs','increments') ! number of increments - newLoadCase%incs = IO_intValue(line,chunkPos,i+1) - case('logincs','logincrements') ! number of increments (switch to log time scaling) - newLoadCase%incs = IO_intValue(line,chunkPos,i+1) - newLoadCase%logscale = 1 - case('freq','frequency','outputfreq') ! frequency of result writings - newLoadCase%outputfrequency = IO_intValue(line,chunkPos,i+1) - case('r','restart','restartwrite') ! frequency of writing restart information - newLoadCase%restartfrequency = IO_intValue(line,chunkPos,i+1) - case('guessreset','dropguessing') - newLoadCase%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory - case('euler') ! rotation of load case given in euler angles - temp_valueVector = 0.0_pReal - l = 1 ! assuming values given in degrees - k = 1 ! assuming keyword indicating degree/radians present - select case (IO_lc(IO_stringValue(line,chunkPos,i+1))) - case('deg','degree') - case('rad','radian') ! don't convert from degree to radian - l = 0 - case default - k = 0 - end select - do j = 1, 3 - temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j) - enddo - call newLoadCase%rot%fromEulers(temp_valueVector(1:3),degrees=(l==1)) - case('rotation','rot') ! assign values for the rotation matrix - temp_valueVector = 0.0_pReal - do j = 1, 9 - temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) - enddo - call newLoadCase%rot%fromMatrix(math_9to33(temp_valueVector)) + loadCases(l)%stress%mask = transpose(reshape(temp_maskVector,[ 3,3])) + loadCases(l)%stress%values = math_9to33(temp_valueVector) end select - enddo readIn - - newLoadCase%followFormerTrajectory = merge(.true.,.false.,currentLoadCase > 1) ! by default, guess from previous load case - + call loadCases(l)%rot%fromAxisAngle(step_mech%get_asFloats('R', & + defaultVal = real([0.0,0.0,1.0,0.0],pReal)),degrees=.true.) + enddo readMech + if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/F/dot_F missing') + + step_discretization => load_step%get('discretization') + if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing') + if(.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing') + loadCases(l)%time = step_discretization%get_asFloat('t') + loadCases(l)%incs = step_discretization%get_asFloat('N') + loadCases(l)%logscale = step_discretization%get_asBool ('log_timestep', defaultVal= .false.) + loadCases(l)%outputfrequency = step_discretization%get_asInt ('f_out', defaultVal=1) + loadCases(l)%restartfrequency = step_discretization%get_asInt ('f_restart', defaultVal=huge(0)) + + loadCases(l)%followFormerTrajectory = .not. (load_step%get_asBool('drop_guessing',defaultVal=.false.) .or. & + merge(.false.,.true.,l > 1)) ! do not continue to predict deformation along former trajectory + reportAndCheck: if (worldrank == 0) then - write (loadcase_string, '(i0)' ) currentLoadCase - print'(/,a,i0)', ' load case: ', currentLoadCase - if (.not. newLoadCase%followFormerTrajectory) & + write (loadcase_string, '(i0)' ) l + print'(/,a,i0)', ' load case: ', l + if (.not. loadCases(l)%followFormerTrajectory) & print*, ' drop guessing along trajectory' - if (newLoadCase%deformation%myType == 'l') then + if (loadCases(l)%deformation%myType == 'L') then do j = 1, 3 - if (any(newLoadCase%deformation%mask(j,1:3) .eqv. .true.) .and. & - any(newLoadCase%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined + if (any(loadCases(l)%deformation%mask(j,1:3) .eqv. .true.) .and. & + any(loadCases(l)%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined enddo print*, ' velocity gradient:' - else if (newLoadCase%deformation%myType == 'f') then + else if (loadCases(l)%deformation%myType == 'F') then print*, ' deformation gradient at end of load case:' - else + else if (loadCases(l)%deformation%myType == 'dot_F') then print*, ' deformation gradient rate:' endif do i = 1, 3; do j = 1, 3 - if(newLoadCase%deformation%mask(i,j)) then - write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j) + if(loadCases(l)%deformation%mask(i,j)) then + write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j) else - write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' + write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' endif enddo; write(IO_STDOUT,'(/)',advance='no') enddo - if (any(newLoadCase%stress%mask .eqv. newLoadCase%deformation%mask)) errorID = 831 ! exclusive or masking only - if (any(newLoadCase%stress%mask .and. transpose(newLoadCase%stress%mask) .and. (math_I3<1))) & + if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831 ! exclusive or masking only + if (any(loadCases(l)%stress%mask .and. transpose(loadCases(l)%stress%mask) .and. (math_I3<1))) & errorID = 838 ! no rotation is allowed by stress BC print*, ' stress / GPa:' do i = 1, 3; do j = 1, 3 - if(newLoadCase%stress%mask(i,j)) then - write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal + if(loadCases(l)%stress%mask(i,j)) then + write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%stress%values(i,j)*1e-9_pReal else - write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' + write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' endif enddo; write(IO_STDOUT,'(/)',advance='no') enddo - if (any(abs(matmul(newLoadCase%rot%asMatrix(), & - transpose(newLoadCase%rot%asMatrix()))-math_I3) > & - reshape(spread(tol_math_check,1,9),[ 3,3]))) errorID = 846 ! given rotation matrix contains strain - if (any(dNeq(newLoadCase%rot%asMatrix(), math_I3))) & + if (any(dNeq(loadCases(l)%rot%asMatrix(), math_I3))) & write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& - transpose(newLoadCase%rot%asMatrix()) - if (newLoadCase%time < 0.0_pReal) errorID = 834 ! negative time increment - print'(a,f0.3)', ' time: ', newLoadCase%time - if (newLoadCase%incs < 1) errorID = 835 ! non-positive incs count - print'(a,i0)', ' increments: ', newLoadCase%incs - if (newLoadCase%outputfrequency < 1) errorID = 836 ! non-positive result frequency - print'(a,i0)', ' output frequency: ', newLoadCase%outputfrequency - if (newLoadCase%restartfrequency < 1) errorID = 839 ! non-positive restart frequency - if (newLoadCase%restartfrequency < huge(0)) & - print'(a,i0)', ' restart frequency: ', newLoadCase%restartfrequency + transpose(loadCases(l)%rot%asMatrix()) + if (loadCases(l)%time < 0.0_pReal) errorID = 834 ! negative time increment + print'(a,f0.3)', ' time: ', loadCases(l)%time + if (loadCases(l)%incs < 1) errorID = 835 ! non-positive incs count + print'(a,i0)', ' increments: ', loadCases(l)%incs + if (loadCases(l)%outputfrequency < 1) errorID = 836 ! non-positive result frequency + print'(a,i0)', ' output frequency: ', loadCases(l)%outputfrequency + if (loadCases(l)%restartfrequency < 1) errorID = 839 ! non-positive restart frequency + if (loadCases(l)%restartfrequency < huge(0)) & + print'(a,i0)', ' restart frequency: ', loadCases(l)%restartfrequency if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message endif reportAndCheck - loadCases = [loadCases,newLoadCase] ! load case is ok, append it enddo !-------------------------------------------------------------------------------------------------- @@ -322,13 +277,13 @@ program DAMASK_grid select case (loadCases(1)%ID(field)) case(FIELD_MECH_ID) call mech_init - + case(FIELD_THERMAL_ID) call grid_thermal_spectral_init - + case(FIELD_DAMAGE_ID) call grid_damage_spectral_init - + end select enddo @@ -345,45 +300,45 @@ program DAMASK_grid '.sta',form='FORMATTED', position='APPEND', status='OLD') endif writeHeader endif - + writeUndeformed: if (interface_restartInc < 1) then print'(/,a)', ' ... writing initial configuration to file ........................' call CPFEM_results(0,0.0_pReal) endif writeUndeformed - - loadCaseLooping: do currentLoadCase = 1, size(loadCases) + + loadCaseLooping: do l = 1, size(loadCases) time0 = time ! load case start time - guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc - - incLooping: do inc = 1, loadCases(currentLoadCase)%incs + guess = loadCases(l)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc + + incLooping: do inc = 1, loadCases(l)%incs totalIncsCounter = totalIncsCounter + 1 !-------------------------------------------------------------------------------------------------- ! forwarding time timeIncOld = timeinc ! last timeinc that brought former inc to an end - if (loadCases(currentLoadCase)%logscale == 0) then ! linear scale - timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) + if (.not. loadCases(l)%logscale) then ! linear scale + timeinc = loadCases(l)%time/real(loadCases(l)%incs,pReal) else - if (currentLoadCase == 1) then ! 1st load case of logarithmic scale + if (l == 1) then ! 1st load case of logarithmic scale timeinc = loadCases(1)%time*(2.0_pReal**real(max(inc-1,1)-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd else ! not-1st load case of logarithmic scale timeinc = time0 * & - ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal))& - -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1 ,pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal))) + ( (1.0_pReal + loadCases(l)%time/time0 )**(real( inc ,pReal)/& + real(loadCases(l)%incs ,pReal))& + -(1.0_pReal + loadCases(l)%time/time0 )**(real( inc-1 ,pReal)/& + real(loadCases(l)%incs ,pReal))) endif endif timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step - + skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc? time = time + timeinc ! just advance time, skip already performed calculation guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference else skipping stepFraction = 0 ! fraction scaled by stepFactor**cutLevel - + subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) - remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time + remainingLoadCaseTime = loadCases(l)%time+time0 - time time = time + timeinc ! forward target time stepFraction = stepFraction + 1 ! count step @@ -392,9 +347,9 @@ program DAMASK_grid print'(/,a)', ' ###########################################################################' print'(1x,a,es12.5,6(a,i0))', & 'Time', time, & - 's: Increment ', inc,'/',loadCases(currentLoadCase)%incs,& + 's: Increment ', inc,'/',loadCases(l)%incs,& '-', stepFraction,'/',subStepFactor**cutBackLevel,& - ' of load case ', currentLoadCase,'/',size(loadCases) + ' of load case ', l,'/',size(loadCases) write(incInfo,'(4(a,i0))') & 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& '-', stepFraction,'/',subStepFactor**cutBackLevel @@ -403,14 +358,14 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! forward fields do field = 1, nActiveFields - select case(loadCases(currentLoadCase)%ID(field)) + select case(loadCases(l)%ID(field)) case(FIELD_MECH_ID) call mech_forward (& cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & - deformation_BC = loadCases(currentLoadCase)%deformation, & - stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rot) - + deformation_BC = loadCases(l)%deformation, & + stress_BC = loadCases(l)%stress, & + rotation_BC = loadCases(l)%rot) + case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) end select @@ -423,7 +378,7 @@ program DAMASK_grid stagIterate = .true. do while (stagIterate) do field = 1, nActiveFields - select case(loadCases(currentLoadCase)%ID(field)) + select case(loadCases(l)%ID(field)) case(FIELD_MECH_ID) solres(field) = mech_solution(incInfo) case(FIELD_THERMAL_ID) @@ -431,9 +386,9 @@ program DAMASK_grid case(FIELD_DAMAGE_ID) solres(field) = grid_damage_spectral_solution(timeinc) end select - + if (.not. solres(field)%converged) exit ! no solution found - + enddo stagIter = stagIter + 1 stagIterate = stagIter < stagItMax & @@ -467,38 +422,38 @@ program DAMASK_grid if (worldrank == 0) close(statUnit) call quit(0) ! quit endif - + enddo subStepLooping - + cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc - + if (all(solres(:)%converged)) then print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged' else print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged' endif; flush(IO_STDOUT) - - if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency + + if (mod(inc,loadCases(l)%outputFrequency) == 0) then ! at output frequency print'(1/,a)', ' ... writing results to file ......................................' flush(IO_STDOUT) call CPFEM_results(totalIncsCounter,time) endif - if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then + if (mod(inc,loadCases(l)%restartFrequency) == 0) then call mech_restartWrite call CPFEM_restartWrite endif endif skipping - + enddo incLooping - + enddo loadCaseLooping - - + + !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation print'(/,a)', ' ###########################################################################' if (worldrank == 0) close(statUnit) - + call quit(0) ! no complains ;) - + end program DAMASK_grid diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index b680dafa9..dfbd3f2f3 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -342,13 +342,13 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, !----------------------------------------------------------------------------------------------- ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F F_aimDot = F_aimDot & + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + elseif(deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed F_aimDot = F_aimDot & + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) - elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed F_aimDot = F_aimDot & + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif @@ -370,9 +370,9 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc - if (stress_BC%myType=='p') then + if (stress_BC%myType=='P') then P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) - elseif (stress_BC%myType=='pdot') then !UNTESTED + elseif (stress_BC%myType=='dot_P') then !UNTESTED P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) endif diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index e7e544886..d2f6b40da 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -308,13 +308,13 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !----------------------------------------------------------------------------------------------- ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F F_aimDot = F_aimDot & + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + elseif(deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed F_aimDot = F_aimDot & + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) - elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed F_aimDot = F_aimDot & + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif @@ -330,13 +330,13 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc - if (stress_BC%myType=='p') then + if (stress_BC%myType=='P') then P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) - elseif (stress_BC%myType=='pdot') then !UNTESTED + elseif (stress_BC%myType=='dot_P') then !UNTESTED P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) endif - F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index f59b68d7a..3d495ddf0 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -344,13 +344,13 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc !----------------------------------------------------------------------------------------------- ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F F_aimDot = F_aimDot & + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + elseif(deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed F_aimDot = F_aimDot & + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) - elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed F_aimDot = F_aimDot & + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif @@ -370,9 +370,9 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc - if (stress_BC%myType=='p') then + if (stress_BC%myType=='P') then P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) - elseif (stress_BC%myType=='pdot') then !UNTESTED + elseif (stress_BC%myType=='dot_P') then !UNTESTED P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) endif diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 4e01904e4..77047a317 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -82,21 +82,21 @@ module spectral_utilities end type tSolutionState type, public :: tBoundaryCondition !< set of parameters defining a boundary condition - real(pReal), dimension(3,3) :: values = 0.0_pReal - logical, dimension(3,3) :: mask = .false. - character(len=pStringLen) :: myType = 'None' + real(pReal), dimension(3,3) :: values = 0.0_pReal + logical, dimension(3,3) :: mask = .false. + character(len=:), allocatable :: myType end type tBoundaryCondition type, public :: tLoadCase type(rotation) :: rot !< rotation of BC type(tBoundaryCondition) :: stress, & !< stress BC - deformation !< deformation BC (Fdot or L) - real(pReal) :: time = 0.0_pReal !< length of increment - integer :: incs = 0, & !< number of increments - outputfrequency = 1, & !< frequency of result writes - restartfrequency = huge(0), & !< frequency of restart writes - logscale = 0 !< linear/logarithmic time inc flag - logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase + deformation !< deformation BC (dot_F, F, or L) + real(pReal) :: time !< length of increment + integer :: incs, & !< number of increments + outputfrequency, & !< frequency of result writes + restartfrequency !< frequency of restart writes + logical :: followFormerTrajectory, & !< follow trajectory of former loadcase + logscale !< logarithmic time inc flag integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) end type tLoadCase