[cig-commits] [commit] master: Making relaxlite library (7338c32)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri Jan 30 19:49:56 PST 2015
Repository : https://github.com/geodynamics/relax
On branch : master
Link : https://github.com/geodynamics/relax/compare/17d2ea614ff358822f81bb94d3774925a60710fb...7338c32a4753447c8df77940d30cacde6530a669
>---------------------------------------------------------------
commit 7338c32a4753447c8df77940d30cacde6530a669
Author: sagar masuti <sagar.masuti at gmail.com>
Date: Sat Jan 31 11:49:27 2015 +0800
Making relaxlite library
>---------------------------------------------------------------
7338c32a4753447c8df77940d30cacde6530a669
src/input.f90 | 2 +-
src/relax.f90 | 2 +-
src/{curelax.f90 => relaxlite.f90} | 558 +++++++++++++++----------------------
src/types.f90 | 7 +-
wscript | 42 ++-
5 files changed, 265 insertions(+), 346 deletions(-)
diff --git a/src/input.f90 b/src/input.f90
index df3ca72..15c08c2 100644
--- a/src/input.f90
+++ b/src/input.f90
@@ -43,7 +43,7 @@ CONTAINS
USE export
USE getopt_m
- TYPE(SIMULATION_STRUC), INTENT(OUT) :: in
+ TYPE(SIMULATION_STRUCT), INTENT(OUT) :: in
CHARACTER :: ch
CHARACTER(256) :: dataline
diff --git a/src/relax.f90 b/src/relax.f90
index 05dd762..ea6c89d 100644
--- a/src/relax.f90
+++ b/src/relax.f90
@@ -212,7 +212,7 @@ PROGRAM relax
!$ INTEGER :: omp_get_max_threads
#endif
REAL*8 :: maxwell(3)
- TYPE(SIMULATION_STRUC) :: in
+ TYPE(SIMULATION_STRUCT) :: in
#ifdef VTK
CHARACTER(256) :: filename,title,name
CHARACTER(3) :: digit
diff --git a/src/curelax.f90 b/src/relaxlite.f90
similarity index 73%
copy from src/curelax.f90
copy to src/relaxlite.f90
index e689857..fb68b29 100644
--- a/src/curelax.f90
+++ b/src/relaxlite.f90
@@ -1,5 +1,5 @@
!-----------------------------------------------------------------------
-! Copyright 2007-2012, Sylvain Barbot
+! Copyright 2007-2013, Sylvain Barbot
!
! This file is part of RELAX
!
@@ -176,6 +176,8 @@
!! and allow scaling of computed time steps. <br>
!! (04-26-11) - include command-line arguments
!! (11-04-11) - compatible with gfortran <br>
+ !! (07-25-13) - include cylindrical and spherical ductile zones <br>
+ !! (01-03-14) - accelerate processing of ductile zones <br>
!!
!! \todo
!! - homogenize VTK output so that geometry of events match event index
@@ -187,11 +189,11 @@
!! - check the fully-relaxed afterslip for uniform stress change
!! - include topography of parameter interface
!! - export afterslip output in VTK legacy format (binary)
+ !! - export ductile zones for cylindrical and spherical geometries
!------------------------------------------------------------------------
-PROGRAM relax
-
+SUBROUTINE relaxlite(in,gps,isverbose)
USE types
- USE input
+! USE input
USE green
USE green_space
USE elastic3d
@@ -202,98 +204,70 @@ PROGRAM relax
#include "include.f90"
IMPLICIT NONE
+ TYPE(SIMULATION_STRUCT), INTENT(INOUT) :: in
+ TYPE(MANIFOLD_STRUCT), INTENT(OUT) :: gps(in%npts)
+ LOGICAL, INTENT(IN) :: isverbose
- INTEGER, PARAMETER :: ITERATION_MAX = 99999
+ INTEGER, PARAMETER :: ITERATION_MAX = 9999
REAL*8, PARAMETER :: STEP_MAX = 1e7
- INTEGER :: i,k,e,oi,iostatus,mech(3),iTensor,iRet
-
+ INTEGER :: i,k,e,oi,iostatus,mech(3)
+#ifdef FFTW3_THREADS
+ INTEGER :: iret
+!$ INTEGER :: omp_get_max_threads
+#endif
REAL*8 :: maxwell(3)
- TYPE(SIMULATION_STRUC) :: in
-#ifdef VTK
CHARACTER(256) :: filename,title,name
- CHARACTER(3) :: digit
-#endif
+
CHARACTER(4) :: digit4
- REAL*8 :: t,Dt,tm,dAmp,dAmp1
- REAL*4 :: cuC1, cuC2
+ REAL*8 :: t,Dt,tm
! arrays
REAL*4, DIMENSION(:,:,:), ALLOCATABLE :: v1,v2,v3,u1,u2,u3,gamma
- REAL*4, DIMENSION(:,:,:), ALLOCATABLE :: u1r,u2r,u3r
+ REAL*4, DIMENSION(:,:,:), ALLOCATABLE :: lineardgammadot0,nonlineardgammadot0
REAL*4, DIMENSION(:,:), ALLOCATABLE :: t1,t2,t3
- REAL*4, DIMENSION(:,:,:), ALLOCATABLE :: inter1,inter2,inter3
TYPE(TENSOR), DIMENSION(:,:,:), ALLOCATABLE :: tau,sig,moment
-#ifdef PAPI_PROF
- CHARACTER (LEN=16) cTimerName
- cTimerName = 'relax'
+#ifdef FFTW3_THREADS
+ CALL sfftw_init_threads(iret)
+#ifdef _OPENMP
+ CALL sfftw_plan_with_nthreads(omp_get_max_threads())
+#else
+ CALL sfftw_plan_with_nthreads(4)
+#endif
#endif
- ! read input parameters
- CALL init(in)
-
- ! abort calculation after help message
- ! or for dry runs
- IF (in%isdryrun) THEN
- PRINT '("dry run: abort calculation")'
- END IF
- IF (in%isdryrun .OR. in%ishelp) THEN
- ! exit program
- GOTO 100
- END IF
-
- CALL cuinit (%VAL(in%sx1), %VAL(in%sx2), %VAL(in%sx3), %VAL(in%dx1), %VAL(in%dx2), %VAL(in%dx3), iostatus)
- IF (iostatus>0) STOP "could not allocate memory"
+ ! allocate space for time series simulations
+ DO k=1,in%npts
+ gps(k)%nepochs=0
+ ALLOCATE(gps(k)%t(1+ITERATION_MAX), &
+ gps(k)%u1(1+ITERATION_MAX), &
+ gps(k)%u2(1+ITERATION_MAX), &
+ gps(k)%u3(1+ITERATION_MAX), STAT=iostatus)
+ IF (iostatus>0) STOP "could not allocate prediction array"
+ END DO
! allocate memory
- ALLOCATE (v1(1,1,1),v2(1,1,1),v3(1,1,1), &
- u1(1,1,1),u2(1,1,1),u3(1,1,1), &
- tau(in%sx1,in%sx2,in%sx3/2),sig(1,1,1),gamma(1,1,1), &
+ ALLOCATE (v1(in%sx1+2,in%sx2,in%sx3),v2(in%sx1+2,in%sx2,in%sx3),v3(in%sx1+2,in%sx2,in%sx3), &
+ u1(in%sx1+2,in%sx2,in%sx3/2),u2(in%sx1+2,in%sx2,in%sx3/2),u3(in%sx1+2,in%sx2,in%sx3/2), &
+ tau(in%sx1,in%sx2,in%sx3/2),sig(in%sx1,in%sx2,in%sx3/2),gamma(in%sx1+2,in%sx2,in%sx3/2), &
t1(in%sx1+2,in%sx2),t2(in%sx1+2,in%sx2),t3(in%sx1+2,in%sx2), &
STAT=iostatus)
IF (iostatus>0) STOP "could not allocate memory"
- CALL curesetvectors ()
-
- iTensor=3
- CALL cutensormemset (%VAL(iTensor))
+ v1=0;v2=0;v3=0;u1=0;u2=0;u3=0;gamma=0;t1=0;t2=0;t3=0
+ i=0
+ CALL tensorfieldadd(tau,tau,in%sx1,in%sx2,in%sx3/2,c1=0._4,c2=0._4)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! - construct pre-stress structure
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IF (ALLOCATED(in%stresslayer)) THEN
CALL tensorstructure(in%stressstruc,in%stresslayer,in%dx3)
- DEALLOCATE(in%stresslayer)
DO k=1,in%sx3/2
tau(:,:,k)=(-1._4) .times. in%stressstruc(k)%t
END DO
- DEALLOCATE(in%stressstruc)
- END IF
-
- ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- ! - construct linear viscoelastic structure
- ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- IF (ALLOCATED(in%linearlayer)) THEN
- CALL viscoelasticstructure(in%linearstruc,in%linearlayer,in%dx3)
- DEALLOCATE(in%linearlayer)
- END IF
-
- ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- ! - construct nonlinear viscoelastic structure
- ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- IF (ALLOCATED(in%nonlinearlayer)) THEN
- CALL viscoelasticstructure(in%nonlinearstruc,in%nonlinearlayer,in%dx3)
- DEALLOCATE(in%nonlinearlayer)
- END IF
-
- ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- ! - construct nonlinear fault creep structure (rate-strenghtening)
- ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- IF (ALLOCATED(in%faultcreeplayer)) THEN
- CALL viscoelasticstructure(in%faultcreepstruc,in%faultcreeplayer,in%dx3)
- DEALLOCATE(in%faultcreeplayer)
END IF
! first event
@@ -306,17 +280,12 @@ PROGRAM relax
! sources
CALL dislocations(in%events(e),in%lambda,in%mu,in%beta,in%sx1,in%sx2,in%sx3, &
in%dx1,in%dx2,in%dx3,v1,v2,v3,t1,t2,t3,tau)
-
- iTensor = 1
- CALL copytau (tau, %VAL(in%sx1), %VAL(in%sx2), %VAL(in%sx3/2), %VAL(iTensor))
-
- CALL cucopytraction (t3, %VAL(in%sx1), %VAL(in%sx2), iRet)
-
CALL traction(in%mu,in%events(e),in%sx1,in%sx2,in%dx1,in%dx2,t,0.d0,t3)
-
- PRINT '("coseismic event ",I3.3)', e
- PRINT 0990
+ IF (isverbose) THEN
+ PRINT '("coseismic event ",I3.3)', e
+ PRINT 0990
+ END IF
! test the presence of dislocations for coseismic calculation
IF ((in%events(e)%nt .NE. 0) .OR. &
@@ -330,65 +299,90 @@ PROGRAM relax
! add displacement from analytic solutions for small patches (avoid
! aliasing)
+ CALL dislocations_disp(in%events(e),in%lambda,in%mu, &
+ in%sx1,in%sx2,in%sx3, &
+ in%dx1,in%dx2,in%dx3,v1,v2,v3)
END IF
! transfer solution
- CALL cufieldrep (%VAL(in%sx1+2),%VAL(in%sx2),%VAL(in%sx3/2))
-
+ CALL fieldrep(u1,v1,in%sx1+2,in%sx2,in%sx3/2)
+ CALL fieldrep(u2,v2,in%sx1+2,in%sx2,in%sx3/2)
+ CALL fieldrep(u3,v3,in%sx1+2,in%sx2,in%sx3/2)
! evaluate stress
- iTensor=2
- cuc1=0._4
- cuc2=-1._4
- CALL cutensorfieldadd (%VAL(iTensor), %VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2), &
- %VAL(cuc1), %VAL(cuc2), sig, tau)
+ CALL tensorfieldadd(sig,tau,in%sx1,in%sx2,in%sx3/2,c1=0._4,c2=-1._4)
+ CALL stressupdate(u1,u2,u3,in%lambda,in%mu, &
+ in%dx1,in%dx2,in%dx3,in%sx1,in%sx2,in%sx3/2,sig)
+
+ ! export the initial displacements
+ CALL pts2series(u1,u2,u3,in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3,in%opts,0._8,1,gps)
+
+ IF (isverbose) THEN
+ IF (ALLOCATED(in%ptsname)) THEN
+ CALL exportpoints(u1,u2,u3,sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
+ in%opts,in%ptsname,0._8,in%wdir,.true.,in%x0,in%y0,in%rot)
+ END IF
+ END IF
- iTensor = 1
- CALL custressupdatewrapper (%VAL(iTensor), %VAL(in%lambda), %VAL(in%mu), &
- %VAL(in%dx1), %VAL(in%dx2), %VAL(in%dx3), &
- %VAL(in%sx1), %VAL(in%sx2), %VAL(in%sx3/2), u1, u2, u3, sig)
+ WRITE (digit4,'(I4.4)') 0
+
+ IF (isverbose) THEN
+ PRINT 1101,0,0._8,0._8,0._8,0._8,0._8,in%interval,0._8,tensoramplitude(tau,in%dx1,in%dx2,in%dx3)
+ END IF
+ IF (in%interval .LE. 0) THEN
+ GOTO 100 ! no time integration
+ END IF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- ! - export observation points, map view of displacement,
- ! - map view of stress components, Coulomb stress on observation
- ! - patches, and full displacement and stress field.
+ ! - construct linear viscoelastic structure
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ IF (ALLOCATED(in%linearlayer)) THEN
+ CALL viscoelasticstructure(in%linearstruc,in%linearlayer,in%dx3)
+ DEALLOCATE(in%linearlayer)
-#ifdef GRD
- IF (in%isoutputgrd) THEN
- CALL exportgrd(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,in%x0,in%y0,in%wdir,0)
+ IF (0 .LT. in%nlwz) THEN
+ ALLOCATE(lineardgammadot0(in%sx1,in%sx2,in%sx3/2),STAT=iostatus)
+ IF (iostatus.GT.0) STOP "could not allocate lineardgammadot0"
+ CALL builddgammadot0(in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,0._8, &
+ in%nlwz,in%linearweakzone,lineardgammadot0)
+ END IF
END IF
-#endif
- IF (ALLOCATED(in%ptsname)) THEN
- CALL exportpoints(u1,u2,u3,sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
- in%opts,in%ptsname,0._8,in%wdir,.true.,in%x0,in%y0,in%rot)
- END IF
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ ! - construct nonlinear viscoelastic structure
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ IF (ALLOCATED(in%nonlinearlayer)) THEN
+ CALL viscoelasticstructure(in%nonlinearstruc,in%nonlinearlayer,in%dx3)
+ DEALLOCATE(in%nonlinearlayer)
- iTensor=2
- CALL cutensoramp (%VAL(iTensor),%VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2),dAmp)
- dAmp=dAmp*DBLE(in%dx1*in%dx2*in%dx3)
- PRINT 1101,0,0._8,0._8,0._8,0._8,0._8,in%interval,0._8,dAmp
+ IF (0 .LT. in%nnlwz) THEN
+ ALLOCATE(nonlineardgammadot0(in%sx1,in%sx2,in%sx3/2),STAT=iostatus)
+ IF (iostatus.GT.0) STOP "could not allocate nonlineardgammadot0"
+ CALL builddgammadot0(in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,0._8, &
+ in%nnlwz,in%nonlinearweakzone,nonlineardgammadot0)
+ END IF
+ END IF
- IF (in%interval .LE. 0) THEN
- GOTO 100 ! no time integration
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ ! - construct nonlinear fault creep structure (rate-strenghtening)
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ IF (ALLOCATED(in%faultcreeplayer)) THEN
+ CALL viscoelasticstructure(in%faultcreepstruc,in%faultcreeplayer,in%dx3)
+ DEALLOCATE(in%faultcreeplayer)
END IF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! - start the relaxation
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- iTensor=2
- CALL cutensormemset (%VAL(iTensor))
+ ALLOCATE(moment(in%sx1,in%sx2,in%sx3/2),STAT=iostatus)
+ IF (iostatus>0) STOP "could not allocate the mechanical structure"
+
+ CALL tensorfieldadd(moment,moment,in%sx1,in%sx2,in%sx3/2,c1=0._4,c2=0._4)
DO i=1,ITERATION_MAX
IF (t .GE. in%interval) GOTO 100 ! proper exit
-#ifdef PAPI_PROF
- ! start timer
- CALL papistartprofiling (cTimerName)
-#endif
-
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! predictor
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -401,31 +395,36 @@ PROGRAM relax
mech(:)=0
! initialize no forcing term in tensor space
- iTensor=2
- CALL cutensormemset (%VAL(iTensor))
+ CALL tensorfieldadd(moment,moment,in%sx1,in%sx2,in%sx3/2,0._4,0._4)
! power density from three mechanisms (linear and power-law viscosity
! and fault creep)
! 1- linear viscosity
-
-#ifdef PAPI_PROF
- cTimerName = 'Eigenstress'
- CALL papistartprofiling (cTimerName)
-#endif
-
IF (ALLOCATED(in%linearstruc)) THEN
- CALL viscouseigenstress(in%mu,in%linearstruc,in%linearweakzone,in%nlwz, &
- sig,in%sx1,in%sx2,in%sx3/2, &
- in%dx1,in%dx2,in%dx3,moment,0.01_8,MAXWELLTIME=maxwell(1))
+ IF (0 .LT. in%nlwz) THEN
+ CALL viscouseigenstress(in%mu,in%linearstruc, &
+ sig,in%sx1,in%sx2,in%sx3/2, &
+ in%dx1,in%dx2,in%dx3,moment,DGAMMADOT0=lineardgammadot0,MAXWELLTIME=maxwell(1))
+ ELSE
+ CALL viscouseigenstress(in%mu,in%linearstruc, &
+ sig,in%sx1,in%sx2,in%sx3/2, &
+ in%dx1,in%dx2,in%dx3,moment,MAXWELLTIME=maxwell(1))
+ END IF
mech(1)=1
END IF
! 2- powerlaw viscosity
IF (ALLOCATED(in%nonlinearstruc)) THEN
- CALL viscouseigenstress(in%mu,in%nonlinearstruc,in%nonlinearweakzone,in%nnlwz, &
- sig,in%sx1,in%sx2,in%sx3/2, &
- in%dx1,in%dx2,in%dx3,moment,0.01_8,MAXWELLTIME=maxwell(2))
- mech(2)=1
+ IF (0 .LT. in%nnlwz) THEN
+ CALL viscouseigenstress(in%mu,in%nonlinearstruc, &
+ sig,in%sx1,in%sx2,in%sx3/2, &
+ in%dx1,in%dx2,in%dx3,moment,DGAMMADOT0=nonlineardgammadot0,MAXWELLTIME=maxwell(2))
+ ELSE
+ CALL viscouseigenstress(in%mu,in%nonlinearstruc, &
+ sig,in%sx1,in%sx2,in%sx3/2, &
+ in%dx1,in%dx2,in%dx3,moment,MAXWELLTIME=maxwell(2))
+ END IF
+ mech(2)=1
END IF
! 3- nonlinear fault creep with rate-strengthening friction
@@ -441,10 +440,6 @@ PROGRAM relax
mech(3)=1
END IF
-#ifdef PAPI_PROF
- CALL papiendprofiling (cTimerName)
-#endif
-
! identify the required time step
tm=1._8/(REAL(mech(1))/maxwell(1)+ &
REAL(mech(2))/maxwell(2)+ &
@@ -463,112 +458,65 @@ PROGRAM relax
! choose an integration time step
CALL integrationstep(tm,Dt,t,oi,in%odt,in%skip,in%tscale,in%events,e,in%ne)
- iTensor=4
- cuc1=0._4
- cuc2=1._4
- CALL cutensorfieldadd (%VAL(iTensor), %VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2), %VAL(cuc1), %VAL(cuc2), sig, moment)
-
- CALL curesetvectors ()
-
-#ifdef PAPI_PROF
- cTimerName = 'bodyforce'
- CALL papistartprofiling(cTimerName)
-#endif
-
- iTensor=1
- CALL cubodyforceswrapper (%VAL(iTensor), %VAL(in%dx1), %VAL(in%dx2), %VAL(in%dx3), &
- %VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2), v1,v2,v3,t1,t2,t3)
-
-#ifdef PAPI_PROF
- CALL papiendprofiling (cTimerName)
-#endif
+ CALL tensorfieldadd(sig,moment,in%sx1,in%sx2,in%sx3/2,c1=0.0_4,c2=1._4)
-
- CALL cucopytraction (t3, %VAL(in%sx1), %VAL(in%sx2), iRet)
+ v1=0;v2=0;v3=0;t1=0;t2=0;t3=0;
+ CALL equivalentbodyforce(sig,in%dx1,in%dx2,in%dx3,in%sx1,in%sx2,in%sx3/2,v1,v2,v3,t1,t2,t3)
! add time-dependent surface loads
CALL traction(in%mu,in%events(e),in%sx1,in%sx2,in%dx1,in%dx2,t,Dt/2.d8,t3,rate=.TRUE.)
-#ifdef PAPI_PROF
- cTimerName = 'green'
- CALL papistartprofiling(cTimerName)
-#endif
-
CALL greenfunctioncowling(v1,v2,v3,t1,t2,t3,in%dx1,in%dx2,in%dx3,in%lambda,in%mu,in%gam)
-#ifdef PAPI_PROF
- CALL papiendprofiling (cTimerName)
-#endif
-
-#ifdef PAPI_PROF
- cTimerName = 'fieldadd'
- CALL papistartprofiling(cTimerName)
-#endif
-
! v1,v2,v3 contain the predictor displacement
- iTensor = 2
- cuC1 = REAL(Dt/2)
- cuC2 = 1._4
- CALL cufieldadd (%VAL(iTensor), v1, v2, v3, u1, u2, u3, %VAL(in%sx1+2),%VAL(in%sx2),%VAL(in%sx3/2), %VAL(cuC1), %VAL(cuC2))
-
-#ifdef PAPI_PROF
- CALL papiendprofiling(cTimerName)
-#endif
-
- iTensor=2
- cuc1=-REAL(Dt/2)
- cuc2=-1._4
- CALL cutensorfieldadd (%VAL(iTensor), %VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2), %VAL(cuc1), %VAL(cuc2), sig, tau)
-
-#ifdef PAPI_PROF
- cTimerName = 'stress'
- CALL papistartprofiling(cTimerName)
-#endif
-
- ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- ! corrector
- ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- iTensor = 2
- CALL custressupdatewrapper (%VAL(iTensor), %VAL(in%lambda), %VAL(in%mu), &
- %VAL(in%dx1), %VAL(in%dx2), %VAL(in%dx3), &
- %VAL(in%sx1), %VAL(in%sx2), %VAL(in%sx3/2), v1, v2, v3, sig)
+ CALL fieldadd(v1,u1,in%sx1+2,in%sx2,in%sx3/2,c1=REAL(Dt/2))
+ CALL fieldadd(v2,u2,in%sx1+2,in%sx2,in%sx3/2,c1=REAL(Dt/2))
+ CALL fieldadd(v3,u3,in%sx1+2,in%sx2,in%sx3/2,c1=REAL(Dt/2))
+ CALL tensorfieldadd(sig,tau,in%sx1,in%sx2,in%sx3/2,c1=-REAL(Dt/2),c2=-1._4)
-#ifdef PAPI_PROF
- CALL papiendprofiling (cTimerName)
-#endif
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ ! corrector
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ CALL stressupdate(v1,v2,v3,in%lambda,in%mu, &
+ in%dx1,in%dx2,in%dx3,in%sx1,in%sx2,in%sx3/2,sig)
! reinitialize moment density tensor
- iTensor=2
- CALL cutensormemset (%VAL(iTensor))
+ CALL tensorfieldadd(moment,moment,in%sx1,in%sx2,in%sx3/2,0._4,0._4)
IF (ALLOCATED(in%linearstruc)) THEN
! linear viscosity
-! v1=0
- CALL viscouseigenstress(in%mu,in%linearstruc,in%linearweakzone,in%nlwz,sig, &
- in%sx1,in%sx2,in%sx3/2, &
- in%dx1,in%dx2,in%dx3,moment,0.01_8,GAMMA=v1)
+ v1=0
+ IF (0 .LT. in%nlwz) THEN
+ CALL viscouseigenstress(in%mu,in%linearstruc, &
+ sig,in%sx1,in%sx2,in%sx3/2, &
+ in%dx1,in%dx2,in%dx3,moment,DGAMMADOT0=lineardgammadot0,GAMMA=v1)
+ ELSE
+ CALL viscouseigenstress(in%mu,in%linearstruc, &
+ sig,in%sx1,in%sx2,in%sx3/2, &
+ in%dx1,in%dx2,in%dx3,moment,GAMMA=v1)
+ END IF
+ ! update slip history
+ CALL fieldadd(gamma,v1,in%sx1+2,in%sx2,in%sx3/2,c2=REAL(Dt))
END IF
-#ifdef PAPI_PROF
- cTimerName = 'Eigenstress'
- CALL papistartprofiling(cTimerName)
-#endif
-
IF (ALLOCATED(in%nonlinearstruc)) THEN
! powerlaw viscosity
-! v1=0
- CALL viscouseigenstress(in%mu,in%nonlinearstruc,in%nonlinearweakzone,in%nnlwz,sig, &
- in%sx1,in%sx2,in%sx3/2, &
- in%dx1,in%dx2,in%dx3,moment,0.01_8,GAMMA=v1)
+ v1=0
+ IF (0 .LT. in%nnlwz) THEN
+ CALL viscouseigenstress(in%mu,in%nonlinearstruc, &
+ sig,in%sx1,in%sx2,in%sx3/2, &
+ in%dx1,in%dx2,in%dx3,moment,DGAMMADOT0=nonlineardgammadot0,GAMMA=v1)
+ ELSE
+ CALL viscouseigenstress(in%mu,in%nonlinearstruc, &
+ sig,in%sx1,in%sx2,in%sx3/2, &
+ in%dx1,in%dx2,in%dx3,moment,GAMMA=v1)
+ END IF
+ ! update slip history
+ CALL fieldadd(gamma,v1,in%sx1+2,in%sx2,in%sx3/2,c2=REAL(Dt))
END IF
-#ifdef PAPI_PROF
- CALL papiendprofiling (cTimerName)
-#endif
-
-
! nonlinear fault creep with rate-strengthening friction
IF (ALLOCATED(in%faultcreepstruc)) THEN
@@ -599,56 +547,19 @@ PROGRAM relax
in%dx1,in%dx2,in%dx3,v1,v2,v3,t1,t2,t3,tau,eigenstress=moment)
END IF
+ v1=0;v2=0;v3=0;t1=0;t2=0;t3=0;
+ CALL equivalentbodyforce(moment,in%dx1,in%dx2,in%dx3,in%sx1,in%sx2,in%sx3/2,v1,v2,v3,t1,t2,t3)
- CALL curesetvectors ()
-
-#ifdef PAPI_PROF
- cTimerName = 'bodyforce'
- CALL papistartprofiling(cTimerName)
-#endif
-
- iTensor=2
- CALL cubodyforceswrapper (%VAL(iTensor), %VAL(in%dx1), %VAL(in%dx2), %VAL(in%dx3), &
- %VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2), v1,v2,v3,t1,t2,t3)
-
-#ifdef PAPI_PROF
- CALL papiendprofiling (cTimerName)
-#endif
-
-
- CALL cucopytraction (t3, %VAL(in%sx1), %VAL(in%sx2), iRet)
-
-#ifdef PAPI_PROF
- cTimerName = 'green'
- CALL papistartprofiling(cTimerName)
-#endif
+ ! add time-dependent surface loads
+ CALL traction(in%mu,in%events(e),in%sx1,in%sx2,in%dx1,in%dx2,t,Dt,t3,rate=.TRUE.)
CALL greenfunctioncowling(v1,v2,v3,t1,t2,t3,in%dx1,in%dx2,in%dx3,in%lambda,in%mu,in%gam)
-#ifdef PAPI_PROF
- CALL papiendprofiling (cTimerName)
-#endif
-
-#ifdef PAPI_PROF
- cTimerName = 'fieldadd'
- CALL papistartprofiling(cTimerName)
-#endif
-
- iTensor = 1
- cuC1 = 1._4
- cuC2 = REAL(Dt)
- CALL cufieldadd (%VAL(iTensor), u1, u2, u3, v1, v2, v3, %VAL(in%sx1+2),%VAL(in%sx2),%VAL(in%sx3/2), %VAL(cuC1), %VAL(cuC2))
-
-#ifdef PAPI_PROF
- CALL papiendprofiling(cTimerName)
-#endif
-
-
- iTensor=5
- cuc1=1._4
- cuc2=REAL(Dt)
- CALL cutensorfieldadd (%VAL(iTensor), %VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2), %VAL(cuc1), %VAL(cuc2), tau, moment)
-
+ ! update deformation field
+ CALL fieldadd(u1,v1,in%sx1+2,in%sx2,in%sx3/2,c2=REAL(Dt))
+ CALL fieldadd(u2,v2,in%sx1+2,in%sx2,in%sx3/2,c2=REAL(Dt))
+ CALL fieldadd(u3,v3,in%sx1+2,in%sx2,in%sx3/2,c2=REAL(Dt))
+ CALL tensorfieldadd(tau,moment,in%sx1,in%sx2,in%sx3/2,c2=REAL(Dt))
CALL frictionadd(in%np,in%n,Dt)
! time increment
@@ -660,114 +571,88 @@ PROGRAM relax
e=e+1
in%events(e)%i=i
- PRINT '("coseismic event ",I3.3)', e
- PRINT 0990
-
- CALL curesetvectors ()
- iTensor = 0
- CALL copytau (tau, %VAL(in%sx1), %VAL(in%sx2), %VAL(in%sx3/2), %VAL(iTensor))
+ IF (isverbose) THEN
+ PRINT '("coseismic event ",I3.3)', e
+ PRINT 0990
+ END IF
+ v1=0;v2=0;v3=0;t1=0;t2=0;t3=0;
CALL dislocations(in%events(e),in%lambda,in%mu, &
in%beta,in%sx1,in%sx2,in%sx3, &
in%dx1,in%dx2,in%dx3,v1,v2,v3,t1,t2,t3,tau)
-
- iTensor = 1
- CALL copytau (tau, %VAL(in%sx1), %VAL(in%sx2), %VAL(in%sx3/2), %VAL(iTensor))
- CALL cucopytraction (t3, %VAL(in%sx1), %VAL(in%sx2), iRet)
-
CALL traction(in%mu,in%events(e),in%sx1,in%sx2,in%dx1,in%dx2,t,0.d0,t3)
! apply the 3d elastic transfert function
CALL greenfunctioncowling(v1,v2,v3,t1,t2,t3, &
in%dx1,in%dx2,in%dx3,in%lambda,in%mu,in%gam)
- iTensor = 1
- cuC1 = 1._4
- cuC2 = 1._4
- CALL cufieldadd (%VAL(iTensor), u1, u2, u3, v1, v2, v3, %VAL(in%sx1+2),%VAL(in%sx2),%VAL(in%sx3/2), %VAL(cuC1), %VAL(cuC2))
+ ! transfer solution
+ CALL fieldadd(u1,v1,in%sx1+2,in%sx2,in%sx3/2)
+ CALL fieldadd(u2,v2,in%sx1+2,in%sx2,in%sx3/2)
+ CALL fieldadd(u3,v3,in%sx1+2,in%sx2,in%sx3/2)
+
END IF
END IF
- iTensor=2
- cuc1=0._4
- cuc2=-1._4
- CALL cutensorfieldadd (%VAL(iTensor), %VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2), %VAL(cuc1), %VAL(cuc2), sig, tau)
+ CALL tensorfieldadd(sig,tau,in%sx1,in%sx2,in%sx3/2,c1=0._4,c2=-1._4)
+ CALL stressupdate(u1,u2,u3,in%lambda,in%mu, &
+ in%dx1,in%dx2,in%dx3,in%sx1,in%sx2,in%sx3/2,sig)
+ ! export displacements
+ CALL pts2series(u1,u2,u3,in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3,in%opts,t,1+i,gps)
-#ifdef PAPI_PROF
- cTimerName = 'stress'
- CALL papistartprofiling(cTimerName)
-#endif
+ IF (isverbose) THEN
+ IF (ALLOCATED(in%ptsname)) THEN
+ CALL exportpoints(u1,u2,u3,sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
+ in%opts,in%ptsname,t,in%wdir,.FALSE.,in%x0,in%y0,in%rot)
+ END IF
+ END IF
- iTensor = 1
- CALL custressupdatewrapper (%VAL(iTensor), %VAL(in%lambda), %VAL(in%mu), &
- %VAL(in%dx1), %VAL(in%dx2), %VAL(in%dx3), &
- %VAL(in%sx1), %VAL(in%sx2), %VAL(in%sx3/2), u1, u2, u3, sig)
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ ! - export displacement and stress
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-#ifdef PAPI_PROF
- CALL papiendprofiling(cTimerName)
-#endif
- ! points are exported at all time steps
- IF (ALLOCATED(in%ptsname)) THEN
- CALL exportpoints(u1,u2,u3,sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
- in%opts,in%ptsname,t,in%wdir,.FALSE.,in%x0,in%y0,in%rot)
+ ! output only at discrete intervals (skip=0, odt>0),
+ ! or every "skip" computational steps (skip>0, odt<0),
+ ! or anytime a coseismic event occurs
+ IF (isverbose) THEN
+ IF (isoutput(in%skip,t,i,in%odt,oi,in%events(e)%time)) THEN
+
+ WRITE (digit4,'(I4.4)') oi
+ PRINT 1101,i,Dt,maxwell,t,in%interval, &
+ tensoramplitude(moment,in%dx1,in%dx2,in%dx3), &
+ tensoramplitude(tau,in%dx1,in%dx2,in%dx3)
+
+ ! update output counter
+ oi=oi+1
+ ELSE
+ PRINT 1100,i,Dt,maxwell,t,in%interval, &
+ tensoramplitude(moment,in%dx1,in%dx2,in%dx3), &
+ tensoramplitude(tau,in%dx1,in%dx2,in%dx3)
+ END IF
END IF
- iTensor=1
- CALL cutensoramp (%VAL(iTensor),%VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2),dAmp, moment)
- iTensor=2
- CALL cutensoramp (%VAL(iTensor),%VAL(in%sx1),%VAL(in%sx2),%VAL(in%sx3/2),dAmp1, tau)
- dAmp=dAmp*DBLE(in%dx1*in%dx2*in%dx3)
- dAmp1=dAmp1*DBLE(in%dx1*in%dx2*in%dx3)
- PRINT 1101,i,Dt,maxwell,t,in%interval, dAmp, dAmp1
- ! update output counter
- oi=oi+1
-
-#ifdef PAPI_PROF
- cTimerName = 'relax'
- CALL papiendprofiling (cTimerName)
-#endif
-
END DO
100 CONTINUE
- DO i=1,in%ne
- IF (ALLOCATED(in%events(i)%s)) DEALLOCATE(in%events(i)%s,in%events(i)%sc)
- IF (ALLOCATED(in%events(i)%ts)) DEALLOCATE(in%events(i)%ts,in%events(i)%tsc)
- END DO
- IF (ALLOCATED(in%events)) DEALLOCATE(in%events)
-
- CALL cudeinit()
+! DO i=1,in%ne
+! IF (ALLOCATED(in%events(i)%s)) DEALLOCATE(in%events(i)%s,in%events(i)%sc)
+! IF (ALLOCATED(in%events(i)%ts)) DEALLOCATE(in%events(i)%ts,in%events(i)%tsc)
+! END DO
+! IF (ALLOCATED(in%events)) DEALLOCATE(in%events)
! free memory
IF (ALLOCATED(gamma)) DEALLOCATE(gamma)
- IF (ALLOCATED(in%opts)) DEALLOCATE(in%opts)
- IF (ALLOCATED(in%ptsname)) DEALLOCATE(in%ptsname)
- IF (ALLOCATED(in%op)) DEALLOCATE(in%op)
- IF (ALLOCATED(in%sop)) DEALLOCATE(in%sop)
- IF (ALLOCATED(in%n)) DEALLOCATE(in%n)
- IF (ALLOCATED(in%stressstruc)) DEALLOCATE(in%stressstruc)
- IF (ALLOCATED(in%stresslayer)) DEALLOCATE(in%stresslayer)
- IF (ALLOCATED(in%linearstruc)) DEALLOCATE(in%linearstruc)
- IF (ALLOCATED(in%linearlayer)) DEALLOCATE(in%linearlayer)
- IF (ALLOCATED(in%linearweakzone)) DEALLOCATE(in%linearweakzone)
- IF (ALLOCATED(in%nonlinearstruc)) DEALLOCATE(in%nonlinearstruc)
- IF (ALLOCATED(in%nonlinearlayer)) DEALLOCATE(in%nonlinearlayer)
- IF (ALLOCATED(in%nonlinearweakzone)) DEALLOCATE(in%nonlinearweakzone)
- IF (ALLOCATED(in%faultcreepstruc)) DEALLOCATE(in%faultcreepstruc)
- IF (ALLOCATED(in%faultcreeplayer)) DEALLOCATE(in%faultcreeplayer)
IF (ALLOCATED(sig)) DEALLOCATE(sig)
IF (ALLOCATED(tau)) DEALLOCATE(tau)
IF (ALLOCATED(moment)) DEALLOCATE(moment)
- IF (ALLOCATED(in%stresslayer)) DEALLOCATE(in%stresslayer)
- IF (ALLOCATED(in%linearlayer)) DEALLOCATE(in%linearlayer)
- IF (ALLOCATED(in%nonlinearlayer)) DEALLOCATE(in%nonlinearlayer)
- IF (ALLOCATED(in%faultcreeplayer)) DEALLOCATE(in%faultcreeplayer)
IF (ALLOCATED(v1)) DEALLOCATE(v1,v2,v3,t1,t2,t3)
IF (ALLOCATED(u1)) DEALLOCATE(u1,u2,u3)
- IF (ALLOCATED(inter1)) DEALLOCATE(inter1,inter2,inter3)
+#ifdef FFTW3_THREADS
+ CALL sfftw_cleanup_threads()
+#endif
0990 FORMAT (" I | Dt | tm(ve) | tm(pl) | tm(as) | t/tmax | power | C:E^i | ")
1100 FORMAT (I3.3," ",ES9.2E2,3ES9.2E2,ES9.2E2,"/",ES7.2E1,2ES9.2E2)
@@ -1066,5 +951,4 @@ CONTAINS
END SUBROUTINE integrationstep
-
-END PROGRAM relax
+END SUBROUTINE relaxlite
diff --git a/src/types.f90 b/src/types.f90
index 8c982dc..b733b03 100644
--- a/src/types.f90
+++ b/src/types.f90
@@ -20,6 +20,7 @@
#include "include.f90"
MODULE types
+ USE ISO_C_BINDING
TYPE SOURCE_STRUCT
SEQUENCE
@@ -33,7 +34,7 @@ MODULE types
TYPE WEAK_STRUCT
SEQUENCE
- REAL*8 :: dgammadot0,x,y,z,width,length,thickness,strike,dip
+ REAL*8 :: dgammadot0,x,y,z,width,length,thickness,strike,dip,age
END TYPE WEAK_STRUCT
TYPE VECTOR_STRUCT
@@ -97,7 +98,7 @@ MODULE types
REAL*8, DIMENSION(:), ALLOCATABLE :: s1,s2,s3
END TYPE MANIFOLD_STRUCT
- TYPE, PUBLIC :: SIMULATION_STRUC
+ TYPE, PUBLIC :: SIMULATION_STRUCT
! grid dimension
INTEGER :: sx1,sx2,sx3
@@ -225,6 +226,6 @@ MODULE types
TYPE(MANIFOLD_STRUCT), DIMENSION(:), ALLOCATABLE :: gps,sim
- END TYPE SIMULATION_STRUC
+ END TYPE SIMULATION_STRUCT
END MODULE types
diff --git a/wscript b/wscript
index 3c2b1e6..8504cad 100644
--- a/wscript
+++ b/wscript
@@ -135,7 +135,7 @@ def configure(cnf):
if not cnf.options.cuda_incdir:
cnf.options.cuda_incdir=cnf.options.cuda_dir + "/include"
if not cnf.options.cuda_libdir:
- cnf.options.cuda_libdir=cnf.options.cuda_dir + "/lib64"
+ cnf.options.cuda_libdir=cnf.options.cuda_dir + "/lib"
if cnf.options.cuda_incdir:
includedirs=[cnf.options.cuda_incdir]
else:
@@ -280,6 +280,40 @@ def configure(cnf):
cnf.write_config_header('config.h')
+
+def relaxlite(ctx) :
+ ctx.shlib(features='c fc fcprogram',
+ source=['src/relaxlite.f90',
+ 'src/types.f90',
+ 'src/ctfft.f',
+ 'src/fourier.f90',
+ 'src/green.f90',
+ 'src/okada/green_space.f90',
+ 'src/okada/dc3d.f',
+ 'src/elastic3d.f90',
+ 'src/friction3d.f90',
+ 'src/viscoelastic3d.f90',
+ 'src/writevtk.c',
+ 'src/writegrd4.2.c',
+ 'src/proj.c',
+ 'src/export.f90',
+ 'src/getdata.f',
+ 'src/getopt_m.f90',
+ 'src/input.f90',
+ 'src/mkl_dfti.f90',
+ 'src/papi_prof.c'],
+ install_path='${PREFIX}/bin',
+ includes=['build'],
+ use=['gmt','proj','openmp','fftw','imkl','cpp','length','papi','stdc++'],
+ target='relaxlite.so'
+ )
+
+from waflib.Build import BuildContext
+
+class miracle(BuildContext):
+ cmd = 'relaxlite'
+ fun = 'relaxlite'
+
def build(bld):
if bld.env.CUDA:
bld.program(features='c fc fcprogram cxx',
@@ -301,9 +335,9 @@ def build(bld):
'src/getopt_m.f90',
'src/input.f90',
'src/mkl_dfti.f90',
- 'src/papi_prof.c',
- 'src/cu_fft.cu',
- 'src/cu_elastic.cu'],
+ 'src/papi_prof.c',
+ 'src/cugreen.cu',
+ 'src/cuelastic.cu'],
install_path='${PREFIX}/bin',
includes=['build'],
use=['gmt','proj','openmp','fftw','imkl','cpp','length','cuda','papi','stdc++'],
More information about the CIG-COMMITS
mailing list