consistent variants including space
This commit is contained in:
Martin Diehl
parent
5a284648ed
commit
78838c2f75
26
src/CLI.f90
26
src/CLI.f90
|
|
@ -156,15 +156,15 @@ subroutine CLI_init
|
|||
if (CLI_restartInc < 0 .or. stat /=0) then
|
||||
print'(/,a)', ' ERROR: Could not parse restart increment: '//trim(arg)
|
||||
call quit(1)
|
||||
endif
|
||||
end if
|
||||
end select
|
||||
if (err /= 0) call quit(1)
|
||||
enddo
|
||||
end do
|
||||
|
||||
if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then
|
||||
print'(/,a)', ' ERROR: Please specify geometry AND load case (-h for help)'
|
||||
call quit(1)
|
||||
endif
|
||||
end if
|
||||
|
||||
if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg))
|
||||
CLI_geomFile = getGeometryFile(geometryArg)
|
||||
|
|
@ -205,14 +205,14 @@ subroutine setWorkingDirectory(workingDirectoryArg)
|
|||
else absolutePath
|
||||
workingDirectory = getCWD()
|
||||
workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg
|
||||
endif absolutePath
|
||||
end if absolutePath
|
||||
|
||||
workingDirectory = trim(rectifyPath(workingDirectory))
|
||||
error = setCWD(trim(workingDirectory))
|
||||
if(error) then
|
||||
print*, 'ERROR: Invalid Working directory: '//trim(workingDirectory)
|
||||
call quit(1)
|
||||
endif
|
||||
end if
|
||||
|
||||
end subroutine setWorkingDirectory
|
||||
|
||||
|
|
@ -256,7 +256,7 @@ function getGeometryFile(geometryParameter)
|
|||
if (.not. file_exists) then
|
||||
print*, 'ERROR: Geometry file does not exists: '//trim(getGeometryFile)
|
||||
call quit(1)
|
||||
endif
|
||||
end if
|
||||
|
||||
end function getGeometryFile
|
||||
|
||||
|
|
@ -279,7 +279,7 @@ function getLoadCaseFile(loadCaseParameter)
|
|||
if (.not. file_exists) then
|
||||
print*, 'ERROR: Load case file does not exists: '//trim(getLoadCaseFile)
|
||||
call quit(1)
|
||||
endif
|
||||
end if
|
||||
|
||||
end function getLoadCaseFile
|
||||
|
||||
|
|
@ -300,14 +300,14 @@ function rectifyPath(path)
|
|||
l = len_trim(rectifyPath)
|
||||
do i = l,3,-1
|
||||
if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
|
||||
enddo
|
||||
end do
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! remove // from path
|
||||
l = len_trim(rectifyPath)
|
||||
do i = l,2,-1
|
||||
if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' '
|
||||
enddo
|
||||
end do
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! remove ../ and corresponding directory from rectifyPath
|
||||
|
|
@ -321,9 +321,9 @@ function rectifyPath(path)
|
|||
k = len_trim(rectifyPath)
|
||||
rectifyPath(j+1:k-1) = rectifyPath(j+2:k)
|
||||
rectifyPath(k:k) = ' '
|
||||
endif
|
||||
end if
|
||||
i = j+index(rectifyPath(j+1:l),'../')
|
||||
enddo
|
||||
end do
|
||||
if(len_trim(rectifyPath) == 0) rectifyPath = '/'
|
||||
|
||||
rectifyPath = trim(rectifyPath)
|
||||
|
|
@ -349,10 +349,10 @@ function makeRelativePath(a,b)
|
|||
do i = 1, min(len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned)))
|
||||
if (a_cleaned(i:i) /= b_cleaned(i:i)) exit
|
||||
if (a_cleaned(i:i) == '/') posLastCommonSlash = i
|
||||
enddo
|
||||
end do
|
||||
do i = posLastCommonSlash+1,len_trim(a_cleaned)
|
||||
if (a_cleaned(i:i) == '/') remainingSlashes = remainingSlashes + 1
|
||||
enddo
|
||||
end do
|
||||
|
||||
makeRelativePath = repeat('..'//'/',remainingSlashes)//b_cleaned(posLastCommonSlash+1:len_trim(b_cleaned))
|
||||
|
||||
|
|
|
|||
|
|
@ -509,7 +509,7 @@ subroutine HDF5_addAttribute_str_array(loc_id,attrLabel,attrValue,path)
|
|||
do i=1,size(attrValue)
|
||||
attrValue_(i) = attrValue(i)//C_NULL_CHAR
|
||||
ptr(i) = c_loc(attrValue_(i))
|
||||
enddo
|
||||
end do
|
||||
|
||||
call H5Screate_simple_f(1,shape(attrValue_,kind=HSIZE_T),space_id,hdferr,shape(attrValue_,kind=HSIZE_T))
|
||||
if(hdferr < 0) error stop 'HDF5 error'
|
||||
|
|
|
|||
34
src/IO.f90
34
src/IO.f90
|
|
@ -92,15 +92,15 @@ function IO_readlines(fileName) result(fileContent)
|
|||
if (.not. warned) then
|
||||
call IO_warning(207,trim(fileName),label1='line',ID1=l)
|
||||
warned = .true.
|
||||
endif
|
||||
end if
|
||||
else
|
||||
line = rawData(startPos:endpos)
|
||||
endif
|
||||
end if
|
||||
startPos = endPos + 2 ! jump to next line start
|
||||
|
||||
fileContent(l) = trim(line)//''
|
||||
l = l + 1
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function IO_readlines
|
||||
|
||||
|
|
@ -129,7 +129,7 @@ function IO_read(fileName) result(fileContent)
|
|||
if (fileLength==0) then
|
||||
close(fileUnit)
|
||||
return
|
||||
endif
|
||||
end if
|
||||
|
||||
read(fileUnit,iostat=myStat) fileContent
|
||||
if (myStat /= 0) call IO_error(102,trim(fileName))
|
||||
|
|
@ -183,8 +183,8 @@ pure function IO_stringPos(string)
|
|||
endOfString: if (right < left) then
|
||||
IO_stringPos(IO_stringPos(1)*2+1) = len_trim(string)
|
||||
exit
|
||||
endif endOfString
|
||||
enddo
|
||||
end if endOfString
|
||||
end do
|
||||
|
||||
end function IO_stringPos
|
||||
|
||||
|
|
@ -204,7 +204,7 @@ function IO_stringValue(string,chunkPos,myChunk)
|
|||
call IO_error(110,'IO_stringValue: "'//trim(string)//'"',label1='chunk',ID1=myChunk)
|
||||
else validChunk
|
||||
IO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1))
|
||||
endif validChunk
|
||||
end if validChunk
|
||||
|
||||
end function IO_stringValue
|
||||
|
||||
|
|
@ -257,8 +257,8 @@ pure function IO_lc(string)
|
|||
IO_lc(i:i) = LOWER(n:n)
|
||||
else
|
||||
IO_lc(i:i) = string(i:i)
|
||||
endif
|
||||
enddo
|
||||
end if
|
||||
end do
|
||||
|
||||
end function IO_lc
|
||||
|
||||
|
|
@ -280,7 +280,7 @@ function IO_rmComment(line)
|
|||
IO_rmComment = trim(line)
|
||||
else
|
||||
IO_rmComment = trim(line(:split-1))
|
||||
endif
|
||||
end if
|
||||
|
||||
end function IO_rmComment
|
||||
|
||||
|
|
@ -302,7 +302,7 @@ integer function IO_stringAsInt(string)
|
|||
else valid
|
||||
IO_stringAsInt = 0
|
||||
call IO_error(111,string)
|
||||
endif valid
|
||||
end if valid
|
||||
|
||||
end function IO_stringAsInt
|
||||
|
||||
|
|
@ -324,7 +324,7 @@ real(pReal) function IO_stringAsFloat(string)
|
|||
else valid
|
||||
IO_stringAsFloat = 0.0_pReal
|
||||
call IO_error(112,string)
|
||||
endif valid
|
||||
end if valid
|
||||
|
||||
end function IO_stringAsFloat
|
||||
|
||||
|
|
@ -344,7 +344,7 @@ logical function IO_stringAsBool(string)
|
|||
else
|
||||
IO_stringAsBool = .false.
|
||||
call IO_error(113,string)
|
||||
endif
|
||||
end if
|
||||
|
||||
end function IO_stringAsBool
|
||||
|
||||
|
|
@ -598,7 +598,7 @@ pure function CRLF2LF(string)
|
|||
CRLF2LF(c-n:c-n) = string(c:c)
|
||||
if (c == len_trim(string)) exit
|
||||
if (string(c:c+1) == CR//LF) n = n + 1
|
||||
enddo
|
||||
end do
|
||||
|
||||
CRLF2LF = CRLF2LF(:c-n)
|
||||
|
||||
|
|
@ -639,17 +639,17 @@ subroutine panel(paneltype,ID,msg,ext_msg,label1,ID1,label2,ID2)
|
|||
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
|
||||
max(1,panelwidth+3-len_trim(ext_msg)-4),'x,a)'
|
||||
write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│'
|
||||
endif
|
||||
end if
|
||||
if (present(label1)) then
|
||||
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label1)),',i9,',&
|
||||
max(1,panelwidth+3-len_trim(label1)-9-7),'x,a)'
|
||||
write(IO_STDERR,formatString) '│ at ',trim(label1),ID1, '│'
|
||||
endif
|
||||
end if
|
||||
if (present(label2)) then
|
||||
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label2)),',i9,',&
|
||||
max(1,panelwidth+3-len_trim(label2)-9-7),'x,a)'
|
||||
write(IO_STDERR,formatString) '│ at ',trim(label2),ID2, '│'
|
||||
endif
|
||||
end if
|
||||
write(formatString,'(a,i2.2,a)') '(a,',max(1,panelwidth),'x,a)'
|
||||
write(IO_STDERR,formatString) ' │', '│'
|
||||
write(IO_STDERR,'(a)') ' └'//DIVIDER//'┘'
|
||||
|
|
|
|||
|
|
@ -103,7 +103,7 @@ recursive function parse_flow(YAML_flow) result(node)
|
|||
class is (tDict)
|
||||
call node%set(key,myVal)
|
||||
end select
|
||||
enddo
|
||||
end do
|
||||
elseif (flow_string(1:1) == '[') then ! start of a list
|
||||
e = 1
|
||||
allocate(tList::node)
|
||||
|
|
@ -116,7 +116,7 @@ recursive function parse_flow(YAML_flow) result(node)
|
|||
class is (tList)
|
||||
call node%append(myVal)
|
||||
end select
|
||||
enddo
|
||||
end do
|
||||
else ! scalar value
|
||||
allocate(tScalar::node)
|
||||
select type (node)
|
||||
|
|
@ -156,7 +156,7 @@ integer function find_end(str,e_char)
|
|||
N_sq = N_sq - merge(1,0,str(i:i) == ']')
|
||||
N_cu = N_cu - merge(1,0,str(i:i) == '}')
|
||||
i = i + 1
|
||||
enddo
|
||||
end do
|
||||
find_end = i
|
||||
|
||||
end function find_end
|
||||
|
|
@ -333,7 +333,7 @@ subroutine skip_empty_lines(blck,s_blck)
|
|||
do while(empty .and. len_trim(blck(s_blck:)) /= 0)
|
||||
empty = len_trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == 0
|
||||
if(empty) s_blck = s_blck + index(blck(s_blck:),IO_EOL)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end subroutine skip_empty_lines
|
||||
|
||||
|
|
@ -387,7 +387,7 @@ logical function flow_is_closed(str,e_char)
|
|||
N_cu = N_cu + merge(1,0,line(i:i) == '{')
|
||||
N_sq = N_sq - merge(1,0,line(i:i) == ']')
|
||||
N_cu = N_cu - merge(1,0,line(i:i) == '}')
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function flow_is_closed
|
||||
|
||||
|
|
@ -410,7 +410,7 @@ subroutine remove_line_break(blck,s_blck,e_char,flow_line)
|
|||
flow_line = flow_line//IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))//' '
|
||||
line_end = flow_is_closed(flow_line,e_char)
|
||||
s_blck = s_blck + index(blck(s_blck:),IO_EOL)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end subroutine remove_line_break
|
||||
|
||||
|
|
@ -439,7 +439,7 @@ subroutine list_item_inline(blck,s_blck,inline,offset)
|
|||
inline = inline//' '//trim(adjustl(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))))
|
||||
s_blck = s_blck + index(blck(s_blck:),IO_EOL)
|
||||
indent_next = indentDepth(blck(s_blck:))
|
||||
enddo
|
||||
end do
|
||||
|
||||
if(scan(inline,",") > 0) inline = '"'//inline//'"'
|
||||
|
||||
|
|
@ -481,7 +481,7 @@ recursive subroutine line_isFlow(flow,s_flow,line)
|
|||
flow(s_flow:s_flow+1) = ', '
|
||||
s_flow = s_flow +2
|
||||
s = s + find_end(line(s+1:),']')
|
||||
enddo
|
||||
end do
|
||||
s_flow = s_flow - 1
|
||||
if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
|
||||
flow(s_flow:s_flow) = ']'
|
||||
|
|
@ -498,7 +498,7 @@ recursive subroutine line_isFlow(flow,s_flow,line)
|
|||
flow(s_flow:s_flow+1) = ', '
|
||||
s_flow = s_flow +2
|
||||
s = s + find_end(line(s+1:),'}')
|
||||
enddo
|
||||
end do
|
||||
s_flow = s_flow -1
|
||||
if(flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow -1
|
||||
flow(s_flow:s_flow) = '}'
|
||||
|
|
@ -646,7 +646,7 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
|
|||
s_flow = s_flow + 2
|
||||
end if
|
||||
|
||||
enddo
|
||||
end do
|
||||
|
||||
s_flow = s_flow - 1
|
||||
if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
|
||||
|
|
@ -733,7 +733,7 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset)
|
|||
flow(s_flow:s_flow) = ' '
|
||||
s_flow = s_flow + 1
|
||||
offset = 0
|
||||
enddo
|
||||
end do
|
||||
|
||||
s_flow = s_flow - 1
|
||||
if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
|
||||
|
|
|
|||
|
|
@ -411,7 +411,7 @@ function tNode_get_byIndex(self,i) result(node)
|
|||
|
||||
do j = 2,i
|
||||
item => item%next
|
||||
enddo
|
||||
end do
|
||||
node => item%node
|
||||
|
||||
end function tNode_get_byIndex
|
||||
|
|
@ -681,7 +681,7 @@ function tNode_contains(self,k) result(exists)
|
|||
exists = .true.
|
||||
return
|
||||
end if
|
||||
enddo
|
||||
end do
|
||||
class is(tList)
|
||||
list => self%asList()
|
||||
do j=1, list%length
|
||||
|
|
@ -689,7 +689,7 @@ function tNode_contains(self,k) result(exists)
|
|||
exists = .true.
|
||||
return
|
||||
end if
|
||||
enddo
|
||||
end do
|
||||
class default
|
||||
call IO_error(706,ext_msg='Expected list or dict')
|
||||
end select
|
||||
|
|
@ -731,7 +731,7 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
|
|||
end if
|
||||
item => item%next
|
||||
j = j + 1
|
||||
enddo
|
||||
end do
|
||||
|
||||
if (.not. found) then
|
||||
call IO_error(143,ext_msg=k)
|
||||
|
|
@ -1333,7 +1333,7 @@ function tList_as1dString(self)
|
|||
scalar => item%node%asScalar()
|
||||
tList_as1dString(i) = scalar%asString()
|
||||
item => item%next
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function tList_as1dString
|
||||
|
||||
|
|
@ -1384,7 +1384,7 @@ subroutine tDict_set(self,key,node)
|
|||
searchExisting: do while (associated(item%next))
|
||||
if (item%key == key) exit
|
||||
item => item%next
|
||||
enddo searchExisting
|
||||
end do searchExisting
|
||||
if (item%key /= key) then
|
||||
allocate(item%next)
|
||||
item => item%next
|
||||
|
|
|
|||
|
|
@ -68,7 +68,7 @@ subroutine discretization_init(materialAt,&
|
|||
discretization_sharedNodesBegin = sharedNodesBegin
|
||||
else
|
||||
discretization_sharedNodesBegin = size(discretization_NodeCoords0,2)
|
||||
endif
|
||||
end if
|
||||
|
||||
end subroutine discretization_init
|
||||
|
||||
|
|
|
|||
|
|
@ -365,7 +365,7 @@ subroutine homogenization_results
|
|||
call thermal_results(ho,group)
|
||||
end if
|
||||
|
||||
enddo
|
||||
end do
|
||||
|
||||
end subroutine homogenization_results
|
||||
|
||||
|
|
@ -383,7 +383,7 @@ subroutine homogenization_forward
|
|||
homogState (ho)%state0 = homogState (ho)%state
|
||||
if(damageState_h(ho)%sizeState > 0) &
|
||||
damageState_h(ho)%state0 = damageState_h(ho)%state
|
||||
enddo
|
||||
end do
|
||||
|
||||
end subroutine homogenization_forward
|
||||
|
||||
|
|
@ -408,7 +408,7 @@ subroutine homogenization_restartWrite(fileHandle)
|
|||
|
||||
call HDF5_closeGroup(groupHandle(2))
|
||||
|
||||
enddo
|
||||
end do
|
||||
|
||||
call HDF5_closeGroup(groupHandle(1))
|
||||
|
||||
|
|
@ -435,7 +435,7 @@ subroutine homogenization_restartRead(fileHandle)
|
|||
|
||||
call HDF5_closeGroup(groupHandle(2))
|
||||
|
||||
enddo
|
||||
end do
|
||||
|
||||
call HDF5_closeGroup(groupHandle(1))
|
||||
|
||||
|
|
@ -476,7 +476,7 @@ subroutine parseHomogenization
|
|||
case default
|
||||
call IO_error(500,ext_msg=homogThermal%get_asString('type'))
|
||||
end select
|
||||
endif
|
||||
end if
|
||||
|
||||
if (homog%contains('damage')) then
|
||||
homogDamage => homog%get('damage')
|
||||
|
|
@ -486,8 +486,8 @@ subroutine parseHomogenization
|
|||
case default
|
||||
call IO_error(500,ext_msg=homogDamage%get_asString('type'))
|
||||
end select
|
||||
endif
|
||||
enddo
|
||||
end if
|
||||
end do
|
||||
|
||||
end subroutine parseHomogenization
|
||||
|
||||
|
|
|
|||
|
|
@ -65,9 +65,9 @@ module subroutine damage_init()
|
|||
allocate(damageState_h(ho)%state (1,Nmembers), source=1.0_pReal)
|
||||
else
|
||||
prm%output = emptyStringArray
|
||||
endif
|
||||
end if
|
||||
end associate
|
||||
enddo
|
||||
end do
|
||||
|
||||
call pass_init()
|
||||
|
||||
|
|
@ -174,7 +174,7 @@ module subroutine damage_results(ho,group)
|
|||
call results_writeDataset(damagestate_h(ho)%state(1,:),group,prm%output(o),&
|
||||
'damage indicator','-')
|
||||
end select
|
||||
enddo outputsLoop
|
||||
end do outputsLoop
|
||||
end associate
|
||||
|
||||
end subroutine damage_results
|
||||
|
|
|
|||
|
|
@ -561,7 +561,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
|
|||
*cosh(prm%c_alpha*nDefNorm) &
|
||||
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
|
||||
*tanh(nDefNorm/num%xSmoo)
|
||||
end do; end do;enddo; end do
|
||||
end do; end do;end do; end do
|
||||
end do interfaceLoop
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -484,8 +484,8 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
|
|||
case default
|
||||
call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice))
|
||||
end select
|
||||
enddo mySystems
|
||||
enddo myFamilies
|
||||
end do mySystems
|
||||
end do myFamilies
|
||||
|
||||
end function lattice_characteristicShear_Twin
|
||||
|
||||
|
|
@ -523,7 +523,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
|
|||
do i = 1, sum(Ntwin)
|
||||
call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg?
|
||||
lattice_C66_twin(1:6,1:6,i) = R%rotStiffness(C66)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function lattice_C66_twin
|
||||
|
||||
|
|
@ -572,19 +572,19 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
|
|||
C_target_unrotated66 = C_parent66
|
||||
else
|
||||
call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target))
|
||||
endif
|
||||
end if
|
||||
|
||||
do i = 1,6
|
||||
if (abs(C_target_unrotated66(i,i))<tol_math_check) &
|
||||
call IO_error(135,'matrix diagonal in transformation',label1='entry',ID1=i)
|
||||
enddo
|
||||
end do
|
||||
|
||||
call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_cF,a_cI)
|
||||
|
||||
do i = 1,sum(Ntrans)
|
||||
call R%fromMatrix(Q(1:3,1:3,i))
|
||||
lattice_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function lattice_C66_trans
|
||||
|
||||
|
|
@ -632,7 +632,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
|
|||
math_cross(normal, direction))
|
||||
if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
|
||||
+ nonSchmidCoefficients(6) * math_outer(direction, direction)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function lattice_nonSchmidMatrix
|
||||
|
||||
|
|
@ -1432,7 +1432,7 @@ function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
|
|||
SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
|
||||
if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) &
|
||||
error stop 'dilatational Schmid matrix for slip'
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function lattice_SchmidMatrix_slip
|
||||
|
||||
|
|
@ -1479,7 +1479,7 @@ function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
|
|||
SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
|
||||
if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) &
|
||||
error stop 'dilatational Schmid matrix for twin'
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function lattice_SchmidMatrix_twin
|
||||
|
||||
|
|
@ -1552,7 +1552,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMa
|
|||
SchmidMatrix(1:3,1:3,1,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
|
||||
SchmidMatrix(1:3,1:3,2,i) = math_outer(coordinateSystem(1:3,3,i),coordinateSystem(1:3,2,i))
|
||||
SchmidMatrix(1:3,1:3,3,i) = math_outer(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,i))
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function lattice_SchmidMatrix_cleavage
|
||||
|
||||
|
|
@ -1719,8 +1719,8 @@ pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym)
|
|||
do i = 1, 6
|
||||
do j = i+1, 6
|
||||
C66_sym(j,i) = C66_sym(i,j)
|
||||
enddo
|
||||
enddo
|
||||
end do
|
||||
end do
|
||||
|
||||
end function lattice_symmetrize_C66
|
||||
|
||||
|
|
@ -1782,7 +1782,7 @@ function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
|
|||
|
||||
do i=1, sum(Nslip); do j=1, sum(Nslip)
|
||||
projection(i,j) = abs(math_inner(n(:,i),t(:,j)))
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
|
||||
end function slipProjection_transverse
|
||||
|
||||
|
|
@ -1806,7 +1806,7 @@ function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
|
|||
|
||||
do i=1, sum(Nslip); do j=1, sum(Nslip)
|
||||
projection(i,j) = abs(math_inner(n(:,i),d(:,j)))
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
|
||||
end function slipProjection_direction
|
||||
|
||||
|
|
@ -1890,8 +1890,8 @@ function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,valu
|
|||
|
||||
buildInteraction(l,k) = values(matrix(i,j))
|
||||
|
||||
enddo; enddo
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
end do; end do
|
||||
|
||||
end function buildInteraction
|
||||
|
||||
|
|
@ -1957,8 +1957,8 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
|
|||
buildCoordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),&
|
||||
normal /norm2(normal))
|
||||
|
||||
enddo activeSystems
|
||||
enddo activeFamilies
|
||||
end do activeSystems
|
||||
end do activeFamilies
|
||||
|
||||
end function buildCoordinateSystem
|
||||
|
||||
|
|
@ -2066,7 +2066,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
|
|||
U = (a_cI/a_cF) * (math_outer(x,x) + (math_outer(y,y)+math_outer(z,z)) * sqrt(2.0_pReal))
|
||||
Q(1:3,1:3,i) = matmul(R%asMatrix(),B%asMatrix())
|
||||
S(1:3,1:3,i) = matmul(R%asMatrix(),U) - MATH_I3
|
||||
enddo
|
||||
end do
|
||||
else if (present(cOverA)) then
|
||||
ss = MATH_I3
|
||||
sd = MATH_I3
|
||||
|
|
@ -2125,7 +2125,7 @@ function getlabels(active,potential,system) result(labels)
|
|||
write(label(i+1:i+2),'(I2.1)') int(system(j,p))
|
||||
label(i+3:i+3) = ' '
|
||||
i = i + 3
|
||||
enddo direction
|
||||
end do direction
|
||||
label(i:i) = ']'
|
||||
|
||||
i = i +1
|
||||
|
|
@ -2134,13 +2134,13 @@ function getlabels(active,potential,system) result(labels)
|
|||
write(label(i+1:i+2),'(I2.1)') int(system(j,p))
|
||||
label(i+3:i+3) = ' '
|
||||
i = i + 3
|
||||
enddo normal
|
||||
end do normal
|
||||
label(i:i) = ')'
|
||||
|
||||
labels(a) = label
|
||||
|
||||
enddo activeSystems
|
||||
enddo activeFamilies
|
||||
end do activeSystems
|
||||
end do activeFamilies
|
||||
|
||||
end function getlabels
|
||||
|
||||
|
|
@ -2170,7 +2170,7 @@ pure function lattice_equivalent_nu(C,assumption) result(nu)
|
|||
/ (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3)))
|
||||
else
|
||||
error stop 'invalid assumption'
|
||||
endif
|
||||
end if
|
||||
|
||||
mu = lattice_equivalent_mu(C,assumption)
|
||||
nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu)
|
||||
|
|
@ -2202,7 +2202,7 @@ pure function lattice_equivalent_mu(C,assumption) result(mu)
|
|||
/ (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6)))
|
||||
else
|
||||
error stop 'invalid assumption'
|
||||
endif
|
||||
end if
|
||||
|
||||
end function lattice_equivalent_mu
|
||||
|
||||
|
|
@ -2266,7 +2266,7 @@ subroutine selfTest
|
|||
if (any(dNeq(T(1,1),[T_hP(1,1),T_hP(2,2)]))) error stop 'Symmetry33_11-22/hP'
|
||||
if (any(dNeq(T(1,1),[T_tI(1,1),T_tI(2,2)]))) error stop 'Symmetry33_11-22/tI'
|
||||
|
||||
enddo
|
||||
end do
|
||||
|
||||
call random_number(C)
|
||||
C(1,1) = C(1,1) + C(1,2) + 0.1_pReal
|
||||
|
|
|
|||
68
src/math.f90
68
src/math.f90
|
|
@ -135,25 +135,25 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
|
|||
s = istart
|
||||
else
|
||||
s = lbound(a,2)
|
||||
endif
|
||||
end if
|
||||
|
||||
if (present(iend)) then
|
||||
e = iend
|
||||
else
|
||||
e = ubound(a,2)
|
||||
endif
|
||||
end if
|
||||
|
||||
if (present(sortDim)) then
|
||||
d = sortDim
|
||||
else
|
||||
d = 1
|
||||
endif
|
||||
end if
|
||||
|
||||
if (s < e) then
|
||||
call qsort_partition(a,ipivot, s,e, d)
|
||||
call math_sort(a, s, ipivot-1, d)
|
||||
call math_sort(a, ipivot+1, e, d)
|
||||
endif
|
||||
end if
|
||||
|
||||
|
||||
contains
|
||||
|
|
@ -175,11 +175,11 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
|
|||
! find the first element on the right side less than or equal to the pivot point
|
||||
do j = iend, istart, -1
|
||||
if (a(sort,j) <= a(sort,istart)) exit
|
||||
enddo
|
||||
end do
|
||||
! find the first element on the left side greater than the pivot point
|
||||
do i = istart, iend
|
||||
if (a(sort,i) > a(sort,istart)) exit
|
||||
enddo
|
||||
end do
|
||||
cross: if (i >= j) then ! exchange left value with pivot and return with the partition index
|
||||
tmp = a(:,istart)
|
||||
a(:,istart) = a(:,j)
|
||||
|
|
@ -190,8 +190,8 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
|
|||
tmp = a(:,i)
|
||||
a(:,i) = a(:,j)
|
||||
a(:,j) = tmp
|
||||
endif cross
|
||||
enddo
|
||||
end if cross
|
||||
end do
|
||||
|
||||
end subroutine qsort_partition
|
||||
|
||||
|
|
@ -216,7 +216,7 @@ pure function math_expand(what,how)
|
|||
|
||||
do i = 1, size(how)
|
||||
math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function math_expand
|
||||
|
||||
|
|
@ -251,7 +251,7 @@ pure function math_eye(d)
|
|||
math_eye = 0.0_pReal
|
||||
do i=1,d
|
||||
math_eye(i,i) = 1.0_pReal
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function math_eye
|
||||
|
||||
|
|
@ -270,7 +270,7 @@ pure function math_identity4th()
|
|||
#ifndef __INTEL_COMPILER
|
||||
do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
|
||||
math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
|
||||
enddo
|
||||
end do
|
||||
#else
|
||||
forall(i=1:3, j=1:3, k=1:3, l=1:3) &
|
||||
math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
|
||||
|
|
@ -298,7 +298,7 @@ real(pReal) pure function math_LeviCivita(i,j,k)
|
|||
math_LeviCivita = -1.0_pReal
|
||||
else
|
||||
math_LeviCivita = 0.0_pReal
|
||||
endif
|
||||
end if
|
||||
|
||||
end function math_LeviCivita
|
||||
|
||||
|
|
@ -348,7 +348,7 @@ pure function math_outer(A,B)
|
|||
#ifndef __INTEL_COMPILER
|
||||
do concurrent(i=1:size(A,1), j=1:size(B,1))
|
||||
math_outer(i,j) = A(i)*B(j)
|
||||
enddo
|
||||
end do
|
||||
#else
|
||||
forall(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j)
|
||||
#endif
|
||||
|
|
@ -398,7 +398,7 @@ pure function math_mul3333xx33(A,B)
|
|||
#ifndef __INTEL_COMPILER
|
||||
do concurrent(i=1:3, j=1:3)
|
||||
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
|
||||
enddo
|
||||
end do
|
||||
#else
|
||||
forall (i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
|
||||
#endif
|
||||
|
|
@ -421,7 +421,7 @@ pure function math_mul3333xx3333(A,B)
|
|||
#ifndef __INTEL_COMPILER
|
||||
do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
|
||||
math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
|
||||
enddo
|
||||
end do
|
||||
#else
|
||||
forall(i=1:3, j=1:3, k=1:3, l=1:3) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
|
||||
#endif
|
||||
|
|
@ -446,7 +446,7 @@ pure function math_exp33(A,n)
|
|||
n_ = n
|
||||
else
|
||||
n_ = 5
|
||||
endif
|
||||
end if
|
||||
|
||||
invFac = 1.0_pReal ! 0!
|
||||
B = math_I3
|
||||
|
|
@ -456,7 +456,7 @@ pure function math_exp33(A,n)
|
|||
invFac = invFac/real(i,pReal) ! invfac = 1/(i!)
|
||||
B = matmul(B,A)
|
||||
math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/(i!)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function math_exp33
|
||||
|
||||
|
|
@ -514,7 +514,7 @@ pure subroutine math_invert33(InvA, DetA, error, A)
|
|||
|
||||
InvA = InvA/DetA
|
||||
error = .false.
|
||||
endif
|
||||
end if
|
||||
|
||||
end subroutine math_invert33
|
||||
|
||||
|
|
@ -541,7 +541,7 @@ pure function math_invSym3333(A)
|
|||
error stop 'matrix inversion error'
|
||||
else
|
||||
math_invSym3333 = math_66toSym3333(temp66)
|
||||
endif
|
||||
end if
|
||||
|
||||
end function math_invSym3333
|
||||
|
||||
|
|
@ -696,7 +696,7 @@ pure function math_9to33(v9)
|
|||
|
||||
do i = 1, 9
|
||||
math_9to33(MAPPLAIN(1,i),MAPPLAIN(2,i)) = v9(i)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function math_9to33
|
||||
|
||||
|
|
@ -721,7 +721,7 @@ pure function math_sym33to6(m33,weighted)
|
|||
w = merge(NRMMANDEL,1.0_pReal,weighted)
|
||||
else
|
||||
w = NRMMANDEL
|
||||
endif
|
||||
end if
|
||||
|
||||
math_sym33to6 = [(w(i)*m33(MAPNYE(1,i),MAPNYE(2,i)),i=1,6)]
|
||||
|
||||
|
|
@ -748,12 +748,12 @@ pure function math_6toSym33(v6,weighted)
|
|||
w = merge(INVNRMMANDEL,1.0_pReal,weighted)
|
||||
else
|
||||
w = INVNRMMANDEL
|
||||
endif
|
||||
end if
|
||||
|
||||
do i=1,6
|
||||
math_6toSym33(MAPNYE(1,i),MAPNYE(2,i)) = w(i)*v6(i)
|
||||
math_6toSym33(MAPNYE(2,i),MAPNYE(1,i)) = w(i)*v6(i)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function math_6toSym33
|
||||
|
||||
|
|
@ -772,7 +772,7 @@ pure function math_3333to99(m3333)
|
|||
#ifndef __INTEL_COMPILER
|
||||
do concurrent(i=1:9, j=1:9)
|
||||
math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
|
||||
enddo
|
||||
end do
|
||||
#else
|
||||
forall(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
|
||||
#endif
|
||||
|
|
@ -794,7 +794,7 @@ pure function math_99to3333(m99)
|
|||
#ifndef __INTEL_COMPILER
|
||||
do concurrent(i=1:9, j=1:9)
|
||||
math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
|
||||
enddo
|
||||
end do
|
||||
#else
|
||||
forall(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
|
||||
#endif
|
||||
|
|
@ -827,7 +827,7 @@ pure function math_sym3333to66(m3333,weighted)
|
|||
#ifndef __INTEL_COMPILER
|
||||
do concurrent(i=1:6, j=1:6)
|
||||
math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
|
||||
enddo
|
||||
end do
|
||||
#else
|
||||
forall(i=1:6, j=1:6) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
|
||||
#endif
|
||||
|
|
@ -1080,8 +1080,8 @@ pure subroutine math_eigh33(w,v,m)
|
|||
else fallback2
|
||||
v(1:3,2) = v(1:3, 2) / norm
|
||||
v(1:3,3) = math_cross(v(1:3,1),v(1:3,2))
|
||||
endif fallback2
|
||||
endif fallback1
|
||||
end if fallback2
|
||||
end if fallback1
|
||||
|
||||
end subroutine math_eigh33
|
||||
|
||||
|
|
@ -1120,7 +1120,7 @@ pure function math_rotationalPart(F) result(R)
|
|||
lambda = sqrt(math_clip(lambda,0.0_pReal)/3.0_pReal)
|
||||
else
|
||||
lambda = sqrt(I_C(1)/3.0_pReal)
|
||||
endif
|
||||
end if
|
||||
|
||||
I_U = [sum(lambda), lambda(1)*lambda(2)+lambda(2)*lambda(3)+lambda(3)*lambda(1), product(lambda)]
|
||||
|
||||
|
|
@ -1188,7 +1188,7 @@ pure function math_eigvalsh33(m)
|
|||
cos((phi+2.0_pReal*TAU)/3.0_pReal) &
|
||||
] &
|
||||
+ I(1)/3.0_pReal
|
||||
endif
|
||||
end if
|
||||
|
||||
end function math_eigvalsh33
|
||||
|
||||
|
|
@ -1238,7 +1238,7 @@ integer pure function math_binomial(n,k)
|
|||
do i = 1, k_
|
||||
math_binomial = (math_binomial * n_)/i
|
||||
n_ = n_ -1
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function math_binomial
|
||||
|
||||
|
|
@ -1302,7 +1302,7 @@ real(pReal) pure elemental function math_clip(a, left, right)
|
|||
if (present(right)) math_clip = min(right,math_clip)
|
||||
if (present(left) .and. present(right)) then
|
||||
if (left>right) error stop 'left > right'
|
||||
endif
|
||||
end if
|
||||
|
||||
end function math_clip
|
||||
|
||||
|
|
@ -1402,7 +1402,7 @@ subroutine selfTest()
|
|||
|
||||
do while(abs(math_det33(t33))<1.0e-9_pReal)
|
||||
call random_number(t33)
|
||||
enddo
|
||||
end do
|
||||
if (any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) &
|
||||
error stop 'math_inv33'
|
||||
|
||||
|
|
@ -1418,7 +1418,7 @@ subroutine selfTest()
|
|||
|
||||
do while(math_det33(t33)<1.0e-2_pReal) ! O(det(F)) = 1
|
||||
call random_number(t33)
|
||||
enddo
|
||||
end do
|
||||
t33_2 = math_rotationalPart(transpose(t33))
|
||||
t33 = math_rotationalPart(t33)
|
||||
if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) &
|
||||
|
|
|
|||
|
|
@ -97,7 +97,7 @@ subroutine parallelization_init
|
|||
open(OUTPUT_UNIT,file='/dev/null',status='replace') ! close() alone will leave some temp files in cwd
|
||||
else
|
||||
open(OUTPUT_UNIT,encoding='UTF-8') ! for special characters in output
|
||||
endif
|
||||
end if
|
||||
#endif
|
||||
|
||||
print'(/,1x,a)', '<<<+- parallelization init -+>>>'
|
||||
|
|
@ -142,8 +142,8 @@ subroutine parallelization_init
|
|||
!$ if (OMP_NUM_THREADS < 1_pI32) then
|
||||
!$ print'(1x,a)', 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default'
|
||||
!$ OMP_NUM_THREADS = 4_pI32
|
||||
!$ endif
|
||||
!$ endif
|
||||
!$ end if
|
||||
!$ end if
|
||||
!$ print'(1x,a,i0)', 'OMP_NUM_THREADS: ',OMP_NUM_THREADS
|
||||
!$ call omp_set_num_threads(OMP_NUM_THREADS)
|
||||
|
||||
|
|
|
|||
|
|
@ -255,7 +255,7 @@ module subroutine mechanical_init(phases)
|
|||
#else
|
||||
output_mechanical(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
|
||||
#endif
|
||||
enddo
|
||||
end do
|
||||
|
||||
do ce = 1, size(material_phaseID,2)
|
||||
ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
|
||||
|
|
@ -267,14 +267,14 @@ module subroutine mechanical_init(phases)
|
|||
phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(material_V_e_0(ma)%data(1:3,1:3,co), &
|
||||
transpose(phase_mechanical_Fp(ph)%data(1:3,1:3,en)))
|
||||
phase_mechanical_Fi(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%rotate(math_inv33(material_V_e_0(ma)%data(1:3,1:3,co)))
|
||||
enddo
|
||||
enddo
|
||||
end do
|
||||
end do
|
||||
|
||||
do ph = 1, phases%length
|
||||
phase_mechanical_F0(ph)%data = phase_mechanical_F(ph)%data
|
||||
phase_mechanical_Fp0(ph)%data = phase_mechanical_Fp(ph)%data
|
||||
phase_mechanical_Fi0(ph)%data = phase_mechanical_Fi(ph)%data
|
||||
enddo
|
||||
end do
|
||||
|
||||
|
||||
call elastic_init(phases)
|
||||
|
|
@ -284,7 +284,7 @@ module subroutine mechanical_init(phases)
|
|||
call plastic_init()
|
||||
do ph = 1,phases%length
|
||||
plasticState(ph)%state0 = plasticState(ph)%state
|
||||
enddo
|
||||
end do
|
||||
|
||||
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
|
||||
|
||||
|
|
@ -473,25 +473,25 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
|
|||
Lpguess = Lpguess_old &
|
||||
+ deltaLp * stepLengthLp
|
||||
cycle LpLoop
|
||||
endif
|
||||
end if
|
||||
|
||||
calculateJacobiLp: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then
|
||||
jacoCounterLp = jacoCounterLp + 1
|
||||
|
||||
do o=1,3; do p=1,3
|
||||
dFe_dLp(o,1:3,p,1:3) = - Delta_t * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j)
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
dRLp_dLp = math_eye(9) &
|
||||
- math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp))
|
||||
temp_9 = math_33to9(residuumLp)
|
||||
call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
|
||||
if (ierr /= 0) return ! error
|
||||
deltaLp = - math_9to33(temp_9)
|
||||
endif calculateJacobiLp
|
||||
end if calculateJacobiLp
|
||||
|
||||
Lpguess = Lpguess &
|
||||
+ deltaLp * steplengthLp
|
||||
enddo LpLoop
|
||||
end do LpLoop
|
||||
|
||||
call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
|
||||
S, Fi_new, ph,en)
|
||||
|
|
@ -513,7 +513,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
|
|||
Liguess = Liguess_old &
|
||||
+ deltaLi * steplengthLi
|
||||
cycle LiLoop
|
||||
endif
|
||||
end if
|
||||
|
||||
calculateJacobiLi: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then
|
||||
jacoCounterLi = jacoCounterLi + 1
|
||||
|
|
@ -522,10 +522,10 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
|
|||
do o=1,3; do p=1,3
|
||||
dFe_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j)
|
||||
dFi_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*invFi_current
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
do o=1,3; do p=1,3
|
||||
dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new)
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
dRLi_dLi = math_eye(9) &
|
||||
- math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) &
|
||||
+ math_mul3333xx3333(dS_dFi, dFi_dLi))) &
|
||||
|
|
@ -534,11 +534,11 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
|
|||
call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
|
||||
if (ierr /= 0) return ! error
|
||||
deltaLi = - math_9to33(temp_9)
|
||||
endif calculateJacobiLi
|
||||
end if calculateJacobiLi
|
||||
|
||||
Liguess = Liguess &
|
||||
+ deltaLi * steplengthLi
|
||||
enddo LiLoop
|
||||
end do LiLoop
|
||||
|
||||
invFp_new = matmul(invFp_current,B)
|
||||
call math_invert33(Fp_new,devNull,error,invFp_new)
|
||||
|
|
@ -613,9 +613,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
|
|||
if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then
|
||||
broken = plastic_deltaState(ph,en)
|
||||
exit iteration
|
||||
endif
|
||||
end if
|
||||
|
||||
enddo iteration
|
||||
end do iteration
|
||||
|
||||
|
||||
contains
|
||||
|
|
@ -638,7 +638,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
|
|||
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
|
||||
else
|
||||
damper = 1.0_pReal
|
||||
endif
|
||||
end if
|
||||
|
||||
end function damper
|
||||
|
||||
|
|
@ -844,7 +844,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
|
|||
#else
|
||||
dotState = IEEE_FMA(A(n,stage),plastic_RKdotState(1:sizeDotState,n),dotState)
|
||||
#endif
|
||||
enddo
|
||||
end do
|
||||
|
||||
#ifndef __INTEL_LLVM_COMPILER
|
||||
plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t
|
||||
|
|
@ -858,7 +858,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
|
|||
dotState = plastic_dotState(Delta_t*C(stage), ph,en)
|
||||
if (any(IEEE_is_NaN(dotState))) exit
|
||||
|
||||
enddo
|
||||
end do
|
||||
if(broken) return
|
||||
|
||||
|
||||
|
|
@ -950,7 +950,7 @@ subroutine results(group,ph)
|
|||
|
||||
do i = 1, size(dataset,1)
|
||||
to_quaternion(:,i) = dataset(i)%asQuaternion()
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function to_quaternion
|
||||
|
||||
|
|
@ -974,7 +974,7 @@ module subroutine mechanical_forward()
|
|||
phase_mechanical_Lp0(ph) = phase_mechanical_Lp(ph)
|
||||
phase_mechanical_S0(ph) = phase_mechanical_S(ph)
|
||||
plasticState(ph)%state0 = plasticState(ph)%state
|
||||
enddo
|
||||
end do
|
||||
|
||||
end subroutine mechanical_forward
|
||||
|
||||
|
|
@ -1037,7 +1037,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
|
|||
subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en)
|
||||
subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en)
|
||||
subState0 = plasticState(ph)%state(:,en)
|
||||
endif
|
||||
end if
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! cut back (reduced time and restore)
|
||||
else
|
||||
|
|
@ -1048,10 +1048,10 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
|
|||
if (subStep < 1.0_pReal) then ! actual (not initial) cutback
|
||||
phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0
|
||||
phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0
|
||||
endif
|
||||
end if
|
||||
plasticState(ph)%state(:,en) = subState0
|
||||
todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair)
|
||||
endif
|
||||
end if
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! prepare for integration
|
||||
|
|
@ -1060,9 +1060,9 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
|
|||
subF = subF0 &
|
||||
+ subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
|
||||
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en)
|
||||
endif
|
||||
end if
|
||||
|
||||
enddo cutbackLooping
|
||||
end do cutbackLooping
|
||||
|
||||
end function phase_mechanical_constitutive
|
||||
|
||||
|
|
@ -1086,14 +1086,14 @@ module subroutine mechanical_restore(ce,includeL)
|
|||
if (includeL) then
|
||||
phase_mechanical_Lp(ph)%data(1:3,1:3,en) = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
|
||||
phase_mechanical_Li(ph)%data(1:3,1:3,en) = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
|
||||
endif ! maybe protecting everything from overwriting makes more sense
|
||||
end if ! maybe protecting everything from overwriting makes more sense
|
||||
|
||||
phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
|
||||
phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
|
||||
phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)
|
||||
|
||||
plasticState(ph)%state(:,en) = plasticState(ph)%State0(:,en)
|
||||
enddo
|
||||
end do
|
||||
|
||||
end subroutine mechanical_restore
|
||||
|
||||
|
|
@ -1164,7 +1164,7 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
|
|||
lhs_3333(1:3,o,1:3,p) = IEEE_FMA(invFi,invFi(p,o),lhs_3333(1:3,o,1:3,p))
|
||||
rhs_3333(1:3,1:3,o,p) = IEEE_FMA(matmul(invSubFi0,dLidS(1:3,1:3,o,p)),-Delta_t,rhs_3333(1:3,1:3,o,p))
|
||||
#endif
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
call math_invert(temp_99,error,math_3333to99(lhs_3333))
|
||||
if (error) then
|
||||
call IO_warning(600,'inversion error in analytic tangent calculation', &
|
||||
|
|
@ -1172,9 +1172,9 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
|
|||
dFidS = 0.0_pReal
|
||||
else
|
||||
dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
|
||||
endif
|
||||
end if
|
||||
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
|
||||
endif
|
||||
end if
|
||||
|
||||
call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
|
||||
phase_mechanical_S(ph)%data(1:3,1:3,en), &
|
||||
|
|
@ -1191,7 +1191,7 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
|
|||
rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
|
||||
temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
|
||||
+ matmul(temp_33_3,dLidS(1:3,1:3,p,o))
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
#ifndef __INTEL_LLVM_COMPILER
|
||||
lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t &
|
||||
+ math_mul3333xx3333(dSdFi,dFidS)
|
||||
|
|
@ -1206,14 +1206,14 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
|
|||
dSdF = rhs_3333
|
||||
else
|
||||
dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
|
||||
endif
|
||||
end if
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! calculate dFpinvdF
|
||||
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
|
||||
do o=1,3; do p=1,3
|
||||
dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * Delta_t
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! assemble dPdF
|
||||
|
|
@ -1224,13 +1224,13 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
|
|||
dPdF = 0.0_pReal
|
||||
do p=1,3
|
||||
dPdF(p,1:3,p,1:3) = transpose(matmul(invFp,temp_33_1))
|
||||
enddo
|
||||
end do
|
||||
do o=1,3; do p=1,3
|
||||
dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
|
||||
+ matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
|
||||
+ matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) &
|
||||
+ matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
|
||||
enddo; enddo
|
||||
end do; end do
|
||||
|
||||
end function phase_mechanical_dPdF
|
||||
|
||||
|
|
|
|||
|
|
@ -61,7 +61,7 @@ module subroutine eigen_init(phases)
|
|||
|
||||
if (maxval(Nmodels) /= 0) then
|
||||
where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID
|
||||
endif
|
||||
end if
|
||||
|
||||
allocate(model_damage(phases%length), source = EIGEN_UNDEFINED_ID)
|
||||
|
||||
|
|
|
|||
|
|
@ -435,7 +435,7 @@ function plastic_active(plastic_label) result(active_plastic)
|
|||
mech => phase%get('mechanical')
|
||||
pl => mech%get('plastic',defaultVal = emptyDict)
|
||||
active_plastic(ph) = pl%get_asString('type',defaultVal='none') == plastic_label
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function plastic_active
|
||||
|
||||
|
|
|
|||
|
|
@ -364,13 +364,13 @@ module function plastic_dislotwin_init() result(myPlasticity)
|
|||
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), &
|
||||
phase_lattice(ph))
|
||||
if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
|
||||
endif slipAndTwinActive
|
||||
end if slipAndTwinActive
|
||||
|
||||
slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
|
||||
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dFloat('h_sl-tr'), &
|
||||
phase_lattice(ph))
|
||||
if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
|
||||
endif slipAndTransActive
|
||||
end if slipAndTransActive
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! allocate state arrays
|
||||
|
|
|
|||
|
|
@ -122,25 +122,25 @@ module subroutine thermal_init(phases)
|
|||
|
||||
allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
|
||||
|
||||
enddo
|
||||
end do
|
||||
|
||||
allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
|
||||
|
||||
if (maxval(thermal_Nsources) /= 0) then
|
||||
where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
|
||||
where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
|
||||
endif
|
||||
end if
|
||||
|
||||
thermal_source_maxSizeDotState = 0
|
||||
do ph = 1,phases%length
|
||||
|
||||
do so = 1,thermal_Nsources(ph)
|
||||
thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0
|
||||
enddo
|
||||
end do
|
||||
|
||||
thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, &
|
||||
maxval(thermalState(ph)%p%sizeDotState))
|
||||
enddo
|
||||
end do
|
||||
|
||||
end subroutine thermal_init
|
||||
|
||||
|
|
@ -170,7 +170,7 @@ module function phase_f_T(ph,en) result(f)
|
|||
|
||||
end select
|
||||
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function phase_f_T
|
||||
|
||||
|
|
@ -195,7 +195,7 @@ function phase_thermal_collectDotState(ph,en) result(broken)
|
|||
|
||||
broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))
|
||||
|
||||
enddo SourceLoop
|
||||
end do SourceLoop
|
||||
|
||||
end function phase_thermal_collectDotState
|
||||
|
||||
|
|
@ -262,7 +262,7 @@ function integrateThermalState(Delta_t, ph,en) result(broken)
|
|||
sizeDotState = thermalState(ph)%p(so)%sizeDotState
|
||||
thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) &
|
||||
+ thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function integrateThermalState
|
||||
|
||||
|
|
@ -277,7 +277,7 @@ module subroutine thermal_restartWrite(groupHandle,ph)
|
|||
|
||||
do so = 1,thermal_Nsources(ph)
|
||||
call HDF5_write(thermalState(ph)%p(so)%state,groupHandle,'omega_thermal')
|
||||
enddo
|
||||
end do
|
||||
|
||||
end subroutine thermal_restartWrite
|
||||
|
||||
|
|
@ -292,7 +292,7 @@ module subroutine thermal_restartRead(groupHandle,ph)
|
|||
|
||||
do so = 1,thermal_Nsources(ph)
|
||||
call HDF5_read(thermalState(ph)%p(so)%state0,groupHandle,'omega_thermal')
|
||||
enddo
|
||||
end do
|
||||
|
||||
end subroutine thermal_restartRead
|
||||
|
||||
|
|
@ -305,8 +305,8 @@ module subroutine thermal_forward()
|
|||
do ph = 1, size(thermalState)
|
||||
do so = 1, size(thermalState(ph)%p)
|
||||
thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state
|
||||
enddo
|
||||
enddo
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine thermal_forward
|
||||
|
||||
|
|
@ -380,8 +380,8 @@ function thermal_active(source_label,src_length) result(active_source)
|
|||
do s = 1, sources%length
|
||||
src => sources%get(s)
|
||||
active_source(s,p) = src%get_asString('type') == source_label
|
||||
enddo
|
||||
enddo
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
end function thermal_active
|
||||
|
|
|
|||
|
|
@ -112,7 +112,7 @@ pure function eval(self,x) result(y)
|
|||
#else
|
||||
y = IEEE_FMA(y,x-self%x_ref,self%coef(o))
|
||||
#endif
|
||||
enddo
|
||||
end do
|
||||
|
||||
end function eval
|
||||
|
||||
|
|
|
|||
|
|
@ -123,7 +123,7 @@ function getCWD()
|
|||
getCWD = c_f_string(getCWD_Cstring)
|
||||
else
|
||||
error stop 'invalid working directory'
|
||||
endif
|
||||
end if
|
||||
|
||||
end function getCWD
|
||||
|
||||
|
|
@ -145,7 +145,7 @@ function getHostName()
|
|||
getHostName = c_f_string(getHostName_Cstring)
|
||||
else
|
||||
getHostName = 'n/a (Error!)'
|
||||
endif
|
||||
end if
|
||||
|
||||
end function getHostName
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue