223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Need only logical mask 

'merge' substitutes multiplication with float mask
This commit is contained in:
Martin Diehl 2020-09-20 16:24:08 +02:00
parent 6367cb8fcb
commit 0a7d4f61ac
5 changed files with 96 additions and 106 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

105
src/grid/DAMASK_grid.f90
View File

@ -24,7 +24,7 @@ program DAMASK_grid
use grid_damage_spectral
use grid_thermal_spectral
use results
implicit none
!--------------------------------------------------------------------------------------------------
@ -88,7 +88,7 @@ program DAMASK_grid
mech_updateCoords
procedure(grid_mech_spectral_basic_restartWrite), pointer :: &
mech_restartWrite
external :: &
quit
class (tNode), pointer :: &
@ -97,10 +97,10 @@ program DAMASK_grid
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
call CPFEM_initAll
call CPFEM_initAll
print'(/,a)', ' <<<+- DAMASK_spectral init -+>>>'; flush(6)
print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
@ -123,16 +123,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 +141,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,24 +150,24 @@ 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
! 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')
@ -180,7 +180,7 @@ program DAMASK_grid
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'
field = 1
newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default
@ -192,7 +192,7 @@ program DAMASK_grid
field = field + 1
newLoadCase%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)))
@ -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
@ -258,9 +256,9 @@ program DAMASK_grid
call newLoadCase%rot%fromMatrix(math_9to33(temp_valueVector))
end select
enddo readIn
newLoadCase%followFormerTrajectory = merge(.true.,.false.,currentLoadCase > 1) ! by default, guess from previous load case
reportAndCheck: if (worldrank == 0) then
write (loadcase_string, '(i0)' ) currentLoadCase
print'(/,a,i0)', ' load case: ', currentLoadCase
@ -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(6,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j)
else
write(6,'(2x,12a)',advance='no') ' * '
endif
enddo; write(6,'(/)',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(6,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal
else
write(6,'(2x,12a)',advance='no') ' * '
@ -325,13 +322,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
@ -348,16 +345,16 @@ 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)
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
totalIncsCounter = totalIncsCounter + 1
@ -382,13 +379,13 @@ program DAMASK_grid
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
time = time + timeinc ! forward target time
@ -417,7 +414,7 @@ program DAMASK_grid
deformation_BC = loadCases(currentLoadCase)%deformation, &
stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rot)
case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack)
end select
@ -438,9 +435,9 @@ program DAMASK_grid
case(FIELD_DAMAGE_ID)
solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld)
end select
if (.not. solres(field)%converged) exit ! no solution found
enddo
stagIter = stagIter + 1
stagIterate = stagIter < stagItMax &
@ -474,17 +471,17 @@ 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(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
print'(1/,a)', ' ... writing results to file ......................................'
flush(6)
@ -495,17 +492,17 @@ program DAMASK_grid
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

26
src/grid/grid_mech_FEM.f90
View File

@ -25,8 +25,6 @@ module grid_mech_FEM
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams) :: params
type :: tNumerics
@ -282,7 +280,7 @@ function grid_mech_FEM_solution(incInfoIn) result(solution)
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
!--------------------------------------------------------------------------------------------------
! solve BVP
@ -326,7 +324,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
@ -340,20 +338,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
@ -374,9 +372,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
P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc
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 + stress_BC%maskFloat*stress_BC%values*timeinc
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)
@ -543,7 +541,7 @@ subroutine formResidual(da_local,x_local, &
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs(params%stress_mask * (P_av - P_aim))) ! mask = 0.0 when no stress bc
err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask)))
!--------------------------------------------------------------------------------------------------
! constructing residual

26
src/grid/grid_mech_spectral_basic.f90
View File

@ -24,8 +24,6 @@ module grid_mech_spectral_basic
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams) :: params
type :: tNumerics
@ -247,7 +245,7 @@ function grid_mech_spectral_basic_solution(incInfoIn) result(solution)
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
!--------------------------------------------------------------------------------------------------
@ -291,7 +289,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
@ -305,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, &
@ -333,9 +331,9 @@ 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
P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc
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 + stress_BC%maskFloat*stress_BC%values*timeinc
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
@ -491,7 +489,7 @@ subroutine formResidual(in, F, &
! stress BC handling
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 - P_aim))) ! 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

30
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -24,8 +24,6 @@ module grid_mech_spectral_polarisation
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams) :: params
type :: tNumerics
@ -275,7 +273,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution)
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg)
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
@ -325,7 +323,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
@ -341,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, &
@ -373,9 +371,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
P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc
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 + stress_BC%maskFloat*stress_BC%values*timeinc
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
@ -593,9 +591,9 @@ subroutine formResidual(in, FandF_tau, &
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! 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-P_aim))) ! mask = 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

15
src/grid/spectral_utilities.f90
View File

@ -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,11 +100,11 @@ 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 :: tNumerics