[cig-commits] commit: fixes bugs w/ gfortran, include c function writevtk.c
Mercurial
hg at geodynamics.org
Mon Oct 31 15:55:17 PDT 2011
changeset: 46:8ee30887f819
tag: tip
user: Sylvain Barbot <sbarbot at caltech.edu>
date: Mon Oct 31 15:55:10 2011 -0700
files: INSTALL examples/run1.sh export.f90 fourier.f90 include.f90 relax.f90 viscoelastic3d.f90 waf wscript
description:
fixes bugs w/ gfortran, include c function writevtk.c
diff -r 9df56652aa63 -r 8ee30887f819 INSTALL
--- a/INSTALL Wed Oct 26 13:37:50 2011 -0700
+++ b/INSTALL Mon Oct 31 15:55:10 2011 -0700
@@ -19,11 +19,11 @@ and then recompile the executable
and then recompile the executable
cd build
- /sw/bin/gfortran relax.f90.1.o types.f90.1.o ctfft.f.1.o fourier.f90.1.o green.f90.1.o elastic3d.f90.1.o friction3d.f90.1.o viscoelastic3d.f90.1.o writegrd4.2.c.1.o proj.c.1.o export.f90.1.o getdata.f.1.o getopt_m.f90.1.o input.f90.1.o mkl_dfti.f90.1.o -o /Users/sbarbot/Documents/src/relax/relax -L.. -lgmt -lnetcdf -lproj -lfftw3f -lfftw3f_threads -lgomp -lquadmath -lgfortran -Wl,-search_paths_first -static-libgcc
+ /sw/bin/gfortran relax.f90.1.o types.f90.1.o ctfft.f.1.o fourier.f90.1.o green.f90.1.o elastic3d.f90.1.o friction3d.f90.1.o viscoelastic3d.f90.1.o writevtk.c.1.o writegrd4.2.c.1.o proj.c.1.o export.f90.1.o getdata.f.1.o getopt_m.f90.1.o input.f90.1.o mkl_dfti.f90.1.o -o /Users/sbarbot/Documents/src/relax/relax -L.. -lgmt -lnetcdf -lproj -lfftw3f -lfftw3f_threads -lgomp -lquadmath -lgfortran -Wl,-search_paths_first -static-libgcc
Then the output of "otool -L ../relax" will give
../relax:
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 125.2.11)
-which is a standard system library available on all Macs.
\ No newline at end of file
+which is a standard system library available on all Macs.
diff -r 9df56652aa63 -r 8ee30887f819 examples/run1.sh
--- a/examples/run1.sh Wed Oct 26 13:37:50 2011 -0700
+++ b/examples/run1.sh Mon Oct 31 15:55:10 2011 -0700
@@ -66,7 +66,7 @@ if [ ! -e $WDIR ]; then
mkdir $WDIR
fi
-time ../relax $* <<EOF | tee $WDIR/in.param
+OMP_NUM_THREADS=4 time ../relax $* <<EOF | tee $WDIR/in.param
# use '#' character to include comments in your input file
# grid size (sx1,sx2,sx3)
256 256 256
diff -r 9df56652aa63 -r 8ee30887f819 export.f90
--- a/export.f90 Wed Oct 26 13:37:50 2011 -0700
+++ b/export.f90 Mon Oct 31 15:55:10 2011 -0700
@@ -1878,8 +1878,9 @@ END SUBROUTINE exportcreep
REAL*8, INTENT(IN) :: dx1,dx2,dx3
CHARACTER(80), INTENT(IN) :: vcfilename
- INTEGER :: iostatus,idum,i1,i2
+ INTEGER :: iostatus,idum,i1,i2,i3
CHARACTER :: q
+ INTEGER :: k1,k2,k3
REAL*8 :: x1,x2,x3
! double-quote character
@@ -1903,12 +1904,15 @@ END SUBROUTINE exportcreep
" format=",a,"ascii",a,">")') q,q,q,q,q,q
! write first component values
- DO x3=0,sx3-1,j3
- DO x2=-sx2/2,sx2/2-1,j2
- DO x1=-sx1/2,sx1/2-1,j1
+ DO k3=0,sx3-1,j3
+ x3=REAL(k3,8)
+ DO k2=-sx2/2,sx2/2-1,j2
+ x2=REAL(k2,8)
+ DO k1=-sx1/2,sx1/2-1,j1
+ x1=REAL(k1,8)
CALL shiftedindex(x1,x2,1._8,sx1,sx2,sx3,1._8,1._8,1._8,i1,i2,idum)
- WRITE (15,'(ES12.2)') u1(i1,i2,x3+1)
+ WRITE (15,'(ES12.2)') u1(i1,i2,k3+1)
END DO
END DO
END DO
@@ -1919,12 +1923,15 @@ END SUBROUTINE exportcreep
" format=",a,"ascii",a,">")') q,q,q,q,q,q
! write second component values
- DO x3=0,sx3-1,j3
- DO x2=-sx2/2,sx2/2-1,j2
- DO x1=-sx1/2,sx1/2-1,j1
+ DO k3=0,sx3-1,j3
+ x3=REAL(k3,8)
+ DO k2=-sx2/2,sx2/2-1,j2
+ x2=REAL(k2,8)
+ DO k1=-sx1/2,sx1/2-1,j1
+ x1=REAL(k1,8)
CALL shiftedindex(x1,x2,1._8,sx1,sx2,sx3,1._8,1._8,1._8,i1,i2,idum)
- WRITE (15,'(ES12.2)') u2(i1,i2,x3+1)
+ WRITE (15,'(ES12.2)') u2(i1,i2,k3+1)
END DO
END DO
@@ -1936,12 +1943,15 @@ END SUBROUTINE exportcreep
" format=",a,"ascii",a,">")') q,q,q,q,q,q
! write third component values
- DO x3=0,sx3-1,j3
- DO x2=-sx2/2,sx2/2-1,j2
- DO x1=-sx1/2,sx1/2-1,j1
+ DO k3=0,sx3-1,j3
+ x3=REAL(k3,8)
+ DO k2=-sx2/2,sx2/2-1,j2
+ x2=REAL(k2,8)
+ DO k1=-sx1/2,sx1/2-1,j1
+ x1=REAL(k1,8)
CALL shiftedindex(x1,x2,1._8,sx1,sx2,sx3,1._8,1._8,1._8,i1,i2,idum)
- WRITE (15,'(ES12.2)') u3(i1,i2,x3+1)
+ WRITE (15,'(ES12.2)') u3(i1,i2,k3+1)
END DO
END DO
@@ -1957,7 +1967,8 @@ END SUBROUTINE exportcreep
" format=",a,"ascii",a, &
" RangeMin=",a,ES12.2,a, &
" RangeMax=",a,ES12.2,a,">")') q,q,q,q,q,q,q,-sx1/2*dx1,q,q,(sx1/2-1)*dx1,q
- DO x1=-sx1/2,sx1/2-1,j1
+ DO k1=-sx1/2,sx1/2-1,j1
+ x1=REAL(k1,8)
WRITE (15,'(ES12.2)') x1*dx1
END DO
WRITE (15,'(" </DataArray>")')
@@ -1966,7 +1977,8 @@ END SUBROUTINE exportcreep
" format=",a,"ascii",a, &
" RangeMin=",a,ES12.2,a, &
" RangeMax=",a,ES12.2,a,">")') q,q,q,q,q,q,q,-sx2/2*dx2,q,q,(sx2/2-1)*dx2,q
- DO x2=-sx2/2,sx2/2-1,j2
+ DO k2=-sx2/2,sx2/2-1,j2
+ x2=REAL(k2,8)
WRITE (15,'(ES12.2)') x2*dx2
END DO
WRITE (15,'(" </DataArray>")')
@@ -1975,7 +1987,8 @@ END SUBROUTINE exportcreep
" format=",a,"ascii",a, &
" RangeMin=",a,ES12.2,a, &
" RangeMax=",a,ES12.2,a,">")') q,q,q,q,q,q,q,0,q,q,(sx3-1)*dx3,q
- DO x3=0,sx3-1,j3
+ DO k3=0,sx3-1,j3
+ x3=REAL(k3,8)
WRITE (15,'(ES12.2)') x3*dx3
END DO
WRITE (15,'(" </DataArray>")')
@@ -1988,150 +2001,6 @@ END SUBROUTINE exportcreep
CLOSE(15)
END SUBROUTINE exportvtk_vectors
-
- !------------------------------------------------------------------
- !> subroutine ExportVTK_Vectors_Legacy
- !! creates a .vtk file in the VTK Legacy binary format with
- !! structured points containing a vector field.
- !!
- !! \author sylvain barbot 05/23/11 - original form
- !------------------------------------------------------------------
- SUBROUTINE exportvtk_vectors_legacy(u1,u2,u3,sx1,sx2,sx3,dx1,dx2,dx3, &
- j1,j2,j3,vcfilename,title,name)
- INTEGER, INTENT(IN) :: sx1,sx2,sx3,j1,j2,j3
-#ifdef ALIGN_DATA
- REAL*4, INTENT(IN), DIMENSION(sx1+2,sx2,sx3) :: u1,u2,u3
-#else
- REAL*4, INTENT(IN), DIMENSION(sx1,sx2,sx3) :: u1,u2,u3
-#endif
- REAL*8, INTENT(IN) :: dx1,dx2,dx3
- CHARACTER(80), INTENT(IN) :: vcfilename,title,name
-
- INTEGER :: iostatus,idum,i1,i2
- CHARACTER :: q,end_rec
- CHARACTER(LEN=180) :: s_buffer
- REAL*8 :: x1,x2,x3
-
- ! double-quote character
- q=char(34)
- ! end-character for binary-record finalize
- end_rec=char(10)
-
- ! open file for mixed binary/ascii writing in VTK legacy format
- OPEN(UNIT=15,FILE=vcfilename,form='UNFORMATTED',ACCESS='SEQUENTIAL', &
- ACTION='WRITE',CONVERT='BIG_ENDIAN', &
- IOSTAT=iostatus)
- IF (iostatus>0) THEN
- WRITE_DEBUG_INFO
- PRINT '(a)', vcfilename
- STOP "could not open file in mixed binary/ascii for export in legacy format"
- END IF
-
- ! writing header of file (INIT)
- WRITE(UNIT=15,IOSTAT=iostatus) '# vtk DataFile Version 3.0'//end_rec
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(title)//end_rec
- WRITE(UNIT=15,IOSTAT=iostatus) 'BINARY'//end_rec
- WRITE(UNIT=15,IOSTAT=iostatus) 'DATASET '//'STRUCTURED_POINTS'//end_rec
-
- ! structured points grid (GEO)
- WRITE(s_buffer,FMT='(A,3I12)',IOSTAT=iostatus) 'DIMENSIONS ',sx1/j1,sx2/j2,sx3/j3
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
- WRITE(s_buffer,FMT='(A,3E14.6E2)',IOSTAT=iostatus) 'ORIGIN ',-dx1*(sx1/2),-dx2*(sx2/2),0.d0
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
- WRITE(s_buffer,FMT='(A,3E14.6E2)',IOSTAT=iostatus) 'SPACING ',dx1*j1,dx2*j2,dx3*j3
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
-
- ! data header for this grid (DAT)
- WRITE(s_buffer,FMT='(A,I12)',IOSTAT=iostatus) 'POINT_DATA ', (sx1/j1)*(sx2/j2)*(sx3/j3)
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
-
- ! data array (VAR)
- WRITE(UNIT=15,IOSTAT=iostatus) 'VECTORS '//TRIM(name)//' float'//end_rec
- DO x3=0,sx3-1,j3
- DO x2=-sx2/2,sx2/2-1,j2
- DO x1=-sx1/2,sx1/2-1,j1
- CALL shiftedindex(x1,x2,1._8,sx1,sx2,sx3,1._8,1._8,1._8,i1,i2,idum)
- WRITE (UNIT=15,IOSTAT=iostatus) REAL(u1(i1,i2,x3+1)),REAL(u2(i1,i2,x3+1)),REAL(u3(i1,i2,x3+1))
- END DO
- END DO
- END DO
- WRITE(UNIT=15,IOSTAT=iostatus) end_rec
-
- ! close binary file (END)
- CLOSE(15)
-
- END SUBROUTINE exportvtk_vectors_legacy
-
- !------------------------------------------------------------------
- !> subroutine ExportVTK_tensors_Legacy
- !! creates a .vtk file in the VTK Legacy binary format with
- !! structured points containing a vector field.
- !!
- !! \author sylvain barbot 05/23/11 - original form
- !------------------------------------------------------------------
- SUBROUTINE exportvtk_tensors_legacy(sig,sx1,sx2,sx3,dx1,dx2,dx3, &
- j1,j2,j3,vcfilename,title,name)
- INTEGER, INTENT(IN) :: sx1,sx2,sx3,j1,j2,j3
- TYPE(TENSOR), INTENT(IN), DIMENSION(sx1,sx2,sx3) :: sig
- REAL*8, INTENT(IN) :: dx1,dx2,dx3
- CHARACTER(80), INTENT(IN) :: vcfilename,title,name
-
- INTEGER :: iostatus,idum,i1,i2
- CHARACTER :: q,end_rec
- CHARACTER(LEN=180) :: s_buffer
- REAL*8 :: x1,x2,x3
-
- ! double-quote character
- q=char(34)
- ! end-character for binary-record finalize
- end_rec=char(10)
-
- ! open file for mixed binary/ascii writing in VTK legacy format
- OPEN(UNIT=15,FILE=vcfilename,form='UNFORMATTED',ACCESS='SEQUENTIAL', &
- ACTION='WRITE',CONVERT='BIG_ENDIAN', &
- IOSTAT=iostatus)
- IF (iostatus>0) THEN
- WRITE_DEBUG_INFO
- PRINT '(a)', vcfilename
- STOP "could not open file in mixed binary/ascii for export in legacy format"
- END IF
-
- ! writing header of file (INIT)
- WRITE(UNIT=15,IOSTAT=iostatus) '# vtk DataFile Version 3.0'//end_rec
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(title)//end_rec
- WRITE(UNIT=15,IOSTAT=iostatus) 'BINARY'//end_rec
- WRITE(UNIT=15,IOSTAT=iostatus) 'DATASET '//'STRUCTURED_POINTS'//end_rec
-
- ! structured points grid (GEO)
- WRITE(s_buffer,FMT='(A,3I12)',IOSTAT=iostatus) 'DIMENSIONS ',sx1/j1,sx2/j2,sx3/j3
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
- WRITE(s_buffer,FMT='(A,3E14.6E2)',IOSTAT=iostatus) 'ORIGIN ',-dx1*(sx1/2),-dx2*(sx2/2),0.d0
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
- WRITE(s_buffer,FMT='(A,3E14.6E2)',IOSTAT=iostatus) 'SPACING ',dx1*j1,dx2*j2,dx3*j3
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
-
- ! data header for this grid (DAT)
- WRITE(s_buffer,FMT='(A,I12)',IOSTAT=iostatus) 'POINT_DATA ', (sx1/j1)*(sx2/j2)*(sx3/j3)
- WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
-
- ! data array (VAR)
- WRITE(UNIT=15,IOSTAT=iostatus) 'TENSORS '//TRIM(name)//' float'//end_rec
- DO x3=0,sx3-1,j3
- DO x2=-sx2/2,sx2/2-1,j2
- DO x1=-sx1/2,sx1/2-1,j1
- CALL shiftedindex(x1,x2,1._8,sx1,sx2,sx3,1._8,1._8,1._8,i1,i2,idum)
- WRITE (UNIT=15,IOSTAT=iostatus) sig(i1,i2,x3+1)%s11,sig(i1,i2,x3+1)%s12,sig(i1,i2,x3+1)%s13, &
- sig(i1,i2,x3+1)%s12,sig(i1,i2,x3+1)%s22,sig(i1,i2,x3+1)%s23, &
- sig(i1,i2,x3+1)%s13,sig(i1,i2,x3+1)%s23,sig(i1,i2,x3+1)%s33
- END DO
- END DO
- END DO
- WRITE(UNIT=15,IOSTAT=iostatus) end_rec
-
- ! close binary file (END)
- CLOSE(15)
-
- END SUBROUTINE exportvtk_tensors_legacy
!------------------------------------------------------------------
!> subroutine ExportVTK_Vectors_Slice
@@ -2152,6 +2021,7 @@ END SUBROUTINE exportcreep
INTEGER :: iostatus,idum,i1,i2
CHARACTER :: q
+ INTEGER :: k1,k2,k3
REAL*8 :: x1,x2,x3
! double-quote character
@@ -2176,11 +2046,13 @@ END SUBROUTINE exportcreep
! write first component values
x3=oz/dx3
- DO x2=-sx2/2,sx2/2-1,j2
- DO x1=-sx1/2,sx1/2-1,j1
+ DO k2=-sx2/2,sx2/2-1,j2
+ x2=REAL(k2,8)
+ DO k1=-sx1/2,sx1/2-1,j1
+ x1=REAL(k1,8)
CALL shiftedindex(x1,x2,1._8,sx1,sx2,sx3,1._8,1._8,1._8,i1,i2,idum)
- WRITE (15,'(ES12.2)') u1(i1,i2,x3+1)
+ WRITE (15,'(ES12.2)') u1(i1,i2,k3+1)
END DO
END DO
WRITE (15,'(" </DataArray>")')
@@ -2191,11 +2063,13 @@ END SUBROUTINE exportcreep
! write second component values
x3=oz/dx3
- DO x2=-sx2/2,sx2/2-1,j2
- DO x1=-sx1/2,sx1/2-1,j1
+ DO k2=-sx2/2,sx2/2-1,j2
+ x2=REAL(k2,8)
+ DO k1=-sx1/2,sx1/2-1,j1
+ x1=REAL(k1,8)
CALL shiftedindex(x1,x2,1._8,sx1,sx2,sx3,1._8,1._8,1._8,i1,i2,idum)
- WRITE (15,'(ES12.2)') u2(i1,i2,x3+1)
+ WRITE (15,'(ES12.2)') u2(i1,i2,k3+1)
END DO
END DO
@@ -2207,11 +2081,13 @@ END SUBROUTINE exportcreep
! write third component values
x3=oz/dx3
- DO x2=-sx2/2,sx2/2-1,j2
- DO x1=-sx1/2,sx1/2-1,j1
+ DO k2=-sx2/2,sx2/2-1,j2
+ x2=REAL(k2,8)
+ DO k1=-sx1/2,sx1/2-1,j1
+ x1=REAL(k1,8)
CALL shiftedindex(x1,x2,1._8,sx1,sx2,sx3,1._8,1._8,1._8,i1,i2,idum)
- WRITE (15,'(ES12.2)') u3(i1,i2,x3+1)
+ WRITE (15,'(ES12.2)') u3(i1,i2,k3+1)
END DO
END DO
@@ -2226,7 +2102,8 @@ END SUBROUTINE exportcreep
" format=",a,"ascii",a, &
" RangeMin=",a,ES12.2,a, &
" RangeMax=",a,ES12.2,a,">")') q,q,q,q,q,q,q,-sx1/2*dx1,q,q,(sx1/2-1)*dx1,q
- DO x1=-sx1/2,sx1/2-1,j1
+ DO k1=-sx1/2,sx1/2-1,j1
+ x1=REAL(k1,8)
WRITE (15,'(ES12.2)') x1*dx1
END DO
WRITE (15,'(" </DataArray>")')
@@ -2235,7 +2112,8 @@ END SUBROUTINE exportcreep
" format=",a,"ascii",a, &
" RangeMin=",a,ES12.2,a, &
" RangeMax=",a,ES12.2,a,">")') q,q,q,q,q,q,q,-sx2/2*dx1,q,q,(sx2/2-1)*dx2,q
- DO x2=-sx2/2,sx2/2-1,j2
+ DO k2=-sx2/2,sx2/2-1,j2
+ x2=REAL(k2,8)
WRITE (15,'(ES12.2)') x2*dx2
END DO
WRITE (15,'(" </DataArray>")')
diff -r 9df56652aa63 -r 8ee30887f819 fourier.f90
--- a/fourier.f90 Wed Oct 26 13:37:50 2011 -0700
+++ b/fourier.f90 Mon Oct 31 15:55:10 2011 -0700
@@ -276,7 +276,13 @@ CONTAINS
iret=DftiCreateDescriptor(desc,DFTI_SINGLE,DFTI_REAL,3,size)
iret=DftiSetValue(desc,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX)
- WRITE_MKL_DEBUG_INFO(iret)
+ IF(iret.NE.0) THEN
+ IF(.NOT.DftiErrorClass(iret,DFTI_NO_ERROR)) THEN
+ WRITE_DEBUG_INFO
+ WRITE (0,*) DftiErrorMessage(iret)
+ STOP 1
+ END IF
+ END IF
IF (FFT_FORWARD == direction) THEN
scale=dx1*dx2*dx3
@@ -294,7 +300,13 @@ CONTAINS
iret=DftiComputeBackward(desc,data)
END IF
iret=DftiFreeDescriptor(desc)
- WRITE_MKL_DEBUG_INFO(iret)
+ IF(iret.NE.0) THEN
+ IF(.NOT.DftiErrorClass(iret,DFTI_NO_ERROR)) THEN
+ WRITE_DEBUG_INFO
+ WRITE (0,*) DftiErrorMessage(iret)
+ STOP 1
+ END IF
+ END IF
END SUBROUTINE fft3
#else
@@ -431,7 +443,13 @@ CONTAINS
iret=DftiCreateDescriptor(desc,DFTI_SINGLE,DFTI_REAL,2,size);
iret=DftiSetValue(desc,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX)
- WRITE_MKL_DEBUG_INFO(iret)
+ IF(iret.NE.0) THEN
+ IF(.NOT.DftiErrorClass(iret,DFTI_NO_ERROR)) THEN
+ WRITE_DEBUG_INFO
+ WRITE (0,*) DftiErrorMessage(iret)
+ STOP 1
+ END IF
+ END IF
IF (FFT_FORWARD == direction) THEN
scale=dx1*dx2
@@ -449,7 +467,13 @@ CONTAINS
iret=DftiComputeBackward(desc,data)
END IF
iret=DftiFreeDescriptor(desc)
- WRITE_MKL_DEBUG_INFO(iret)
+ IF(iret.NE.0) THEN
+ IF(.NOT.DftiErrorClass(iret,DFTI_NO_ERROR)) THEN
+ WRITE_DEBUG_INFO
+ WRITE (0,*) DftiErrorMessage(iret)
+ STOP 1
+ END IF
+ END IF
END SUBROUTINE fft2
#else
@@ -549,7 +573,13 @@ CONTAINS
REAL*4 :: scale
iret=DftiCreateDescriptor(desc,DFTI_SINGLE,DFTI_COMPLEX,1,sx)
- WRITE_MKL_DEBUG_INFO(iret)
+ IF(iret.NE.0) THEN
+ IF(.NOT.DftiErrorClass(iret,DFTI_NO_ERROR)) THEN
+ WRITE_DEBUG_INFO
+ WRITE (0,*) DftiErrorMessage(iret)
+ STOP 1
+ END IF
+ END IF
IF (FFT_FORWARD == direction) THEN
scale=dx
@@ -563,7 +593,13 @@ CONTAINS
iret=DftiComputeBackward(desc,data)
END IF
iret=DftiFreeDescriptor(desc)
- WRITE_MKL_DEBUG_INFO(iret)
+ IF(iret.NE.0) THEN
+ IF(.NOT.DftiErrorClass(iret,DFTI_NO_ERROR)) THEN
+ WRITE_DEBUG_INFO
+ WRITE (0,*) DftiErrorMessage(iret)
+ STOP 1
+ END IF
+ END IF
END SUBROUTINE fft1
#else
diff -r 9df56652aa63 -r 8ee30887f819 include.f90
--- a/include.f90 Wed Oct 26 13:37:50 2011 -0700
+++ b/include.f90 Mon Oct 31 15:55:10 2011 -0700
@@ -42,7 +42,7 @@
#ifdef IMKL_FFT
-#define WRITE_MKL_DEBUG_INFO(i) IF (i .NE. 0) THEN; IF (.NOT. DftiErrorClass(i,DFTI_NO_ERROR)) THEN; WRITE_DEBUG_INFO; WRITE (0,*) DftiErrorMessage(i); STOP 1; END IF; END IF
+#define WRITE_MKL_DEBUG_INFO(i) IF(i.NE.0)THEN;IF(.NOT.DftiErrorClass(i,DFTI_NO_ERROR))THEN;WRITE_DEBUG_INFO;WRITE (0,*) DftiErrorMessage(i);STOP 1;END IF;END IF
#endif
! adjust data alignment for the Fourier transform
diff -r 9df56652aa63 -r 8ee30887f819 relax.f90
--- a/relax.f90 Wed Oct 26 13:37:50 2011 -0700
+++ b/relax.f90 Mon Oct 31 15:55:10 2011 -0700
@@ -373,10 +373,10 @@ PROGRAM relax
IF (in%isoutputvtk) THEN
!filename=trim(in%wdir)//"/disp-000.vtr"
!CALL exportvtk_vectors(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,filename)
- filename=trim(in%wdir)//"/disp-000.vtk"
- title="coseismic displacement vector field"
- name="displacement"
- CALL exportvtk_vectors_legacy(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3, &
+ filename=trim(in%wdir)//"/disp-000.vtk"//char(0)
+ title="coseismic displacement vector field"//char(0)
+ name="displacement"//char(0)
+ CALL exportvtk_vectors_legacy(u1,u2,u3,in%sx1,in%sx2,in%sx3/8,in%dx1,in%dx2,in%dx3, &
4,4,8,filename,title,name)
!CALL exportvtk_vectors_slice(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,filename)
END IF
@@ -390,9 +390,9 @@ PROGRAM relax
#ifdef VTK
WRITE (digit4,'(I4.4)') oi-1
IF (in%isoutputvtk .AND. in%isoutputstress) THEN
- filename=trim(in%wdir)//"/sigma-"//digit4//".vtk"
- title="stress tensor field"
- name="stress"
+ filename=trim(in%wdir)//"/sigma-"//digit4//".vtk"//char(0)
+ title="stress tensor field"//char(0)
+ name="stress"//char(0)
CALL exportvtk_tensors_legacy(sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
4,4,8,filename,title,name)
END IF
@@ -442,7 +442,7 @@ PROGRAM relax
! and fault creep)
! 1- linear viscosity
IF (ALLOCATED(in%linearstruc)) THEN
- CALL viscouseigenstress(in%mu,in%linearstruc,in%linearweakzone, &
+ 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))
mech(1)=1
@@ -450,7 +450,7 @@ PROGRAM relax
! 2- powerlaw viscosity
IF (ALLOCATED(in%nonlinearstruc)) THEN
- CALL viscouseigenstress(in%mu,in%nonlinearstruc,in%nonlinearweakzone, &
+ 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
@@ -472,9 +472,9 @@ PROGRAM relax
#ifdef VTK
IF (in%isoutputvtk .AND. in%isoutputstress) THEN
WRITE (digit,'(I3.3)') oi-1
- filename=trim(in%wdir)//"/power-"//digit//".vtk"
- title="stress rate tensor field"
- name="power"
+ filename=trim(in%wdir)//"/power-"//digit//".vtk"//char(0)
+ title="stress rate tensor field"//char(0)
+ name="power"//char(0)
CALL exportvtk_tensors_legacy(moment,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
4,4,8,filename,title,name)
END IF
@@ -526,7 +526,7 @@ PROGRAM relax
IF (ALLOCATED(in%linearstruc)) THEN
! linear viscosity
v1=0
- CALL viscouseigenstress(in%mu,in%linearstruc,in%linearweakzone,sig, &
+ 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)
@@ -537,7 +537,7 @@ PROGRAM relax
IF (ALLOCATED(in%nonlinearstruc)) THEN
! powerlaw viscosity
v1=0
- CALL viscouseigenstress(in%mu,in%nonlinearstruc,in%nonlinearweakzone,sig, &
+ 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)
@@ -589,9 +589,9 @@ PROGRAM relax
WRITE (digit,'(I3.3)') oi
!filename=trim(in%wdir)//"/eqbf-"//digit//".vtr"
!CALL exportvtk_vectors(v1,v2,v3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,filename)
- filename=trim(in%wdir)//"/eqbf-"//digit//".vtk"
- title="instantaneous equivalent body-force rate vector field"
- name="body-force-rate"
+ filename=trim(in%wdir)//"/eqbf-"//digit//".vtk"//char(0)
+ title="instantaneous equivalent body-force rate vector field"//char(0)
+ name="body-force-rate"//char(0)
CALL exportvtk_vectors_legacy(v1,v2,v3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3, &
4,4,8,filename,title,name)
END IF
@@ -707,9 +707,9 @@ PROGRAM relax
! export total displacement in VTK XML format
!filename=trim(in%wdir)//"/disp-"//digit//".vtr"
!CALL exportvtk_vectors(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,filename)
- filename=trim(in%wdir)//"/disp-"//digit//".vtk"
- title="cumulative displacement vector field"
- name="displacement"
+ filename=trim(in%wdir)//"/disp-"//digit//".vtk"//char(0)
+ title="cumulative displacement vector field"//char(0)
+ name="displacement"//char(0)
CALL exportvtk_vectors_legacy(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3, &
4,4,8,filename,title,name)
!CALL exportvtk_vectors_slice(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,filename)
@@ -717,9 +717,9 @@ PROGRAM relax
! export instantaneous velocity in VTK XML format
!filename=trim(in%wdir)//"/vel-"//digit//".vtr"
!CALL exportvtk_vectors(v1,v2,v3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,filename)
- filename=trim(in%wdir)//"/vel-"//digit//".vtk"
- title="instantaneous velocity vector field"
- name="velocity"
+ filename=trim(in%wdir)//"/vel-"//digit//".vtk"//char(0)
+ title="instantaneous velocity vector field"//char(0)
+ name="velocity"//char(0)
CALL exportvtk_vectors_legacy(v1,v2,v3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
8,8,16,filename,title,name)
!CALL exportvtk_vectors_slice(v1,v2,v3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,filename)
@@ -742,9 +742,9 @@ PROGRAM relax
#ifdef VTK
WRITE (digit4,'(I4.4)') oi
IF (in%isoutputvtk .AND. in%isoutputstress) THEN
- filename=trim(in%wdir)//"/sigma-"//digit4//".vtk"
- title="stress tensor field"
- name="stress"
+ filename=trim(in%wdir)//"/sigma-"//digit4//".vtk"//char(0)
+ title="stress tensor field"//char(0)
+ name="stress"//char(0)
CALL exportvtk_tensors_legacy(sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
4,4,8,filename,title,name)
END IF
diff -r 9df56652aa63 -r 8ee30887f819 viscoelastic3d.f90
--- a/viscoelastic3d.f90 Wed Oct 26 13:37:50 2011 -0700
+++ b/viscoelastic3d.f90 Mon Oct 31 15:55:10 2011 -0700
@@ -120,11 +120,12 @@ CONTAINS
!!
!! \author sylvain barbot (08/30/08) - original form
!-----------------------------------------------------------------
- SUBROUTINE viscouseigenstress(mu,structure,ductilezones,sig,sx1,sx2,sx3, &
+ SUBROUTINE viscouseigenstress(mu,structure,ductilezones,nz,sig,sx1,sx2,sx3, &
dx1,dx2,dx3,moment,beta,maxwelltime,gamma)
REAL*8, INTENT(IN) :: mu,dx1,dx2,dx3,beta
+ INTEGER, INTENT(IN) :: nz
TYPE(LAYER_STRUCT), DIMENSION(:), INTENT(IN) :: structure
- TYPE(WEAK_STRUCT), DIMENSION(:), INTENT(IN) :: ductilezones
+ TYPE(WEAK_STRUCT), DIMENSION(nz), INTENT(IN) :: ductilezones
INTEGER, INTENT(IN) :: sx1,sx2,sx3
TYPE(TENSOR), INTENT(IN), DIMENSION(sx1,sx2,sx3) :: sig
TYPE(TENSOR), INTENT(OUT), DIMENSION(sx1,sx2,sx3) :: moment
@@ -173,12 +174,12 @@ CONTAINS
gammadot0=structure(i3)%gammadot0
! perturbation from isolated viscous zones
- dg0=dgammadot0(ductilezones,x1,x2,x3,beta)
+ dg0=dgammadot0(ductilezones,nz,x1,x2,x3,beta)
! local fluidity structure
gammadot0=gammadot0+dg0
- IF (1e-9 .GT. gammadot0) CYCLE
+ IF (1.0d-9 .GT. gammadot0) CYCLE
! local deviatoric stress
s=tensordeviatoric(sig(i1,i2,i3))
@@ -190,7 +191,7 @@ CONTAINS
tauc=tau-cohesion
! cohesion test
- IF (tauc .LE. 1e-9) CYCLE
+ IF (tauc .LE. 1.0d-9) CYCLE
! powerlaw viscosity
gammadot=gammadot0*(tauc/mu)**power
@@ -221,17 +222,15 @@ CONTAINS
!!
!! \author sylvain barbot (3/29/10) - original form
!---------------------------------------------------------
- REAL*8 FUNCTION dgammadot0(zones,x1,x2,x3,beta)
- TYPE(WEAK_STRUCT), INTENT(IN), DIMENSION(:) :: zones
+ REAL*8 FUNCTION dgammadot0(zones,n,x1,x2,x3,beta)
+ INTEGER, INTENT(IN) :: n
+ TYPE(WEAK_STRUCT), INTENT(IN), DIMENSION(nz) :: zones
REAL*8, INTENT(IN) :: x1,x2,x3,beta
REAL*8 :: dg,x,y,z,L,W,D,strike,dip,LM
REAL*8 :: cstrike,sstrike,cdip,sdip, &
xr,yr,zr,x2r,Wp,Lp,Dp,x1s,x2s,x3s
- INTEGER :: n,i
-
- ! number of ductile zones
- n=SIZE(zones,1)
+ INTEGER :: i
! default is no change in fluidity
dgammadot0=0._8
@@ -239,9 +238,15 @@ CONTAINS
DO i=1,n
! retrieve weak zone geometry
dg=zones(i)%dgammadot0
- x=zones(i)%x;y=zones(i)%y;z=zones(i)%z
- W=zones(i)%length;L=zones(i)%width;D=zones(i)%thickness
- strike=zones(i)%strike;dip=zones(i)%dip
+
+ x=zones(i)%x
+ y=zones(i)%y
+ z=zones(i)%z
+ W=zones(i)%length
+ L=zones(i)%width
+ D=zones(i)%thickness
+ strike=zones(i)%strike
+ dip=zones(i)%dip
! effective tapered dimensions
Wp=W*(1._8+2._8*beta)/2._8
diff -r 9df56652aa63 -r 8ee30887f819 waf
Binary file waf has changed
diff -r 9df56652aa63 -r 8ee30887f819 wscript
--- a/wscript Wed Oct 26 13:37:50 2011 -0700
+++ b/wscript Mon Oct 31 15:55:10 2011 -0700
@@ -48,13 +48,13 @@ def options(opt):
other.add_option('--use-ctfft', action='store_true', default=False,
help='Use slow internal CTFFT instead of MKL or FFTW')
-def ifort_modifier_darwin(conf):
- from waflib.Tools import fc_config
- fc_config.fortran_modifier_darwin(conf)
+# def ifort_modifier_darwin(conf):
+# from waflib.Tools import fc_config
+# fc_config.fortran_modifier_darwin(conf)
def configure(cnf):
- from waflib.Configure import conf
- conf(ifort_modifier_darwin)
+ # from waflib.Configure import conf
+ # conf(ifort_modifier_darwin)
cnf.load('compiler_c compiler_fc')
# Find Proj
@@ -199,6 +199,7 @@ def build(bld):
'elastic3d.f90',
'friction3d.f90',
'viscoelastic3d.f90',
+ 'writevtk.c',
'writegrd4.2.c',
'proj.c',
'export.f90',
More information about the CIG-COMMITS
mailing list