Skip to content
Snippets Groups Projects
Commit 586c32f8 authored by Dave Whipp's avatar Dave Whipp
Browse files

Updated vermajor to 3 (v0.3), added options for restarting with compaction...

Updated vermajor to 3 (v0.3), added options for restarting with compaction output, modified code to not read/write output for certain fields when they are not used (isostasy, compaction, nest)
parent fa59150c
No related branches found
No related tags found
No related merge requests found
......@@ -106,7 +106,7 @@ integer input_vermajor,input_verminor
integer output_vermajor,output_verminor,output_verstat,output_verrev
logical,dimension(:),allocatable::influid,do_it,subset_leaves,instrain
logical :: nest
logical :: nest,isostasyout,isobcout,nestout,compactionout
character clsf*3,c4*4,cc4*4,dir*128
character cs,nestin
......@@ -121,7 +121,7 @@ double precision,dimension(:),allocatable::s11,s12,s13,s22,s23,s33,str11,str12,s
double precision,dimension(:),allocatable::azimuth1,azimuth3,dip1,dip3
double precision zmax,dz,l1,l2,l3,n11,n12,n13,n21,n22,n23,n31,n32,n33,con,dxy
double precision :: sselemx,sselemy,sselemz,xminls,yminls,zminls,usum,vsum,wsum
double precision :: uvwsum,disp,wmod
double precision :: uvwsum,disp,wmod,wiso,wcompact
double precision dxyz,x,y,z,xcut,dumpdp
integer icut
......@@ -141,9 +141,9 @@ nne=8
! 3=release
! verrev is incremented for minor changes and bugfixes
vermajor=0
verminor=2
verminor=3
verstat=0
verrev=2
verrev=0
write (*,'(a)') '!------------------------------------------------------------------------------|'
write (*,'(a)') '!------------------------------------------------------------------------------|'
......@@ -328,17 +328,19 @@ else
endif
write(*,'(a)') '--------------------------------------------------------------------------'
read (7) isostasyout,isobcout,nestout,compactionout
read (7) ov%noctree,ov%nnode,ov%nleaves,ner,ov%nlsf,npcl,current_time
nn=ov%nnode
allocate(ov%on(nn))
allocate(ov%x(nn),ov%y(nn),ov%z(nn))
allocate(ov%u(nn),ov%v(nn),ov%w(nn),ov%wiso(nn),ov%wcompact(nn))
allocate(ov%u(nn),ov%v(nn),ov%w(nn))
allocate(ov%pressure(ov%nleaves),ov%strain(nn))
allocate(ov%e2d(ov%nleaves),ov%eviscosity(ov%nleaves),ov%is_plastic(ov%nleaves))
allocate(ov%dilatr(ov%nleaves),ov%whole_leaf_in_fluid(ov%nleaves))
allocate(ov%crit(ov%nleaves),ov%matnum(ov%nleaves),ov%yield_ratio(ov%nleaves))
allocate(ov%frict_angle(ov%nleaves),ov%compaction_density(ov%nleaves))
allocate(ov%frict_angle(ov%nleaves))
allocate(ov%temp(nn))
allocate(ov%lsf(nn,ov%nlsf))
allocate(influid(nn))
......@@ -352,19 +354,24 @@ allocate(cl%x(npcl),cl%y(npcl),cl%z(npcl),cl%x0(npcl),cl%y0(npcl),cl%z0(npcl))
allocate(cl%strain(npcl),cl%lsf0(npcl),cl%temp(npcl),cl%press(npcl))
allocate(cl%e2dp(npcl),cl%tag(npcl),cl%matnum(npcl),cl%ematnump(npcl))
if (isostasyout) allocate(ov%wiso(nn))
if (compactionout) allocate(ov%wcompact(nn),ov%compaction_density(ov%nleaves))
!========================
!=====[nodal values]=====
!========================
read(7)(ov%x(i),ov%y(i),ov%z(i), &
ov%u(i),ov%v(i),ov%w(i), &
ov%wiso(i),ov%wcompact(i), &
ov%lsf(i,1:ov%nlsf), &
ov%temp(i), &
ov%nodal_pressure(i), &
ov%strain(i), &
kfx,kfy,kfz,kft,i=1,nn)
if (isostasyout) read (7) (ov%wiso(i),i=1,nn)
if (compactionout) read (7) (ov%wcompact(i),i=1,nn)
!======================
!=====[icon array]=====
!======================
......@@ -380,10 +387,11 @@ read(7) (ov%icon(1:8,ie), &
ov%matnum(ie), &
ov%yield_ratio(ie), &
ov%frict_angle(ie), &
ov%compaction_density(ie), &
ov%whole_leaf_in_fluid(ie), &
ie=1,ov%nleaves)
if (compactionout) read (7) (ov%compaction_density(ie),ie=1,ov%nleaves)
!==============================
!=====[octree information]=====
!==============================
......@@ -506,16 +514,24 @@ read (7) (cl%x(i), &
!=====[isostasy information]=====
!================================
read (7) nb
dxy=1.d0/real(nb)
allocate(zisodisp(nb+1,nb+1))
read (7) ((zisodisp(i,j),j=1,nb+1),i=1,nb+1)
if (isobcout) then
read (7) nb
dxy=1.d0/real(nb)
allocate(zisodisp(nb+1,nb+1))
read (7) ((zisodisp(i,j),j=1,nb+1),i=1,nb+1)
endif
!============================
!=====[nest information]=====
!============================
read(7) sselemx,sselemy,sselemz,xminls,yminls,zminls
if (nestout) read(7) sselemx,sselemy,sselemz,xminls,yminls,zminls
!======================================
!=====[density string information]=====
!======================================
! Not yet implemented
close(7)
......@@ -576,8 +592,8 @@ write(*,'(a,2f30.20)') 'z range : ',minval(ov%z(1:nn)), ma
write(*,'(a,2f30.20)') 'u range : ',minval(ov%u(1:nn)), maxval(ov%u(1:nn))
write(*,'(a,2f30.20)') 'v range : ',minval(ov%v(1:nn)), maxval(ov%v(1:nn))
write(*,'(a,2f30.20)') 'w range : ',minval(ov%w(1:nn)), maxval(ov%w(1:nn))
write(*,'(a,2f30.20)') 'w (iso) range : ',minval(ov%wiso(1:nn)), maxval(ov%wiso(1:nn))
write(*,'(a,2f30.20)') 'w (compaction) range : ',minval(ov%wcompact(1:nn)),maxval(ov%wcompact(1:nn))
if (isostasyout) write(*,'(a,2f30.20)') 'w (iso) range : ',minval(ov%wiso(1:nn)), maxval(ov%wiso(1:nn))
if (compactionout) write(*,'(a,2f30.20)') 'w (compaction) range : ',minval(ov%wcompact(1:nn)),maxval(ov%wcompact(1:nn))
write(*,'(a,2f30.20)') 'pressure range : ',minval(ov%pressure), maxval(ov%pressure)
write(*,'(a,2f30.20)') 'smooth pressure range : ',minval(ov%spressure), maxval(ov%spressure)
write(*,'(a,2f30.20)') 'pressure (nodes) range : ',minval(ov%nodal_pressure),maxval(ov%nodal_pressure)
......@@ -585,7 +601,7 @@ write(*,'(a,2f30.20)') 'e2d range : ',minval(ov%e2d), ma
write(*,'(a,2f30.20)') 'e2d (nodes) range : ',minval(ov_nodee2d), maxval(ov_nodee2d)
write(*,'(a,2f30.20)') 'yield ratio range : ',minval(ov%yield_ratio), maxval(ov%yield_ratio)
write(*,'(a,2f30.20)') 'friction angle range : ',minval(ov%frict_angle), maxval(ov%frict_angle)
write(*,'(a,2f30.20)') 'compact density range : ',minval(ov%compaction_density), maxval(ov%compaction_density)
if (compactionout) write(*,'(a,2f30.20)') 'compact density range : ',minval(ov%compaction_density), maxval(ov%compaction_density)
write(*,'(a,2f30.20)') 'dilatation rate range : ',minval(ov%dilatr), maxval(ov%dilatr)
write(*,'(a,2f30.20)') 'crit range : ',minval(ov%crit), maxval(ov%crit)
write(*,'(a,2f30.20)') 'crit (nodes) range : ',minval(ov_nodecrit), maxval(ov_nodecrit)
......@@ -952,7 +968,6 @@ end if
!enddo
!close(iunit)
iunit=30
open(unit=iunit,file='visu.vtk')
write(iunit,'(a)')'# vtk DataFile Version 3.0'
......@@ -1129,38 +1144,54 @@ if(output_velo_vect==1) then
if (abs(ov%w(rtf(i))).le.1.d-99) ov%w(rtf(i))=0.d0
write(iunit,'(3e17.6)')ov%u(rtf(i)),ov%v(rtf(i)),ov%w(rtf(i))
enddo
endif
endif
if(output_mech_velo_vect==1) then
write(iunit,'(a)')'VECTORS mech-velo double'
do i=1,nnn
wiso=0.d0
wcompact=0.d0
if (abs(ov%u(rtf(i))).ge.1.d98) ov%u(rtf(i))=sign(1.d98,ov%u(rtf(i)))
if (abs(ov%u(rtf(i))).le.1.d-99) ov%u(rtf(i))=0.d0
if (abs(ov%v(rtf(i))).ge.1.d98) ov%v(rtf(i))=sign(1.d98,ov%v(rtf(i)))
if (abs(ov%v(rtf(i))).le.1.d-99) ov%v(rtf(i))=0.d0
wmod=ov%w(rtf(i))-ov%wiso(rtf(i))-ov%wcompact(rtf(i))
if (isostasyout) wiso=ov%wiso(rtf(i))
if (compactionout) wcompact=ov%wcompact(rtf(i))
wmod=ov%w(rtf(i))-wiso-wcompact
if (abs(wmod).ge.1.d98) wmod=sign(1.d98,wmod)
if (abs(wmod).le.1.d-99) wmod=0.d0
write(iunit,'(3e17.6)')ov%u(rtf(i)),ov%v(rtf(i)),wmod
enddo
endif
if(output_iso_velo_vect==1) then
write(iunit,'(a)')'VECTORS iso-velo double'
do i=1,nnn
if (abs(ov%wiso(rtf(i))).ge.1.d98) ov%wiso(rtf(i))=sign(1.d98,ov%wiso(rtf(i)))
if (abs(ov%wiso(rtf(i))).le.1.d-99) ov%wiso(rtf(i))=0.d0
write(iunit,'(3e17.6)')0.d0,0.d0,ov%wiso(rtf(i))
enddo
if(output_iso_velo_vect==1) then
if (isostasyout) then
write(iunit,'(a)')'VECTORS iso-velo double'
do i=1,nnn
if (abs(ov%wiso(rtf(i))).ge.1.d98) ov%wiso(rtf(i))=sign(1.d98,ov%wiso(rtf(i)))
if (abs(ov%wiso(rtf(i))).le.1.d-99) ov%wiso(rtf(i))=0.d0
write(iunit,'(3e17.6)')0.d0,0.d0,ov%wiso(rtf(i))
enddo
else
write(*,'(a)') '--------------------------------------------------------------------------'
write(*,'(a)') 'WARNING: Isostasy vector output requested, but not found in DOUAR output file.'
write(*,'(a)') ' No isostasy velocity vectors written to visu.vtk file.'
endif
endif
if(output_compaction_velo_vect==1) then
write(iunit,'(a)')'VECTORS compaction-velo double'
do i=1,nnn
if (abs(ov%wcompact(rtf(i))).ge.1.d98) ov%wcompact(rtf(i))=sign(1.d98,ov%wcompact(rtf(i)))
if (abs(ov%wcompact(rtf(i))).le.1.d-99) ov%wcompact(rtf(i))=0.d0
write(iunit,'(3e17.6)')0.d0,0.d0,ov%wcompact(rtf(i))
enddo
if(output_compaction_velo_vect==1) then
if (compactionout) then
write(iunit,'(a)')'VECTORS compaction-velo double'
do i=1,nnn
if (abs(ov%wcompact(rtf(i))).ge.1.d98) ov%wcompact(rtf(i))=sign(1.d98,ov%wcompact(rtf(i)))
if (abs(ov%wcompact(rtf(i))).le.1.d-99) ov%wcompact(rtf(i))=0.d0
write(iunit,'(3e17.6)')0.d0,0.d0,ov%wcompact(rtf(i))
enddo
else
write(*,'(a)') '--------------------------------------------------------------------------'
write(*,'(a)') 'WARNING: Compaction vector output requested, but not found in DOUAR output file.'
write(*,'(a)') ' No compaction velocity vectors written to visu.vtk file.'
endif
endif
if (output_strain_tensor/=0) then
......@@ -1288,17 +1319,23 @@ if (output_frict_angle_fielde==1) then
end if
if (output_comp_density_fielde==1) then
write(iunit,'(a)')'SCALARS comp_density double 1'
write(iunit,'(a)')'LOOKUP_TABLE default'
do ie=1,ov%nleaves
if (elvoid(ie).eq.0) then
if (abs(ov%compaction_density(ie)).ge.1.d98) then
ov%compaction_density(ie)=sign(1.d98,ov%compaction_density(ie))
endif
if (abs(ov%compaction_density(ie)).le.1.d-99) ov%compaction_density(ie)=0.d0
write(iunit,'(e13.6)') ov%compaction_density(ie)
endif
enddo
if (compactionout) then
write(iunit,'(a)')'SCALARS comp_density double 1'
write(iunit,'(a)')'LOOKUP_TABLE default'
do ie=1,ov%nleaves
if (elvoid(ie).eq.0) then
if (abs(ov%compaction_density(ie)).ge.1.d98) then
ov%compaction_density(ie)=sign(1.d98,ov%compaction_density(ie))
endif
if (abs(ov%compaction_density(ie)).le.1.d-99) ov%compaction_density(ie)=0.d0
write(iunit,'(e13.6)') ov%compaction_density(ie)
endif
enddo
else
write(*,'(a)') '--------------------------------------------------------------------------'
write(*,'(a)') 'WARNING: Compaction density output requested, but not found in DOUAR output file.'
write(*,'(a)') ' Compaction density will not be written to visu.vtk file.'
endif
end if
if (output_dilatr_fielde==1) then
......@@ -1609,58 +1646,64 @@ endif
!==============================================================================
if (output_zisodisp==1) then
open(unit=30,file='zisodisp.vtk')
write(30,'(a)')'# vtk DataFile Version 3.0'
write(30,'(a)')'zisodisp'
write(30,'(a)')'ASCII'
write(30,'(a)')'DATASET UNSTRUCTURED_GRID'
write(30,'(a7,i10,a6)')'POINTS ',int((nb+1)**2.d0),' float'
do j=1,nb+1
do i=1,nb+1
write(30,'(3f16.11)') (i-1)*dxy,(j-1)*dxy,zisodisp(i,j)
if (isobcout) then
open(unit=30,file='zisodisp.vtk')
write(30,'(a)')'# vtk DataFile Version 3.0'
write(30,'(a)')'zisodisp'
write(30,'(a)')'ASCII'
write(30,'(a)')'DATASET UNSTRUCTURED_GRID'
write(30,'(a7,i10,a6)')'POINTS ',int((nb+1)**2.d0),' float'
do j=1,nb+1
do i=1,nb+1
write(30,'(3f16.11)') (i-1)*dxy,(j-1)*dxy,zisodisp(i,j)
enddo
enddo
enddo
write(30,'(a6, 2I10)') 'CELLS ',int(nb**2.d0),int((4+1)*nb**2.d0)
ptcnt=0
do j=1,nb
do i=1,nb
write(30,'(9I10)')4,ptcnt,ptcnt+1,ptcnt+1+(nb+1),ptcnt+(nb+1)
write(30,'(a6, 2I10)') 'CELLS ',int(nb**2.d0),int((4+1)*nb**2.d0)
ptcnt=0
do j=1,nb
do i=1,nb
write(30,'(9I10)')4,ptcnt,ptcnt+1,ptcnt+1+(nb+1),ptcnt+(nb+1)
ptcnt=ptcnt+1
enddo
ptcnt=ptcnt+1
enddo
ptcnt=ptcnt+1
enddo
write(30,'(A11, I10)') 'CELL_TYPES ',int(nb**2.d0)
do j=1,nb
do i=1,nb
write(30,'(I2)')9 ! VTK quad (4 nodes)
write(30,'(A11, I10)') 'CELL_TYPES ',int(nb**2.d0)
do j=1,nb
do i=1,nb
write(30,'(I2)')9 ! VTK quad (4 nodes)
enddo
enddo
enddo
write(30,'(a11,i10)') 'POINT_DATA ',int((nb+1)**2.d0)
write(30,'(a)')'SCALARS zisodisp float 1'
write(30,'(a)')'LOOKUP_TABLE default'
do j=1,nb+1
do i=1,nb+1
write(30,'(e13.6)') zisodisp(i,j)
write(30,'(a11,i10)') 'POINT_DATA ',int((nb+1)**2.d0)
write(30,'(a)')'SCALARS zisodisp float 1'
write(30,'(a)')'LOOKUP_TABLE default'
do j=1,nb+1
do i=1,nb+1
write(30,'(e13.6)') zisodisp(i,j)
enddo
enddo
enddo
!write(30,'(a11,i10)') 'POINT_DATA ',int((nb+1)**2.d0)
write(30,'(a)')'SCALARS zisoslx float 1'
write(30,'(a)')'LOOKUP_TABLE default'
do j=1,nb+1
do i=1,nb+1
if (i==1) then
write(30,'(e13.6)') (zisodisp(i+1,j)-zisodisp(i,j))/dxy
elseif (i==nb+1) then
write(30,'(e13.6)') (zisodisp(i,j)-zisodisp(i-1,j))/dxy
else
write(30,'(e13.6)') ((zisodisp(i,j)-zisodisp(i-1,j))/dxy+&
(zisodisp(i+1,j)-zisodisp(i,j))/dxy)/2.d0
endif
!write(30,'(a11,i10)') 'POINT_DATA ',int((nb+1)**2.d0)
write(30,'(a)')'SCALARS zisoslx float 1'
write(30,'(a)')'LOOKUP_TABLE default'
do j=1,nb+1
do i=1,nb+1
if (i==1) then
write(30,'(e13.6)') (zisodisp(i+1,j)-zisodisp(i,j))/dxy
elseif (i==nb+1) then
write(30,'(e13.6)') (zisodisp(i,j)-zisodisp(i-1,j))/dxy
else
write(30,'(e13.6)') ((zisodisp(i,j)-zisodisp(i-1,j))/dxy+&
(zisodisp(i+1,j)-zisodisp(i,j))/dxy)/2.d0
endif
enddo
enddo
enddo
close(30)
write(*,*) '--------------------------------------------------------------------------'
write(*,*)'opla I am done producing zisodisp.vtk'
close(30)
write(*,*) '--------------------------------------------------------------------------'
write(*,*)'opla I am done producing zisodisp.vtk'
else
write(*,'(a)') '--------------------------------------------------------------------------'
write(*,'(a)') 'WARNING: Isostasy displacement output requested, but not found in DOUAR output file.'
write(*,'(a)') ' No isostasy displacement output written to zisodisp.vtk file.'
endif
end if
!==============================================================================
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment