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

regridding is now working, changed the subroutine into a function 

changed order of arrays in nearest neighbor search to make it fortran fast
constitutive.f90 and homogenization.f90 write state size out during initialization
setup/setup_processing.py is using byterecl to be compatible with binary files written out by solver
This commit is contained in:
Krishna Komerla 2012-05-08 14:57:06 +00:00
parent a9a72cee97
commit c752dd5474
7 changed files with 194 additions and 128 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
code/DAMASK_spectral_interface.f90
View File

@ -183,12 +183,12 @@ subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn)
call GET_ENVIRONMENT_VARIABLE('HOST',hostName)
call GET_ENVIRONMENT_VARIABLE('USER',userName)
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
write(6,'(a,a)') 'Host Name: ', trim(hostName)
write(6,'(a,a)') 'User Name: ', trim(userName)
write(6,'(a,a)') 'Path Separator: ', getPathSep()

8
code/constitutive.f90
View File

@ -29,6 +29,7 @@
MODULE constitutive
use prec, only: pInt, p_vec
use IO, only: IO_write_jobBinaryFile
implicit none
type(p_vec), dimension(:,:,:), allocatable :: &
@ -356,6 +357,13 @@ endif
enddo
!$OMP END PARALLEL DO
!----- write out state size file----------------
open(777)
call IO_write_jobBinaryFile(777,'sizeStateConst', size(constitutive_sizeState))
write (777,rec=1) constitutive_sizeState
close(777)
!-----------------------------------------------
constitutive_maxSizeState = maxval(constitutive_sizeState)
constitutive_maxSizeDotState = maxval(constitutive_sizeDotState)
constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)

11
code/damask.core.pyf
View File

@ -167,8 +167,8 @@ python module core ! in
real*8, dimension(3), intent(in) :: geomdim
integer, intent(in) :: queryPoints
integer, intent(in) :: domainPoints
real*8, dimension(queryPoints,spatialDim), intent(in), depend(queryPoints,spatialDim) :: querySet
real*8, dimension(domainPoints,spatialDim), intent(in), depend(domainPoints,spatialDim) :: domainSet
real*8, dimension(spatialDim,queryPoints), intent(in), depend(queryPoints,spatialDim) :: querySet
real*8, dimension(spatialDim,domainPoints), intent(in), depend(domainPoints,spatialDim) :: domainSet
! output variable
integer, dimension(queryPoints), intent(out), depend(queryPoints) :: indices
! depending on input
@ -178,9 +178,10 @@ python module core ! in
module mesh ! in :mesh:mesh.f90
subroutine mesh_regrid(resNew) ! in :mesh:mesh.f90
integer, dimension(3), intent(in,out), optional :: resNew
end subroutine mesh_regrid
function mesh_regrid(resNewInput) ! in :mesh:mesh.f90
integer, dimension(3) :: mesh_regrid
integer, dimension(3), intent(in), optional :: resNewInput
end function mesh_regrid
end module mesh

11
code/homogenization.f90
View File

@ -34,6 +34,7 @@ MODULE homogenization
!*** Include other modules ***
use prec, only: pInt,pReal,p_vec
use IO, only: IO_write_jobBinaryFile
implicit none
! ****************************************************************
@ -161,7 +162,7 @@ forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems)
end forall
! --- ALLOCATE AND INITIALIZE GLOBAL STATE AND POSTRESULTS VARIABLES ---
! --- ALLOCATE AND INITIALIZE GLOBAL STATE AND POSTRESULTS VARIABLES ----------
!$OMP PARALLEL DO PRIVATE(myInstance)
do e = 1,mesh_NcpElems ! loop over elements
@ -192,6 +193,14 @@ end forall
enddo
enddo
!$OMP END PARALLEL DO
!---write state size file out---------------------------------------
open(777)
call IO_write_jobBinaryFile(777,'sizeStateHomog',size(homogenization_sizeState))
write (777,rec=1) homogenization_sizeState
close(777)
!--------------------------------------------------------------
homogenization_maxSizeState = maxval(homogenization_sizeState)
homogenization_maxSizePostResults = maxval(homogenization_sizePostResults)
materialpoint_sizeResults = 1 & ! grain count

16
code/math.f90
View File

@ -3832,19 +3832,19 @@ subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, do
real(pReal), dimension(3), intent(in) :: geomdim
integer(pInt), intent(in) :: domainPoints
integer(pInt), intent(in) :: queryPoints
real(pReal), dimension(queryPoints,spatialDim), intent(in) :: querySet
real(pReal), dimension(domainPoints,spatialDim), intent(in) :: domainSet
real(pReal), dimension(spatialDim,queryPoints), intent(in) :: querySet
real(pReal), dimension(spatialDim,domainPoints), intent(in) :: domainSet
! output variable
integer(pInt), dimension(queryPoints), intent(out) :: indices
! other variables depending on input
real(pReal), dimension(3,(3_pInt**spatialDim)*domainPoints) :: domainSetLarge
real(pReal), dimension(spatialDim,(3_pInt**spatialDim)*domainPoints) :: domainSetLarge
! other variables
integer(pInt) :: i,j, l,m,n
type(kdtree2), pointer :: tree
type(kdtree2_result), dimension(1) :: Results
if (size(querySet(1,:)) /= spatialDim) call IO_error(407_pInt,ext_msg='query set')
if (size(domainSet(1,:)) /= spatialDim) call IO_error(407_pInt,ext_msg='domain set')
if (size(querySet(:,1)) /= spatialDim) call IO_error(407_pInt,ext_msg='query set')
if (size(domainSet(:,1)) /= spatialDim) call IO_error(407_pInt,ext_msg='domain set')
i = 0_pInt
@ -3852,14 +3852,14 @@ subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, do
do j = 1_pInt, domainPoints
do l = -1_pInt, 1_pInt; do m = -1_pInt, 1_pInt
i = i + 1_pInt
domainSetLarge(1:3,i) = domainSet(j,1:3) + math_mul33x3(Favg,real([l,m,0_pInt],pReal)*geomdim)
domainSetLarge(1:2,i) = domainSet(1:2,j) +matmul(Favg(1:2,1:2),real([l,m],pReal)*geomdim(1:2))
enddo; enddo
enddo
else
do j = 1_pInt, domainPoints
do l = -1_pInt, 1_pInt; do m = -1_pInt, 1_pInt; do n = -1_pInt, 1_pInt
i = i + 1_pInt
domainSetLarge(1:3,i) = domainSet(j,1:3) + math_mul33x3(Favg,real([l,m,n],pReal)*geomdim)
domainSetLarge(1:3,i) = domainSet(1:3,j) + math_mul33x3(Favg,real([l,m,n],pReal)*geomdim)
enddo; enddo; enddo
enddo
endif
@ -3867,7 +3867,7 @@ subroutine math_nearestNeighborSearch(spatialDim, Favg, geomdim, queryPoints, do
tree => kdtree2_create(domainSetLarge,sort=.true.,rearrange=.true.)
do j = 1_pInt, queryPoints
call kdtree2_n_nearest(tp=tree, qv=querySet(j,1:spatialDim),nn=1_pInt, results = Results)
call kdtree2_n_nearest(tp=tree, qv=querySet(1:spatialDim,j),nn=1_pInt, results = Results)
indices(j) = Results(1)%idx
enddo

258
code/mesh.f90
View File

@ -3528,7 +3528,7 @@ end subroutine mesh_tell_statistics
subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automatic determination of new grid
function mesh_regrid(resNewInput) !use new_res=[0.0,0.0,0.0] for automatic determination of new grid
use prec, only: &
pInt, &
pReal
@ -3542,86 +3542,122 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati
deformed_FFT, &
math_mul33x3
integer(pInt), dimension(3), optional, intent(inout) :: resNew
integer(pInt):: maxsize, i, j, k, ielem, Npoints, NpointsNew
integer(pInt), dimension(3) :: res
integer(pInt), dimension(:), allocatable :: indices, outputSize
real(pReal), dimension(3) :: geomdim = 0.0_pReal
integer(pInt), dimension(3) :: mesh_regrid
integer(pInt), dimension(3), optional, intent(in) :: resNewInput
integer(pInt):: maxsize, i, j, k, ielem, Npoints, NpointsNew, spatialDim
integer(pInt), dimension(3) :: res, resNew
integer(pInt), dimension(:), allocatable :: indices
real(pReal), dimension(3) :: geomdim
real(pReal), dimension(3,3) :: Favg
real(pReal), dimension(:,:), allocatable :: &
coordinatesNew
coordinatesNew, &
coordinatesLinear
real(pReal), dimension(:,:,:), allocatable :: &
material_phase, material_phaseNew
real(pReal), dimension(:,:,:,:,:), allocatable :: &
F, FNew, &
Fp, FpNew, &
Lp, LpNew, &
dcsdE, dcsdENew
real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: &
dPdF, dPdFNew
real(pReal), dimension (:,:,:,:), allocatable :: &
coordinates, &
Tstar, TstarNew, &
stateConst, stateConstNew, &
stateHomog, stateHomogNew
integer(pInt), dimension (:,:), allocatable :: SizeConst
coordinates, &
Tstar, TstarNew, &
stateHomog, &
stateConst
integer(pInt), dimension(:,:), allocatable :: &
sizeStateConst, sizeStateHomog
res = mesh_spectral_getResolution()
geomdim = mesh_spectral_getDimension()
if (.not. present(resNew)) then
resNew=res
endif
Npoints = res(1)*res(2)*res(3)
NpointsNew = resNew(1)*resNew(2)*resNew(3)
print*, 'resolution ', res
print*, 'new resolution', resNew
allocate(coordinates(res(1),res(2),res(3),3))
allocate(coordinatesNew(resNew(1)*resNew(2)*resNew(3),3))
allocate(indices(Npoints))
allocate(F(res(1),res(2),res(3),3,3))
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
! ----Calculate average deformation for remesh--------
! ----Calculate deformed configuration and average--------
do i= 1_pInt,3_pInt; do j = 1_pInt,3_pInt
Favg(i,j) = sum(F(1:res(1),1:res(2),1:res(3),i,j)) / Npoints
enddo; enddo
allocate(coordinates(res(1),res(2),res(3),3))
call deformed_fft(res,geomdim,Favg,1.0_pReal,F,coordinates)
deallocate(F)
! ----Store coordinates into a linear list--------
allocate(coordinatesLinear(3,Npoints))
ielem = 0_pInt
do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1)
ielem = ielem + 1_pInt
coordinatesLinear(1:3,ielem) = coordinates(i,j,k,1:3)
enddo; enddo; enddo
deallocate(coordinates)
! ----For 2D /3D case----------------------------------
if (res(3)== 1_pInt) then
spatialDim = 2_pInt
else
spatialDim = 3_pInt
endif
! ----determine/calculate new resolution--------------
if (.not. present(resNewInput)) then
resNew=res
elseif(all(resNewInput==0.0_pReal)) then
!do automatic determination
else
resNew=resNewInput
endif
NpointsNew = resNew(1)*resNew(2)*resNew(3)
! ----Calculate regular new coordinates------------
allocate(coordinatesNew(3,NpointsNew))
ielem = 0_pInt
do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1)
ielem = ielem + 1_pInt
coordinatesNew(ielem,1:3) = math_mul33x3(Favg, geomdim/real(resNew,pReal)*real([i,j,k],pReal) &
coordinatesNew(1:3,ielem) = math_mul33x3(Favg, geomdim/real(resNew,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*resNew,pReal))
enddo; enddo; enddo
!----- Nearest neighbour search -----------------------
call math_nearestNeighborSearch(res, Favg, geomdim, NpointsNew, &
coordinates, coordinatesNew, indices)
indices = indices /3_pInt**3_pInt +1_pInt ! adjust it for 2D case
allocate(indices(Npoints))
call math_nearestNeighborSearch(spatialDim, Favg, geomdim, NpointsNew, Npoints, &
coordinatesNew, coordinatesLinear, indices)
deallocate(coordinatesNew)
indices = indices / 3_pInt**spatialDim +1_pInt
deallocate(F)
deallocate(coordinates)
deallocate(coordinatesNew)
!---------------------------------------------------------------------------
allocate(material_phase (1,1, Npoints))
allocate(material_phaseNew (1,1, NpointsNew))
call IO_read_jobBinaryFile(777,'recordedPhase',trim(getSolverJobName()),size(material_phase))
read (777,rec=1) material_phase
close (777)
do i = 1, NpointsNew
material_phaseNew(1,1,i) = material_phase(1,1,indices(i))
enddo
! ---------------------------------------------------------------------------
allocate(F (3,3,1,1, Npoints))
allocate(FNew(3,3,1,1, NpointsNew))
call IO_write_jobBinaryFile(777,'recordedPhase',size(material_phaseNew))
write (777,rec=1) material_phaseNew
close (777)
deallocate(material_phase)
deallocate(material_phaseNew)
!---------------------------------------------------------------------------
allocate(F (3,3,1,1, Npoints))
allocate(FNew (3,3,1,1, NpointsNew))
call IO_read_jobBinaryFile(777,'convergedF',trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
call flush(6)
do i = 1, NpointsNew
FNew(1:3,1:3,1,1,i) = F(1:3,1:3,1,1,indices(i))
enddo
call flush(6)
call IO_write_jobBinaryFile(777,'convergedF',size(FNew))
write (777,rec=1) FNew
@ -3630,15 +3666,14 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati
deallocate(FNew)
!---------------------------------------------------------------------
allocate(Fp (3,3,1,1,Npoints))
allocate(FpNew(3,3,1,1,NpointsNew))
allocate(Fp (3,3,1,1,Npoints))
allocate(FpNew (3,3,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergedFp',trim(getSolverJobName()),size(Fp))
read (777,rec=1) Fp
close (777)
do i = 1, NpointsNew
FpNew(1:3,1:3,1,1,i) = Fp(1:3,1:3,1,1,indices(i))
enddo
enddo
call IO_write_jobBinaryFile(777,'convergedFp',size(FpNew))
write (777,rec=1) FpNew
@ -3647,50 +3682,44 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati
deallocate(FpNew)
!------------------------------------------------------------------------
allocate(Lp (3,3,1,1,Npoints))
allocate(LpNew(3,3,1,1,NpointsNew))
allocate(Lp (3,3,1,1,Npoints))
allocate(LpNew (3,3,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergedLp',trim(getSolverJobName()),size(Lp))
read (777,rec=1) Lp
close (777)
do i = 1, NpointsNew
LpNew(1:3,1:3,1,1,i) = Lp(1:3,1:3,1,1,indices(i))
enddo
call IO_write_jobBinaryFile(777,'convergedLp',size(LpNew))
write (777,rec=1) LpNew
close (777)
deallocate(Lp)
deallocate(LpNew)
!----------------------------------------------------------------------------
allocate(dcsdE (6,6,1,1,Npoints))
allocate(dcsdENew(6,6,1,1,NpointsNew))
!----------------------------------------------------------------------------
allocate(dcsdE (6,6,1,1,Npoints))
allocate(dcsdENew (6,6,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergeddcsdE',trim(getSolverJobName()),size(dcsdE))
read (777,rec=1) dcsdE
close (777)
do i = 1, NpointsNew
dcsdENew(1:6,1:6,1,1,i) = dcsdE(1:6,1:6,1,1,indices(i))
enddo
call IO_write_jobBinaryFile(777,'convergeddcsdE',size(dcsdENew))
write (777,rec=1) dcsdENew
close (777)
deallocate(dcsdE)
deallocate(dcsdENew)
! ---------------------------------------------------------------------------
allocate(dPdF (3,3,3,3,1,1,Npoints))
allocate(dPdFNew(3,3,3,3,1,1,NpointsNew))
!---------------------------------------------------------------------------
allocate(dPdF (3,3,3,3,1,1,Npoints))
allocate(dPdFNew (3,3,3,3,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergeddPdF',trim(getSolverJobName()),size(dPdF))
read (777,rec=1) dPdF
close (777)
do i = 1, NpointsNew
dPdFNew(1:3,1:3,1:3,1:3,1,1,i) = dPdF(1:3,1:3,1:3,1:3,1,1,indices(i))
enddo
call IO_write_jobBinaryFile(777,'convergeddPdF',size(dPdFNew))
write (777,rec=1) dPdFNew
close (777)
@ -3698,66 +3727,85 @@ subroutine mesh_regrid(resNew) !use new_res=[0.0,0.0,0.0] for automati
deallocate(dPdFNew)
!---------------------------------------------------------------------------
allocate(Tstar (6,1,1,Npoints))
allocate(TstarNew(6,1,1,NpointsNew))
allocate(Tstar (6,1,1,Npoints))
allocate(TstarNew (6,1,1,NpointsNew))
call IO_read_jobBinaryFile(777,'convergedTstar',trim(getSolverJobName()),size(Tstar))
read (777,rec=1) Tstar
close (777)
do i = 1, NpointsNew
TstarNew(1:6,1,1,i) = Tstar(1:6,1,1,indices(i))
enddo
call IO_write_jobBinaryFile(777,'convergedTstar',size(TstarNew))
write (777,rec=1) TstarNew
close (777)
deallocate(Tstar)
deallocate(TstarNew)
!----------------------------------------------------------------------------
!for homog and state const
! allocate(SizeConst (1,Npoints))
! allocate(StateConst (1,Npoints))
! allocate(StateConstNew (1,Npoints))
! call IO_read_jobBinaryFile(777,'SizeConst',trim(getSolverJobName()))
! read (777,rec=1) SizeConst
! close(777)
! maxsize = maxval(SizeConst)
! call IO_read_jobBinaryFile(777,'convergedStateConst',trim(getSolverJobName()))
! do i =1, Npoints
! do l = 1,SizeConst(1,i)
! read(777,rec=i) StateConst(1,1,i,l)
! enddo
! enddo
! close(777)
allocate(sizeStateConst(1,Npoints))
call IO_read_jobBinaryFile(777,'sizeStateConst',trim(getSolverJobName()),size(sizeStateConst))
read (777,rec=1) sizeStateConst
close (777)
maxsize = maxval(sizeStateConst)
allocate(StateConst (1,1,Npoints,maxsize))
call IO_read_jobBinaryFile(777,'convergedStateConst',trim(getSolverJobName()))
k = 0_pInt
do i =1, Npoints
do j = 1,sizeStateConst(1,i)
k = k+1_pInt
read(777,rec=k) StateConst(1,1,i,j)
enddo
enddo
close(777)
call IO_write_jobBinaryFile(777,'convergedStateConst')
k = 0_pInt
do i = 1,NpointsNew
do j = 1,sizeStateConst(1,i)
k=k+1_pInt
write(777,rec=k) StateConst(1,1,indices(i),j)
enddo
enddo
close (777)
deallocate(sizeStateConst)
deallocate(StateConst)
!----------------------------------------------------------------------------
allocate(sizeStateHomog(1,Npoints))
call IO_read_jobBinaryFile(777,'sizeStateHomog',trim(getSolverJobName()),size(sizeStateHomog))
read (777,rec=1) sizeStateHomog
close (777)
maxsize = maxval(sizeStateHomog)
allocate(stateHomog (1,1,Npoints,maxsize))
call IO_read_jobBinaryFile(777,'convergedStateHomog',trim(getSolverJobName()))
k = 0_pInt
do i =1, Npoints
do j = 1,sizeStateHomog(1,i)
k = k+1_pInt
read(777,rec=k) stateHomog(1,1,i,j)
enddo
enddo
close(777)
call IO_write_jobBinaryFile(777,'convergedStateHomog')
k = 0_pInt
do i = 1,NpointsNew
do j = 1,sizeStateHomog(1,i)
k=k+1_pInt
write(777,rec=k) stateHomog(1,1,indices(i),j)
enddo
enddo
close (777)
deallocate(sizeStateHomog)
deallocate(stateHomog)
! do i = 1, NpointsNew
! do l = 1,SizeConst(1,i)
! StateConstNew(1,1,i,l) = StateConst(1,1,indices(i),l)
! enddo
! enddo
deallocate(indices)
! call IO_write_jobBinaryFile(777,'convergedStateConst',trim(getSolverJobName()))
! m = 0_pInt
! do i = 1,NpointsNew
! do l = 1,SizeConst(1,i)
! write(777,rec=m) StateConstNew(1,1,i,l)
! enddo
! enddo
! close (777)
!----------------------------------------------------------------------------
! State homog
! call IO_read_jobBinaryFile(777,'convergedStateHomog',FEmodelGeometry)
! m = 0_pInt
! do k = 1,Npoints;
! do l = 1,homogenization_sizeState(j,k)
! m = m+1_pInt
! read(777,rec=m) stateHomog(1,k)%p(l)
! enddo
! enddo; enddo
! close (777)
end subroutine mesh_regrid
mesh_regrid = resNew
end function mesh_regrid
function mesh_spectral_getDimension(fileUnit)
use IO, only: &

6
processing/setup/setup_processing.py
View File

@ -50,9 +50,9 @@ if options.compiler not in compilers:
f2py_compiler = {
'gfortran': 'gnu95 --f90flags="-fno-range-check -xf95-cpp-input -std=f2008 -fall-intrinsics"',
'gnu95': 'gnu95 --f90flags="-fno-range-check -xf95-cpp-input -std=f2008 -fall-intrinsics"',
'intel32': 'intel --f90flags="-fpp -stand f03 -diag-disable 5268"',
'intel': 'intelem --f90flags="-fpp -stand f03 -diag-disable 5268"',
'ifort': 'intelem --f90flags="-fpp -stand f03 -diag-disable 5268"',
'intel32': 'intel --f90flags="-fpp -stand f03 -diag-disable 5268 -assume byterecl"',
'intel': 'intelem --f90flags="-fpp -stand f03 -diag-disable 5268 -assume byterecl"',
'ifort': 'intelem --f90flags="-fpp -stand f03 -diag-disable 5268 -assume byterecl"',
}[options.compiler]
compiler = {
'gfortran': 'gfortran',