diff --git a/PRIVATE b/PRIVATE index 25ce39dd0..05024077b 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 25ce39dd0f5bf49dc5c2bec20767d93d2b76d353 +Subproject commit 05024077bc89350d0313fc33d695a25612ac5130 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 8158996d6..a1bb23375 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -181,7 +181,7 @@ program DAMASK_grid 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='stress' + newLoadCase%stress%myType='p' field = 1 newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then @@ -210,18 +210,16 @@ program DAMASK_grid 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 enddo - newLoadCase%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) ! logical mask in 3x3 notation - newLoadCase%deformation%maskFloat = merge(ones,zeros,newLoadCase%deformation%maskLogical) ! float (1.0/0.0) mask in 3x3 notation - newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation + 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') temp_valueVector = 0.0_pReal 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 enddo - newLoadCase%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) - newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical) - newLoadCase%stress%values = math_9to33(temp_valueVector) + 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 @@ -268,8 +266,8 @@ program DAMASK_grid print*, ' drop guessing along trajectory' if (newLoadCase%deformation%myType == 'l') then do j = 1, 3 - if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & - any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined + 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 enddo print*, ' velocity gradient:' else if (newLoadCase%deformation%myType == 'f') then @@ -278,20 +276,19 @@ program DAMASK_grid print*, ' deformation gradient rate:' endif do i = 1, 3; do j = 1, 3 - if(newLoadCase%deformation%maskLogical(i,j)) then + if(newLoadCase%deformation%mask(i,j)) then write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j) else write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' endif enddo; write(IO_STDOUT,'(/)',advance='no') enddo - if (any(newLoadCase%stress%maskLogical .eqv. & - newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only - if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) & - .and. (math_I3<1))) errorID = 838 ! no rotation is allowed by stress BC + 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))) & + errorID = 838 ! no rotation is allowed by stress BC print*, ' stress / GPa:' do i = 1, 3; do j = 1, 3 - if(newLoadCase%stress%maskLogical(i,j)) then + if(newLoadCase%stress%mask(i,j)) then write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal else write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' @@ -368,11 +365,7 @@ program DAMASK_grid timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) else if (currentLoadCase == 1) then ! 1st load case of logarithmic scale - if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd - else ! not-1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal)) - endif + 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)/& @@ -432,17 +425,11 @@ program DAMASK_grid do field = 1, nActiveFields select case(loadCases(currentLoadCase)%ID(field)) case(FIELD_MECH_ID) - solres(field) = mech_solution (& - incInfo,timeinc,timeIncOld, & - stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rot) - + solres(field) = mech_solution(incInfo) case(FIELD_THERMAL_ID) - solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld) - + solres(field) = grid_thermal_spectral_solution(timeinc) case(FIELD_DAMAGE_ID) - solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld) - + solres(field) = grid_damage_spectral_solution(timeinc) end select if (.not. solres(field)%converged) exit ! no solution found diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 27dc5c1bd..7e529fc68 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -156,11 +156,10 @@ end subroutine grid_damage_spectral_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral damage scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) +function grid_damage_spectral_solution(timeinc) result(solution) real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment + timeinc !< increment in time for current solution integer :: i, j, k, cell type(tSolutionState) :: solution PetscInt :: devNull @@ -174,7 +173,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%timeinc = timeinc - params%timeincOld = timeinc_old call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index e1ceb42b0..b680dafa9 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -25,63 +25,56 @@ module grid_mech_FEM implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types - type(tSolutionParams), private :: params + type(tSolutionParams) :: params - type, private :: tNumerics + type :: tNumerics integer :: & itmin, & !< minimum number of iterations itmax !< maximum number of iterations real(pReal) :: & - err_div, & - divTol, & - BCTol, & eps_div_atol, & !< absolute tolerance for equilibrium eps_div_rtol, & !< relative tolerance for equilibrium eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC eps_stress_rtol !< relative tolerance for fullfillment of stress BC end type tNumerics - type(tNumerics), private :: num - logical, private:: & - debugRotation + type(tNumerics) :: num ! numerics parameters. Better name? + + logical :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: mech_grid - SNES, private :: mech_snes - Vec, private :: solution_current, solution_lastInc, solution_rate + DM :: mech_grid + SNES :: mech_snes + Vec :: solution_current, solution_lastInc, solution_rate !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc - real(pReal), private :: detJ - real(pReal), private, dimension(3) :: delta - real(pReal), private, dimension(3,8) :: BMat - real(pReal), private, dimension(8,8) :: HGMat - PetscInt, private :: xstart,ystart,zstart,xend,yend,zend + real(pReal), dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc + real(pReal) :: detJ + real(pReal), dimension(3) :: delta + real(pReal), dimension(3,8) :: BMat + real(pReal), dimension(8,8) :: HGMat + PetscInt :: xstart,ystart,zstart,xend,yend,zend !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & + real(pReal), dimension(3,3) :: & F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient - F_aim_lastIter = math_I3, & F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - - character(len=:), allocatable, private :: incInfo !< time and increment information - - real(pReal), private, dimension(3,3,3,3) :: & + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal + character(len=:), allocatable :: incInfo !< time and increment information + real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness S = 0.0_pReal !< current compliance (filled up with zeros) - real(pReal), private :: & + real(pReal) :: & err_BC !< deviation from stress BC - integer, private :: & + integer :: & totalIter = 0 !< total iteration in current increment public :: & @@ -99,7 +92,6 @@ contains subroutine grid_mech_FEM_init real(pReal) :: HGCoeff = 0.0e-2_pReal - PetscInt, dimension(0:worldsize-1) :: localK real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal real(pReal), dimension(4,8) :: & @@ -111,33 +103,34 @@ subroutine grid_mech_FEM_init -1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, & 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) + real(pReal), dimension(3,3,3,3) :: devNull PetscErrorCode :: ierr + PetscScalar, pointer, dimension(:,:,:,:) :: & + u_current,u_lastInc + PetscInt, dimension(0:worldsize-1) :: localK integer(HID_T) :: fileHandle, groupHandle character(len=pStringLen) :: & fileName class(tNode), pointer :: & num_grid, & debug_grid - real(pReal), dimension(3,3,3,3) :: devNull - PetscScalar, pointer, dimension(:,:,:,:) :: & - u_current,u_lastInc print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT) -!----------------------------------------------------------------------------------------------- +!------------------------------------------------------------------------------------------------- ! debugging options debug_grid => config_debug%get('grid', defaultVal=emptyList) debugRotation = debug_grid%contains('rotation') - + !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) - num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) - num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal) - num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) - num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) - num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) + num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal) + num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') @@ -242,6 +235,7 @@ subroutine grid_mech_FEM_init F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restartRead + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 @@ -268,19 +262,12 @@ end subroutine grid_mech_FEM_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the FEM scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) +function grid_mech_FEM_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - type(rotation), intent(in) :: & - rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -292,14 +279,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) -!-------------------------------------------------------------------------------------------------- -! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -341,6 +321,14 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%mask + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) @@ -349,20 +337,20 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, else C_volAvgLastInc = C_volAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + 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 - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + 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 - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif if (guess) then @@ -382,6 +370,12 @@ 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 + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) + endif + call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr) @@ -497,8 +491,6 @@ subroutine formResidual(da_local,x_local, & PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal PetscScalar, dimension(8,3) :: x_elem, f_elem PetscInt :: i, ii, j, jj, k, kk, ctr, ele - real(pReal), dimension(3,3) :: & - deltaF_aim PetscInt :: & PETScIter, & nfuncs @@ -547,10 +539,8 @@ subroutine formResidual(da_local,x_local, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim_lastIter = F_aim - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) - F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc + err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index e7fb9bd19..e7e544886 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -24,11 +24,9 @@ module grid_mech_spectral_basic implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params - type, private :: tNumerics + type :: tNumerics logical :: update_gamma !< update gamma operator with current stiffness integer :: & itmin, & !< minimum number of iterations @@ -42,7 +40,7 @@ module grid_mech_spectral_basic type(tNumerics) :: num ! numerics parameters. Better name? - logical, private :: debugRotation + logical :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -62,10 +60,10 @@ module grid_mech_spectral_basic F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information - real(pReal), private, dimension(3,3,3,3) :: & + real(pReal), 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 @@ -96,18 +94,17 @@ subroutine grid_mech_spectral_basic_init real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal - class (tNode), pointer :: & - num_grid, & - debug_grid - PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & F ! pointer to solution data - PetscInt, dimension(worldsize) :: localK + PetscInt, dimension(0:worldsize-1) :: localK integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: & fileName + class (tNode), pointer :: & + num_grid, & + debug_grid print'(/,a)', ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(IO_STDOUT) @@ -126,13 +123,13 @@ subroutine grid_mech_spectral_basic_init ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) - num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) - num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) - num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) - num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.01_pReal) - num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) - num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) + num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) + num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=1.0e-3_pReal) + num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') @@ -158,7 +155,7 @@ subroutine grid_mech_spectral_basic_init call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) localK = 0 - localK(worldrank+1) = grid3 + localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary @@ -202,8 +199,8 @@ subroutine grid_mech_spectral_basic_init endif restartRead materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateCoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_updateCoords(reshape(F,shape(F_lastInc))) + call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer @@ -231,19 +228,12 @@ end subroutine grid_mech_spectral_basic_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the basic scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) +function grid_mech_spectral_basic_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - type(rotation), intent(in) :: & - rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -255,17 +245,9 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) -!-------------------------------------------------------------------------------------------------- -! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old - !-------------------------------------------------------------------------------------------------- ! solve BVP call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) @@ -303,7 +285,14 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr - PetscScalar, dimension(:,:,:,:), pointer :: F + PetscScalar, pointer, dimension(:,:,:,:) :: F + +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%mask + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -314,20 +303,20 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + 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 - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + 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 - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -341,7 +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 - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + 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 + 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 rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -375,7 +370,7 @@ subroutine grid_mech_spectral_basic_restartWrite call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) - print'(a)', ' writing solver data required for restart to file'; flush(IO_STDOUT) + print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName,'w') @@ -469,6 +464,7 @@ subroutine formResidual(in, F, & call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment + !-------------------------------------------------------------------------------------------------- ! begin of new iteration newIteration: if (totalIter <= PETScIter) then @@ -491,16 +487,16 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) + deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! 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) = residuum ! store fPK field for subsequent FFT forward transform call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" - err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use + err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 5cc7d1602..f59b68d7a 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -15,7 +15,6 @@ module grid_mech_spectral_polarisation use DAMASK_interface use HDF5_utilities use math - use rotations use spectral_utilities use FEsolving use config @@ -25,8 +24,6 @@ module grid_mech_spectral_polarisation implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params type :: tNumerics @@ -48,7 +45,7 @@ module grid_mech_spectral_polarisation type(tNumerics) :: num ! numerics parameters. Better name? - logical, private :: debugRotation + logical :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -71,8 +68,8 @@ module grid_mech_spectral_polarisation 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_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness @@ -108,10 +105,6 @@ subroutine grid_mech_spectral_polarisation_init real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal - class (tNode), pointer :: & - num_grid, & - debug_grid - PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & FandF_tau, & ! overall pointer to solution data @@ -122,13 +115,16 @@ subroutine grid_mech_spectral_polarisation_init integer :: fileUnit character(len=pStringLen) :: & fileName + class (tNode), pointer :: & + num_grid, & + debug_grid print'(/,a)', ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(IO_STDOUT) print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' -!------------------------------------------------------------------------------------------------ +!------------------------------------------------------------------------------------------------- ! debugging options debug_grid => config_debug%get('grid',defaultVal=emptyList) debugRotation = debug_grid%contains('rotation') @@ -136,17 +132,18 @@ subroutine grid_mech_spectral_polarisation_init !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) - num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) - num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) - num%eps_curl_atol = num_grid%get_asFloat ('eps_curl_atol', defaultVal=1.0e-10_pReal) - num%eps_curl_rtol = num_grid%get_asFloat ('eps_curl_rtol', defaultVal=5.0e-4_pReal) - num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal) - num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) - num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) - num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) - num%alpha = num_grid%get_asFloat ('alpha', defaultVal=1.0_pReal) - num%beta = num_grid%get_asFloat ('beta', defaultVal=1.0_pReal) + + num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) + num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_curl_atol = num_grid%get_asFloat('eps_curl_atol', defaultVal=1.0e-10_pReal) + num%eps_curl_rtol = num_grid%get_asFloat('eps_curl_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal) + num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) + num%alpha = num_grid%get_asFloat('alpha', defaultVal=1.0_pReal) + num%beta = num_grid%get_asFloat('beta', defaultVal=1.0_pReal) if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') @@ -228,8 +225,8 @@ subroutine grid_mech_spectral_polarisation_init endif restartRead materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateCoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_updateCoords(reshape(F,shape(F_lastInc))) + call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer @@ -259,19 +256,12 @@ end subroutine grid_mech_spectral_polarisation_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) +function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - type(rotation), intent(in) :: & - rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -283,21 +273,13 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (num%update_gamma) then + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) + if(num%update_gamma) then call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) endif -!-------------------------------------------------------------------------------------------------- -! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old - !-------------------------------------------------------------------------------------------------- ! solve BVP call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) @@ -335,10 +317,17 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr - PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau + PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau integer :: i, j, k real(pReal), dimension(3,3) :: F_lambda33 +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%mask + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) F => FandF_tau(0: 8,:,:,:) F_tau => FandF_tau(9:17,:,:,:) @@ -350,20 +339,20 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + 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 - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + 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 - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -381,6 +370,12 @@ 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 + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) + elseif (stress_BC%myType=='pdot') 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 rotation_BC%rotate(F_aim,active=.true.)),& [9,grid(1),grid(2),grid3]) @@ -595,15 +590,15 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & - -params%rotation_BC%rotate(F_av)) + & - params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc + err_BC = maxval(abs(merge(P_av-P_aim, & + math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),& + params%stress_mask))) ! calculate divergence tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise call utilities_FFTtensorForward - err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress + err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -622,7 +617,7 @@ subroutine formResidual(in, FandF_tau, & 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() + err_curl = utilities_curlRMS() end subroutine formResidual diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 49be5ad7e..b4f9acddb 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -148,11 +148,10 @@ end subroutine grid_thermal_spectral_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) +function grid_thermal_spectral_solution(timeinc) result(solution) real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment + timeinc !< increment in time for current solution integer :: i, j, k, cell type(tSolutionState) :: solution PetscInt :: devNull @@ -166,7 +165,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%timeinc = timeinc - params%timeincOld = timeinc_old call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index dc6430c72..73aaa7789 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -47,15 +47,15 @@ module spectral_utilities complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field fourier representation for fftw real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for fftw complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field fourier representation for fftw - complex(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method - complex(pReal), private, dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives - complex(pReal), private, dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives - real(pReal), private, dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness + complex(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method + complex(pReal), dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives + complex(pReal), dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives + real(pReal), dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness !-------------------------------------------------------------------------------------------------- ! plans for FFTW - type(C_PTR), private :: & + type(C_PTR) :: & planTensorForth, & !< FFTW MPI plan P(x) to P(k) planTensorBack, & !< FFTW MPI plan F(k) to F(x) planVectorForth, & !< FFTW MPI plan v(x) to v(k) @@ -65,7 +65,7 @@ module spectral_utilities !-------------------------------------------------------------------------------------------------- ! variables controlling debugging - logical, private :: & + logical :: & debugGeneral, & !< general debugging of spectral solver debugRotation, & !< also printing out results in lab frame debugPETSc !< use some in debug defined options for more verbose PETSc solution @@ -82,10 +82,9 @@ 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, & - maskFloat = 0.0_pReal - logical, dimension(3,3) :: maskLogical = .false. - character(len=pStringLen) :: myType = 'None' + real(pReal), dimension(3,3) :: values = 0.0_pReal + logical, dimension(3,3) :: mask = .false. + character(len=pStringLen) :: myType = 'None' end type tBoundaryCondition type, public :: tLoadCase @@ -101,21 +100,21 @@ module spectral_utilities integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) end type tLoadCase - type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase - real(pReal), dimension(3,3) :: stress_mask, stress_BC + type, public :: tSolutionParams + real(pReal), dimension(3,3) :: stress_BC + logical, dimension(3,3) :: stress_mask type(rotation) :: rotation_BC - real(pReal) :: timeinc - real(pReal) :: timeincOld + real(pReal) :: timeinc, timeincOld end type tSolutionParams - type, private :: tNumerics + type :: tNumerics integer :: & divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints] logical :: & memory_efficient !< calculate gamma operator on the fly end type tNumerics - type(tNumerics), private :: num ! numerics parameters. Better name? + type(tNumerics) :: num ! numerics parameters. Better name? enum, bind(c); enumerator :: & DERIVATIVE_CONTINUOUS_ID, & diff --git a/src/material.f90 b/src/material.f90 index 6688a0ae5..b6d043183 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -2,7 +2,7 @@ !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Parses material config file, either solverJobName.materialConfig or material.config +!> @brief Defines phase and homogenization !-------------------------------------------------------------------------------------------------- module material use prec @@ -83,7 +83,7 @@ module material material_homogenizationAt !< homogenization ID of each element integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack material_homogenizationMemberAt !< position of the element within its homogenization instance - integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) + integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) material_phaseAt !< phase ID of each element integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,elem) material_phaseMemberAt !< position of the element within its phase instance @@ -326,7 +326,7 @@ subroutine material_parseMicrostructure constituents, & !> list of constituents constituent, & !> constituent definition phases, & - homogenization + homogenizations integer, dimension(:), allocatable :: & counterPhase, & @@ -362,14 +362,14 @@ subroutine material_parseMicrostructure phases => config_material%get('phase') allocate(counterPhase(phases%length),source=0) - homogenization => config_material%get('homogenization') - allocate(counterHomogenization(homogenization%length),source=0) + homogenizations => config_material%get('homogenization') + allocate(counterHomogenization(homogenizations%length),source=0) do e = 1, discretization_nElem microstructure => microstructures%get(discretization_microstructureAt(e)) constituents => microstructure%get('constituents') - material_homogenizationAt(e) = homogenization%getIndex(microstructure%get_asString('homogenization')) + material_homogenizationAt(e) = homogenizations%getIndex(microstructure%get_asString('homogenization')) do i = 1, discretization_nIP counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1 material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e)) diff --git a/src/rotations.f90 b/src/rotations.f90 index e19513866..0b72c1dd5 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -1452,7 +1452,7 @@ subroutine selfTest contains - function quaternion_equal(qu1,qu2) result(ok) + pure recursive function quaternion_equal(qu1,qu2) result(ok) real(pReal), intent(in), dimension(4) :: qu1,qu2 logical :: ok