[cig-commits] r18885 - seismo/3D/SPECFEM3D_GEOTECH/trunk/src
homnath at geodynamics.org
homnath at geodynamics.org
Fri Sep 9 00:13:40 PDT 2011
Author: homnath
Date: 2011-09-09 00:13:40 -0700 (Fri, 09 Sep 2011)
New Revision: 18885
Removed:
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semgeotech.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semslope3d.f90
Log:
replaced with .F90 files
Deleted: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semgeotech.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semgeotech.f90 2011-09-09 07:06:01 UTC (rev 18884)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semgeotech.f90 2011-09-09 07:13:40 UTC (rev 18885)
@@ -1,371 +0,0 @@
-! this is a main program SPECFEM3D_GEOTECH
-! REVISION:
-! HNG, Jul 14,2011; HNG, Jul 11,2011; Apr 09,2010
-program semgeotech
-! import necessary libraries
-use global
-use string_library, only : parse_file
-!use math_constants
-use mesh_spec
-#if (USE_MPI)
-use mpi_library
-use math_library_mpi
-#else
-use serial_library
-#endif
-use visual
-
-implicit none
-integer :: funit,i,ios,istat,j,k,neq
-integer :: i_elmt,i_node,i_inc,i_srf,ielmt,igdof,imat,inode
-
-integer :: gnod(8),map2exodus(8),ngllxy,node_hex8(8)
-
-character(len=250) :: arg1,inp_fname,out_fname,prog
-character(len=150) :: path
-character(len=20), parameter :: wild_char='********************'
-character(len=20) :: ensight_etype
-character(len=80) :: buffer,destag ! this must be 80 characters long
-character(len=20) :: ext,format_str,ptail
-character(len=250) :: case_file,geo_file,sum_file
-integer :: npart,nt,tinc,tstart,twidth,ts ! ts: time set for ensight gold
-
-real(kind=kreal) :: cpu_tstart,cpu_tend,telap,step_telap,max_telap,mean_telap
-
-logical :: ismpi !.true. : MPI, .false. : serial
-integer :: ipart,myid,nproc
-integer :: tot_nelmt,max_nelmt,min_nelmt,tot_nnode,max_nnode,min_nnode
-
-character(len=250) :: errtag ! error message
-integer :: errcode
-logical :: isopen ! flag to check whether the file is opened
-
-ismpi=.false.; ipart=0; myid=1; nproc=1; ptail=''
-errtag=""; errcode=-1
-
-call start_process(myid,nproc,stdout)
-ipart=myid-1 ! partition id starts from 0
-
-#if (USE_MPI)
-ismpi=.true.
-#endif
-
-call get_command_argument(0, prog)
-!----input and initialisation----
-if (command_argument_count() <= 0) then
- call error_stop('ERROR: no input file!',stdout,myid)
-endif
-
-call get_command_argument(1, arg1)
-if(trim(arg1)==('--help'))then
- if(myid==1)then
- write(stdout,'(a)')'Usage: '//trim(prog)//' [Options] [input_file]'
- write(stdout,'(a)')'Options:'
- write(stdout,'(a)')' --help : Display this information.'
- write(stdout,'(a)')' --version : Display version information.'
- endif
- !call sync_process
- call close_process()
-elseif(trim(arg1)==('--version'))then
- if(myid==1)then
- write(stdout,'(a)')'SPECFEM3D_GEOTECH 1.0 Beta'
- write(stdout,'(a)')'This is free software; see the source for copying '
- write(stdout,'(a)')'conditions. There is NO warranty; not even for '
- write(stdout,'(a)')'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.'
- endif
- !call sync_process
- call close_process()
-endif
-
-! starting timer
-call cpu_time(cpu_tstart)
-
-! get input file name
-call get_command_argument(1, inp_fname)
-
-! read input data
-call read_input(ismpi,myid,nproc,inp_fname,errcode,errtag)
-if(errcode/=0)call error_stop(errtag,stdout,myid)
-!call sync_process()
-
-tot_nelmt=nelmt; max_nelmt=nelmt; min_nelmt=nelmt
-tot_nnode=nnode; max_nnode=nnode; min_nnode=nnode
-#if (USE_MPI)
-tot_nelmt=sumscal_par(nelmt); tot_nnode=sumscal_par(nnode)
-max_nelmt=maxscal_par(nelmt); max_nnode=maxscal_par(nnode)
-min_nelmt=minscal_par(nelmt); min_nnode=minscal_par(nnode)
-#endif
-if(myid==1)then
-write(stdout,*)'elements => total:',tot_nelmt,' max:',max_nelmt,' min:',min_nelmt
-write(stdout,*)'nodes => total:',tot_nnode,' max:',max_nnode,' min:',min_nnode
-endif
-
-if (trim(method)/='sem')then
- write(errtag,'(a)')'ERROR: wrong input for sem3d!'
- call error_stop(errtag,stdout,myid)
-endif
-
-write(format_str,*)ceiling(log10(real(nproc)+1.))
-format_str='(a,i'//trim(adjustl(format_str))//'.'//trim(adjustl(format_str))//')'
-
-call parse_file(inp_fname,path,file_head,ext)
-
-#if (USE_MPI)
-write(ptail,fmt=format_str)'_proc',ipart
-#endif
-
-if(nexcav==0)then
- nt=nsrf
-else
- nt=nexcav+1 ! include 0 excavation stage (i.e., initial)
-endif
-
-twidth=ceiling(log10(real(nt)+1.))
-! write original meshes
-! write Ensight gold .case file
-ensight_etype='hexa8'
-ts=1 ! time set
-tstart=1; tinc=1
-!write(*,'(a)',advance='no')'writing Ensight case file...'
-
-case_file=trim(out_path)//trim(file_head)//'_original'//trim(ptail)//'.case'
-geo_file=trim(file_head)//'_original'//trim(ptail)//'.geo'
-
-open(unit=11,file=trim(case_file),status='replace',action='write',iostat = ios)
-if( ios /= 0 ) then
- write(errtag,'(a)')'ERROR: file "'//trim(case_file)//'" cannot be opened!'
- call error_stop(errtag,stdout,myid)
-endif
-
-write(11,'(a)')'FORMAT'
-write(11,'(a,/)')'type: ensight gold'
-
-write(11,'(a)')'GEOMETRY'
-write(11,'(a,a,/)')'model: ',trim(geo_file)
-
-close(11)
-
-! write original mesh with 8-noded hexahedra
-! open Ensight Gold geo file to store mesh data
-geo_file=trim(out_path)//trim(geo_file)
-npart=1
-destag='unstructured meshes'
-call write_ensight_geo(geo_file,ensight_etype,destag,npart,nelmt,nnode, &
-real(g_coord),g_num)
-
-! create spectral elements
-if(myid==1)write(stdout,'(a)',advance='no')'creating spectral elements...'
-call hex2spec(ndim,ngnod,nelmt,nnode,ngllx,nglly,ngllz,errcode,errtag)
-if(errcode/=0)call error_stop(errtag,stdout,myid)
-if(myid==1)write(stdout,*)'complete!'
-
-tot_nelmt=nelmt; max_nelmt=nelmt; min_nelmt=nelmt
-tot_nnode=nnode; max_nnode=nnode; min_nnode=nnode
-#if (USE_MPI)
-tot_nelmt=sumscal_par(nelmt); tot_nnode=sumscal_par(nnode)
-max_nelmt=maxscal_par(nelmt); max_nnode=maxscal_par(nnode)
-min_nelmt=minscal_par(nelmt); min_nnode=minscal_par(nnode)
-#endif
-if(myid==1)then
-write(stdout,*)'elements => total:',tot_nelmt,' max:',max_nelmt,' min:',min_nelmt
-write(stdout,*)'nodes => total:',tot_nnode,' max:',max_nnode,' min:',min_nnode
-endif
-
-!call sync_process
-
-!stop
-nenod=ngll !(ngllx*nglly*ngllz) ! number of elemental nodes (nodes per element)
-! number of elemental degrees of freedom
-nedof=nndof*nenod
-
-
-ngllxy=ngllx*nglly
-
-! geometrical nodes (corner nodes) in EXODUS/CUBIT order
-! bottom nodes
-gnod(1)=1;
-gnod(2)=ngllx
-gnod(3)=ngllxy;
-gnod(4)=gnod(3)-ngllx+1
-! top nodes
-gnod(5)=(ngllz-1)*ngllxy+1;
-gnod(6)=gnod(5)+ngllx-1
-gnod(7)=ngll;
-gnod(8)=gnod(7)-ngllx+1
-
-! map sequential node numbering to exodus/cubit order for 8-noded hexahedra
-map2exodus=(/ 1,2,4,3,5,6,8,7 /)
-
-! write Ensight gold .case file
-ensight_etype='hexa8'
-ts=1 ! time set
-tstart=1; tinc=1
-
-case_file=trim(out_path)//trim(file_head)//trim(ptail)//'.case'
-open(unit=11,file=trim(case_file),status='replace',action='write',iostat = ios)
-if( ios /= 0 ) then
- write(errtag,'(a)')'ERROR: file "'//trim(case_file)//'" cannot be opened!'
- call error_stop(errtag,stdout,myid)
-endif
-
-write(11,'(a)')'FORMAT'
-write(11,'(a,/)')'type: ensight gold'
-
-write(11,'(a)')'GEOMETRY'
-!write(11,'(a,a,/)')'model: ',trim(geo_file)
-if(nexcav==0)then
- write(11,'(a,a/)')'model: ',trim(file_head)//trim(ptail)//'.geo'
-else
- write(11,'(a,i10,a,a/)')'model: ',ts,' ',trim(file_head)//'_step'//wild_char(1:twidth)//trim(ptail)//'.geo'
-endif
-
-write(11,'(a)')'VARIABLE'
-!write(11,'(a,i10,a,a,a,a,/)')'vector per node: ',ts,' displacement ', &
-!trim(file_head)//'_step'//wild_char(1:twidth)//'.dis'
-!write(11,'(a,i10,a,a,a,a,/)')'vector per node: ',ts,' principal_stress ', &
-!trim(file_head)//'_step'//wild_char(1:twidth)//'.sig'
-
-if(savedata%disp)then
- write(11,'(a,i10,a,a,a,a,/)')'vector per node: ',ts,' ','displacement',' ', &
- trim(file_head)//'_step'//wild_char(1:twidth)//trim(ptail)//'.dis'
-endif
-if(savedata%stress)then
- write(11,'(a,i10,a,a,a,a,/)')'tensor symm per node: ',ts,' ','stress',' ', &
- trim(file_head)//'_step'//wild_char(1:twidth)//trim(ptail)//'.sig'
-endif
-if(savedata%psigma)then
- write(11,'(a,i10,a,a,a,a,/)')'vector per node: ',ts,' ','principal_stress',' ', &
- trim(file_head)//'_step'//wild_char(1:twidth)//trim(ptail)//'.psig'
-endif
-if(savedata%porep)then
- write(11,'(a,a,a,a,a,/)')'scalar per node: ',' ','pore_pressure',' ', &
- trim(file_head)//trim(ptail)//'.por'
-endif
-if(savedata%vmeps)then
- write(11,'(a,i10,a,a,a,a,/)')'scalar per node: ',ts,' ','plastic_strain',' ', &
- trim(file_head)//'_step'//wild_char(1:twidth)//trim(ptail)//'.eps'
-endif
-if(savedata%scf)then
- write(11,'(a,i10,a,a,a,a,/)')'scalar per node: ',ts,' ','stress_concentration_factor',' ', &
- trim(file_head)//'_step'//wild_char(1:twidth)//trim(ptail)//'.scf'
-endif
-if(savedata%maxtau)then
- write(11,'(a,i10,a,a,a,a,/)')'scalar per node: ',ts,' ','max_shear_stress',' ', &
- trim(file_head)//'_step'//wild_char(1:twidth)//trim(ptail)//'.mtau'
-endif
-if(savedata%nsigma)then
- write(11,'(a,i10,a,a,a,a,/)')'scalar per node: ',ts,' ','normal_stress',' ', &
- trim(file_head)//'_step'//wild_char(1:twidth)//trim(ptail)//'.nsig'
-endif
-write(11,'(a)')'TIME'
-write(11,'(a,i10)')'time set:',ts
-write(11,'(a,i10)')'number of steps:',nsrf
-write(11,'(a,i10)')'filename start number:',tstart
-write(11,'(a,i10)')'filename increment:',tinc
-write(11,'(a)',advance='no')'time values: '
-
-if(nexcav==0)then
- do i=1,nt
- write(11,'(e12.5)',advance='yes')srf(i)
- enddo
-else
- do i=0,nt-1
- write(11,'(e12.5)',advance='yes')real(i)
- enddo
-endif
-close(11)
-
-! Format string
-write(format_str,*)twidth
-format_str='(a,i'//trim(adjustl(format_str))//'.'//trim(adjustl(format_str))//',a)'
-
-! write geo file for inital stage (original)
-! open Ensight Gold geo file to store mesh data
-if(nexcav==0)then
- write(geo_file,fmt=format_str)trim(out_path)//trim(file_head)//trim(ptail)//'.geo'
-else
- write(geo_file,fmt=format_str)trim(out_path)//trim(file_head)//'_step',0,trim(ptail)//'.geo'
-endif
-npart=1
-destag='unstructured meshes'
-call write_ensight_geocoord(geo_file,destag,npart,nnode,real(g_coord),funit)
-
-! writes element information
-buffer=ensight_etype
-write(funit)buffer
-write(funit)nelmt*(ngllx-1)*(nglly-1)*(ngllz-1)
-
-! do not substract 1 for ensight file
-do i_elmt=1,nelmt
- do k=1,ngllz-1
- do j=1,nglly-1
- do i=1,ngllx-1
- ! corner nodes in a sequential numbering
- node_hex8(1)=(k-1)*ngllxy+(j-1)*ngllx+i
- node_hex8(2)=node_hex8(1)+1
-
- node_hex8(3)=node_hex8(1)+ngllx
- node_hex8(4)=node_hex8(3)+1
-
- node_hex8(5)=node_hex8(1)+ngllxy
- node_hex8(6)=node_hex8(5)+1
-
- node_hex8(7)=node_hex8(5)+ngllx
- node_hex8(8)=node_hex8(7)+1
- ! map to exodus/cubit numbering and write
- write(funit)g_num(node_hex8(map2exodus),i_elmt)
- enddo
- enddo
- enddo
-enddo
-close(funit)
-
-! open summary file
-sum_file = trim(out_path)//trim(file_head)//'_summary'//trim(ptail)
-open(unit=10,file=trim(sum_file),status='replace',action='write',iostat=ios)
-write(10,*)'--------------------------------------------'
-write(10,*)'Result summary produced by SPECFEM3D_GEOTECH'
-write(10,*)'--------------------------------------------'
-close(10)
-
-if(myid==1)write(stdout,'(a)')'--------------------------------------------'
-
-! call main routines
-if(nexcav==0)then
- ! slope stability
- call semslope3d(ismpi,myid,nproc,gnod,sum_file,ptail,format_str)
-else
- ! excavation
- !call semexcav3d(ismpi,myid,nproc,gnod,sum_file,ptail,format_str)
-endif
-!-----------------------------------
-
-! compute elapsed time
-call cpu_time(cpu_tend)
-telap=cpu_tend-cpu_tstart
-max_telap=telap; mean_telap=telap
-#if (USE_MPI)
-max_telap=maxscal_par(telap)
-mean_telap=sumscal_par(telap)/real(nproc,kreal)
-#endif
-
-write(format_str,*)ceiling(log10(real(max_telap)+1.))+5 ! 1 . and 4 decimals
-format_str='(3(f'//trim(adjustl(format_str))//'.4,1X))'
-open(10,file=trim(sum_file),status='old',position='append',action='write')
-write(10,*)'ELAPSED TIME, MAX ELAPSED TIME, MEAN ELAPSED TIME'
-write(10,fmt=format_str)telap,max_telap,mean_telap
-close(10)
-!-----------------------------------
-
-if(myid==1)then
- write(stdout,*) ! write new line
- write(stdout,'(a)')'--------------------------------------------'
- inquire(stdout,opened=isopen)
- if(isopen)close(stdout)
-endif
-
-call sync_process
-call close_process()
-
-end program semgeotech
-!===========================================
Deleted: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semslope3d.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semslope3d.f90 2011-09-09 07:06:01 UTC (rev 18884)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semslope3d.f90 2011-09-09 07:13:40 UTC (rev 18885)
@@ -1,513 +0,0 @@
-! this is a main routine for slope stabiliy analysis
-! REVISION:
-! HNG, Jul 14,2011; HNG, Jul 11,2011; Apr 09,2010
-subroutine semslope3d(ismpi,myid,nproc,gnod,sum_file,ptail,format_str)
-! import necessary libraries
-use global
-use string_library, only : parse_file
-use math_constants
-use gll_library
-!use mesh_spec
-use shape_library
-use math_library
-use preprocess
-!use gauss_library
-!use excavation
-use plastic_library
-#if (USE_MPI)
-use mpi_library
-use ghost_library_mpi
-use math_library_mpi
-use solver_mpi
-#else
-use serial_library
-use solver
-#endif
-use visual
-use postprocess
-
-implicit none
-logical,intent(in) :: ismpi
-integer,intent(in) :: myid,nproc
-integer,intent(in) :: gnod(8)
-character(len=250),intent(in) :: sum_file
-character(len=20),intent(in) :: ptail,format_str
-
-integer :: funit,i,ios,istat,j,k,neq
-integer :: i_elmt,i_node,i_inc,i_srf,ielmt,igdof,imat,inode
-!real(kind=kreal),parameter :: two_third=two/r3
-real(kind=kreal) :: detjac,dq1,dq2,dq3,dsbar,dt,f,fmax,h1,h2,lode_theta,sf,sigm
-
-real(kind=kreal) :: uerr,umax,uxmax
-integer :: cg_iter,cg_tot,nl_iter,nl_tot
-logical :: nl_isconv ! logical variable to check convergence of
-! nonlinear (NL) iterations
-
-real(kind=kreal) :: cmat(nst,nst),devp(nst),eps(nst),erate(nst),evp(nst), &
-flow(nst,nst),m1(nst,nst),m2(nst,nst),m3(nst,nst),effsigma(nst),sigma(nst)
-! dynamic arrays
-integer,allocatable::gdof(:,:),gdof_elmt(:,:),num(:),node_valency(:)
-! factored parameters
-real(kind=kreal),allocatable :: cohf(:),nuf(:),phif(:),psif(:),ymf(:)
-real(kind=kreal),allocatable::bodyload(:),bmat(:,:),bload(:),coord(:,:), &
-der(:,:),deriv(:,:),dprecon(:),eld(:),eload(:),evpt(:,:,:),excavload(:,:), &
-extload(:),jac(:,:),load(:),nodalu(:,:),storkm(:,:,:),oldx(:),x(:), &
-stress_local(:,:,:),stress_global(:,:),scf(:),vmeps(:)
-!,psigma(:,:),psigma0(:,:),taumax(:),nsigma(:)
-integer,allocatable :: egdof(:) ! elemental global degree of freedom
-
-integer :: map2exodus(8),ngllxy,node_hex8(8)
-real(kind=kreal),allocatable :: dshape_hex8(:,:,:)
-real(kind=kreal),parameter :: jacobi_alpha=0.0_kreal,jacobi_beta=0.0_kreal
-!double precision
-real(kind=kreal),allocatable :: xigll(:),wxgll(:) !double precision
-real(kind=kreal),allocatable :: etagll(:),wygll(:) !double precision
-real(kind=kreal),allocatable :: zetagll(:),wzgll(:) !double precision
-real(kind=kreal),allocatable :: gll_weights(:),gll_points(:,:)
-real(kind=kreal),allocatable :: lagrange_gll(:,:),dlagrange_gll(:,:,:)
-
-character(len=250) :: inp_fname,out_fname,prog
-character(len=150) :: path
-character(len=20), parameter :: wild_char='********************'
-character(len=20) :: ensight_etype
-character(len=80) :: buffer,destag ! this must be 80 characters long
-character(len=20) :: ext !,format_str !,ptail
-character(len=250) :: case_file,geo_file !,sum_file
-integer :: npart,tinc,tstart,twidth,ts ! ts: time set for ensight gold
-
-real(kind=kreal) :: cpu_tstart,cpu_tend,telap,step_telap,max_telap,mean_telap
-
-logical :: gravity,pseudoeq ! gravity load and pseudostatic load
-real(kind=kreal),allocatable :: wpressure(:) ! water pressure
-logical,allocatable :: submerged_node(:)
-
-!logical :: ismpi !.true. : MPI, .false. : serial
-integer :: ipart !,myid,nproc
-integer :: tot_nelmt,max_nelmt,min_nelmt,tot_nnode,max_nnode,min_nnode
-integer :: tot_neq,max_neq,min_neq
-integer :: ngpart,maxngnode
-! number of active ghost partitions for a node
-integer,allocatable :: ngpart_node(:)
-character(len=250) :: errtag ! error message
-integer :: errcode
-logical :: isopen ! flag to check whether the file is opened
-
-errtag=""; errcode=-1
-
-ipart=myid-1 ! partition id starts from 0
-
-! apply displacement boundary conditions
-if(myid==1)write(stdout,'(a)',advance='no')'applying BC...'
-allocate(gdof(nndof,nnode),stat=istat)
-if (istat/=0)then
- write(stdout,*)'ERROR: cannot allocate memory!'
- stop
-endif
-gdof=1
-call apply_bc(ismpi,myid,nproc,gdof,neq,errcode,errtag)
-if(errcode/=0)call error_stop(errtag,stdout,myid)
-if(myid==1)write(stdout,*)'complete!'
-!-------------------------------------
-
-allocate(num(nenod),evpt(nst,ngll,nelmt),coord(ngnod,ndim), &
-jac(ndim,ndim),der(ndim,ngnod),deriv(ndim,nenod),bmat(nst,nedof), &
-eld(nedof),bload(nedof),eload(nedof),nodalu(nndof,nnode), &
-egdof(nedof),stat=istat)
-if (istat/=0)then
- write(stdout,*)'ERROR: cannot allocate memory!'
- stop
-endif
-
-!neq=maxval(gdof)
-tot_neq=neq; max_neq=neq; min_neq=neq
-#if (USE_MPI)
-tot_neq=sumscal_par(neq); max_neq=maxscal_par(neq); min_neq=minscal_par(neq)
-#endif
-if(myid==1)then
- write(stdout,*)'degrees of freedoms => total:',tot_neq,' max:',max_neq, &
- ' min:',min_neq
-endif
-
-! get GLL quadrature points and weights
-allocate(xigll(ngllx),wxgll(ngllx),etagll(nglly),wygll(nglly),zetagll(ngllz), &
-wzgll(ngllz))
-call zwgljd(xigll,wxgll,ngllx,jacobi_alpha,jacobi_beta)
-call zwgljd(etagll,wygll,nglly,jacobi_alpha,jacobi_beta)
-call zwgljd(zetagll,wzgll,ngllz,jacobi_alpha,jacobi_beta)
-
-! get derivatives of shape functions for 8-noded hex
-allocate(dshape_hex8(ndim,ngnod,ngll))
-call dshape_function_hex8(ndim,ngnod,ngllx,nglly,ngllz,xigll,etagll,zetagll, &
-dshape_hex8)
-deallocate(xigll,wxgll,etagll,wygll,zetagll,wzgll)
-! compute gauss-lobatto-legendre quadrature information
-allocate(gll_weights(ngll),gll_points(ndim,ngll),lagrange_gll(ngll,ngll), &
-dlagrange_gll(ndim,ngll,ngll))
-call gll_quadrature(ndim,ngllx,nglly,ngllz,ngll,gll_points,gll_weights, &
-lagrange_gll,dlagrange_gll)
-!--------------------------------
-
-! store elemental global degrees of freedoms from nodal gdof
-! this removes the repeated use of reshape later but it has larger size than gdof!!!
-allocate(gdof_elmt(nedof,nelmt))
-gdof_elmt=0
-do i_elmt=1,nelmt
- gdof_elmt(:,i_elmt)=reshape(gdof(:,g_num(:,i_elmt)),(/nedof/)) !g=g_g(:,i_elmt)
-enddo
-!-------------------------------
-
-! compute stiffness and body load
-if(myid==1)write(stdout,'(a)',advance='no')'preprocessing...'
-
-allocate(stress_local(nst,ngll,nelmt))
-! compute initial stress assuming elastic domain
-stress_local=zero
-
-allocate(extload(0:neq),dprecon(0:neq),storkm(nedof,nedof,nelmt),stat=istat)
-! elastic(0:neq),
-if (istat/=0)then
- write(stdout,*)'ERROR: cannot allocate memory!'
- stop
-endif
-extload=zero; gravity=.true.; pseudoeq=iseqload
-call stiffness_bodyload(nelmt,neq,gnod,g_num,gdof_elmt,mat_id,gam,nu,ym, &
-dshape_hex8,lagrange_gll,dlagrange_gll,gll_weights,storkm,dprecon,extload, &
-gravity,pseudoeq)
-
-!print*,minval(dprecon),maxval(dprecon)
-!print*,minval(extload),maxval(extload)
-!print*,minval(storkm),maxval(storkm)
-!stop
-if(myid==1)write(stdout,*)'complete!'
-!-------------------------------
-
-! apply traction boundary conditions
-if(istraction)then
- if(myid==1)write(*,'(a)',advance='no')'applying traction...'
- call apply_traction(ismpi,myid,nproc,gnod,gdof,neq,extload,errcode,errtag)
- if(errcode/=0)call error_stop(errtag,stdout,myid)
- if(myid==1)write(*,*)'complete!'
-endif
-!-------------------------------
-
-! compute water pressure
-if(iswater)then
- if(myid==1)write(stdout,'(a)',advance='no')'computing water pressure...'
- allocate(wpressure(nnode),submerged_node(nnode))
- call compute_pressure(ismpi,myid,nproc,wpressure,submerged_node,errcode,errtag)
- if(errcode/=0)call error_stop(errtag,stdout,myid)
- ! write pore pressure file
-
- ! Open Ensight Gold data file to store data
- out_fname=trim(out_path)//trim(file_head)//trim(ptail)//'.por'
- npart=1;
- destag='Pore pressure'
- call write_ensight_pernode(out_fname,destag,npart,1,nnode,real(wpressure))
- if(myid==1)write(stdout,*)'complete!'
-endif
-!-------------------------------
-
-if(myid==1)write(stdout,'(a)')'--------------------------------------------'
-
-#if (USE_MPI)
-! prepare ghost partitions for the communication
-call prepare_ghost(myid,nproc,gdof,ngpart,maxngnode)
-#endif
-
-#if (USE_MPI)
-! assemble from ghost partitions
-call assemble_ghosts(myid,ngpart,maxngnode,nndof,neq,dprecon,dprecon)
-!fmax=maxvec_par(abs(dprecon))
-!if(myid==1)print*,fmax
-!stop
-#endif
-!print*,minval(dprecon),maxval(dprecon)
-!print*,minval(storkm),maxval(storkm)
-!stop
-dprecon(1:)=one/dprecon(1:); dprecon(0)=zero
-
-allocate(stress_global(nst,nnode),vmeps(nnode))
-allocate(scf(nnode))
-!,psigma(ndim,nnode),psigma0(ndim,nnode),taumax(nnode),nsigma(nnode))
-allocate(node_valency(nnode))
-
-! compute node valency only once
-node_valency=0
-do i_elmt=1,nelmt
- num=g_num(:,i_elmt)
- node_valency(num)=node_valency(num)+1
-enddo
-
-! open summary file
-open(unit=10,file=trim(sum_file),status='old',position='append',action='write',iostat=ios)
-write(10,*)'CG_MAXITER, CG_TOL, NL_MAXITER, NL_TOL'
-write(10,*)cg_maxiter,cg_tol,nl_maxiter,nl_tol
-write(10,*)'Number of SRFs'
-write(10,*)nsrf
-write(10,*)'SRF, CGITER, NLITER, UXMAX, UMAX, fmax'
-close(10)
-
-allocate(bodyload(0:neq),load(0:neq),x(0:neq),oldx(0:neq),stat=istat)
-! elastic(0:neq),
-if (istat/=0)then
- write(stdout,*)'ERROR: cannot allocate memory!'
- stop
-endif
-
-allocate(cohf(nmat),nuf(nmat),phif(nmat),psif(nmat),ymf(nmat))
-
-if(myid==1)then
- write(stdout,'(a,e10.4,a,i5)')'CG_TOL:',cg_tol,' CG_MAXITER:',cg_maxiter
- write(stdout,'(a,e10.4,a,i5)')'NL_TOL:',nl_tol,' NL_MAXITER:',nl_maxiter
- write(stdout,'(a)',advance='no')'SRFs:'
- do i_srf=1,nsrf
- write(stdout,'(f7.4,1x)',advance='no')srf(i_srf)
- enddo
- write(stdout,'(/,a)')'--------------------------------------------'
-endif
-
-! strength reduction (factor of safety) loop
-srf_loop: do i_srf=1,nsrf
- if(myid==1)write(stdout,'(/,a,f7.4)')'SRF:',srf(i_srf)
-
- ! initialize
- nodalu=zero; vmeps=zero
- stress_local=zero; scf=inftol
-
- ! strength reduction
- call strength_reduction(srf(i_srf),phinu,nmat,coh,nu,phi,psi,cohf,nuf,phif, &
- psif,istat)
-
- ! compute minimum pseudo-time step for viscoplasticity
- dt=dt_viscoplas(nmat,nuf,phif,ym)
-
- ! recompute stiffness if either of nu and ym has changed
- if(istat==1)then
- ! in future this should be changed so that only the elements with changed
- ! material properties are involved
- dprecon=zero
- call stiffness_bodyload(nelmt,neq,gnod,g_num,gdof_elmt,mat_id,gam,nuf,ym, &
- dshape_hex8, &
-lagrange_gll,dlagrange_gll,gll_weights,storkm,dprecon)
-
-#if (USE_MPI)
- ! assemble from ghost partitions
- call assemble_ghosts(myid,ngpart,maxngnode,nndof,neq,dprecon,dprecon)
-#endif
- dprecon(1:)=one/dprecon(1:); dprecon(0)=zero
- endif
-
- !print*,nsrf,srf(i_srf),nuf,phif,dt
- !print*,sin(phif*deg2rad),one-two*nuf
-
-#if (USE_MPI)
- ! find global dt
- dt=minscal_par(dt)
-#endif
- cg_tot=0; nl_tot=0
- ! load incremental loop
- !if(myid==1)write(stdout,'(a,i10)')' total load increments:',ninc
- !extload=extload/ninc
- !load_increment: do i_inc=1,ninc
-
- bodyload=zero; evpt=zero
- x=zero; oldx=zero
-
- !print*,maxval(abs(bodyload)),maxval(abs(extload)),maxval(abs(dprecon))
-
- ! plastic iteration loop
- plastic: do nl_iter=1,nl_maxiter
- fmax=zero
-
- load=extload+bodyload
- load(0)=zero
-
- ! pcg solver
- !x=zero
-#if (USE_MPI)
- !call pcg_solver_par(myid,ngpart,maxngnode,neq,storkm,x,load,dprecon_g, &
- !gdof_elmt,cg_iter)
- call pcg_solver_par(myid,ngpart,maxngnode,neq,nelmt,storkm,x,load,dprecon, &
- gdof_elmt,cg_iter,errcode,errtag)
-#else
- call pcg_solver(neq,nelmt,storkm,x,load,dprecon,gdof_elmt,cg_iter,errcode, &
- errtag)
-#endif
- if(errcode/=0)call error_stop(errtag,stdout,myid)
- cg_tot=cg_tot+cg_iter
- x(0)=zero
-
- if(allelastic)then
- call elastic_stress(nelmt,neq,gnod,g_num,gdof_elmt,mat_id,dshape_hex8, &
- dlagrange_gll,x,stress_local)
-
- exit plastic
- endif
-
- ! check plastic convergence
-#if (USE_MPI)
- uerr=maxvec_par(abs(x-oldx))/maxvec_par(abs(x))
-#else
- uerr=maxval(abs(x-oldx))/maxval(abs(x))
-#endif
- oldx=x
- nl_isconv=uerr.le.nl_tol
-
- ! compute stress and check failure
- do i_elmt=1,nelmt
- ielmt=i_elmt
- imat=mat_id(ielmt)
-
- call compute_cmat(cmat,ym(imat),nuf(imat))
- num=g_num(:,ielmt)
- coord=transpose(g_coord(:,num(gnod))) !transpose(g_coord(:,num(1:ngnod)))
- egdof=gdof_elmt(:,ielmt)
- !reshape(gdof(:,g_num(:,ielmt)),(/nedof/)) !g=g_g(:,i_elmt)
- eld=x(egdof)
-
- bload=zero
- do i=1,ngll ! loop over integration points
- jac=matmul(dshape_hex8(:,:,i),coord)
- detjac=determinant(jac)
- call invert(jac)
-
- deriv=matmul(jac,dlagrange_gll(:,i,:))
- call compute_bmat(bmat,deriv)
- eps=matmul(bmat,eld)
- eps=eps-evpt(:,i,ielmt)
- sigma=matmul(cmat,eps)
-
- ! compute effective stress
- effsigma=sigma+stress_local(:,i,ielmt)
- if(iswater)then
- if(submerged_node(num(i)))then
- ! water pressure is compressive (negative)
- effsigma(1:3)=effsigma(1:3)+wpressure(num(i))
- endif
- endif
-
- !effsigma=effsigma+stress_local(:,i,ielmt)
- call stress_invariant(effsigma,sigm,dsbar,lode_theta)
- ! check whether yield is violated
- call mohcouf(phif(imat),cohf(imat),sigm,dsbar,lode_theta,f)
- if(f>fmax)fmax=f
-
- if(f>=zero)then !.or.(nl_isconv.or.nl_iter==nl_maxiter))then
- call mohcouq(psif(imat),dsbar,lode_theta,dq1,dq2,dq3)
- call formm(effsigma,m1,m2,m3)
- flow=f*(m1*dq1+m2*dq2+m3*dq3)
-
- erate=matmul(flow,effsigma)
- evp=erate*dt
- evpt(:,i,ielmt)=evpt(:,i,ielmt)+evp
- devp=matmul(cmat,evp)
- ! if not converged we need body load for next iteration
- if(.not.nl_isconv .and. nl_iter/=nl_maxiter)then
- !devp(1:3)=devp(1:3)-wpressure(num(i))
- eload=matmul(devp,bmat)
- bload=bload+eload*detjac*gll_weights(i)
- end if
- end if
- if(nl_isconv.or.nl_iter==nl_maxiter)then
- devp=sigma
- ! compute von Mises effective plastic strain
- !vmeps(num(i))=vmeps(num(i))+sqrt(two_third* &
- !dot_product(evpt(:,i,ielmt),evpt(:,i,ielmt)))
- ! update stresses
- stress_local(:,i,ielmt)=effsigma
- !phifr=atan(tnph/srf(i_srf))
- !sf=(sigm*sin(phifr)-cohf*cos(phifr))/(-dsbar*(cos(lode_theta)/ &
- !sqrt(r3)-sin(lode_theta)*sin(phifr)/r3))
- !if(sf<scf(num(i)))scf(num(i))=sf
- endif
- end do ! i_gll
-
- if(nl_isconv .or. nl_iter==nl_maxiter)cycle
- ! compute total body load vector
- bodyload(egdof)=bodyload(egdof)+bload
- end do ! i_elmt
- bodyload(0)=zero
-#if (USE_MPI)
- fmax=maxscal_par(fmax)
- uxmax=maxvec_par(abs(x))
-#else
- uxmax=maxval(abs(x))
-#endif
- if(myid==1)then
- write(stdout,'(a,a,i4,a,f12.6,a,f12.6,a,f12.6)',advance='no')CR, &
- ' nl_iter:',nl_iter,' f_max:',fmax,' uerr:',uerr,' umax:',uxmax
- endif
- if(nl_isconv.or.nl_iter==nl_maxiter)exit
- end do plastic ! plastic iteration
- ! check if the plastic iteration did not converge
- if(nl_iter>=nl_maxiter .and. .not.nl_isconv)then
- write(stdout,*)'WARNING: nonconvergence in nonlinear iterations!'
- write(stdout,*)'desired tolerance:',nl_tol,' achieved tolerance:',uerr
- endif
-
- nl_tot=nl_tot+nl_iter
- !if(myid==1)print*,cg_tot,nl_tot
- ! nodal displacement
- do i=1,nndof
- do j=1,nnode
- if(gdof(i,j)/=0)then
- nodalu(i,j)=nodalu(i,j)+x(gdof(i,j)) !-elastic(gdof(i,j))
- endif
- enddo
- enddo
- !enddo load_increment ! load increment loop
-
- ! write summary
-#if (USE_MPI)
- uxmax=maxvec_par(abs(reshape(nodalu,(/nndof*nnode/))))
- umax=maxvec_par(sqrt(nodalu(1,:)*nodalu(1,:)+ &
- nodalu(2,:)*nodalu(2,:)+nodalu(3,:)*nodalu(3,:)))
-#else
- uxmax=maxval(abs(reshape(nodalu,(/nndof*nnode/))))
- umax=maxval(sqrt(nodalu(1,:)*nodalu(1,:)+ &
- nodalu(2,:)*nodalu(2,:)+nodalu(3,:)*nodalu(3,:)))
-#endif
- open(10,file=trim(sum_file),status='old',position='append',action='write')
- write(10,*)srf(i_srf),cg_tot,nl_tot,uxmax,umax,fmax
- close(10)
-
- ! compute average effective strain
- do i_node=1,nnode
- inode=i_node
- vmeps(inode)=vmeps(inode)/real(node_valency(inode),kreal)
- enddo
-
- ! compute stress_global
- stress_global=zero
- do i_elmt=1,nelmt
- ielmt=i_elmt
- num=g_num(:,ielmt)
- stress_global(:,num)=stress_global(:,num)+stress_local(:,:,ielmt)
- enddo
-
- ! compute average stress at sharing nodes
- do i_node=1,nnode
- inode=i_node
- stress_global(:,inode)=stress_global(:,inode)/real(node_valency(inode),kreal)
- enddo
-
- call save_data(ptail,format_str,i_srf,nnode,nelmt,g_num, &
- nodalu,scf,vmeps,stress_global)
-
- if(nl_iter==nl_maxiter)exit
-
-enddo srf_loop ! i_srf safety factor loop
-deallocate(mat_id,gam,ym,coh,nu,phi,psi,srf)
-deallocate(g_coord,g_num)
-deallocate(load,bodyload,extload,oldx,x,dprecon,storkm,stat=istat)
-#if (USE_MPI)
- do i=1,ngpart
- deallocate(gpart(i)%node,gpart(i)%gdof)
- enddo
- deallocate(gpart)
-#endif
-!-----------------------------------
-
-return
-end subroutine semslope3d
-!===========================================
More information about the CIG-COMMITS
mailing list