[cig-commits] r19206 - seismo/2D/SPECFEM2D/trunk/UTILS
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Nov 16 03:57:38 PST 2011
Author: dkomati1
Date: 2011-11-16 03:57:37 -0800 (Wed, 16 Nov 2011)
New Revision: 19206
Added:
seismo/2D/SPECFEM2D/trunk/UTILS/create_sinus_Paul_Cristini_ficsinmod.py
seismo/2D/SPECFEM2D/trunk/UTILS/process_DATA_Par_files_to_update_their_format_when_new_parameters_are_added.py
seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul_Cristini.py
Removed:
seismo/2D/SPECFEM2D/trunk/UTILS/create_sinus_Paul_ficsinmod.py
seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul.python
Log:
added Paul Cristini's script to process all DATA/Par_files to update their format automatically when new parameters are added
Copied: seismo/2D/SPECFEM2D/trunk/UTILS/create_sinus_Paul_Cristini_ficsinmod.py (from rev 19205, seismo/2D/SPECFEM2D/trunk/UTILS/create_sinus_Paul_ficsinmod.py)
===================================================================
--- seismo/2D/SPECFEM2D/trunk/UTILS/create_sinus_Paul_Cristini_ficsinmod.py (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/UTILS/create_sinus_Paul_Cristini_ficsinmod.py 2011-11-16 11:57:37 UTC (rev 19206)
@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+from pylab import *
+from os import *
+from sys import *
+import math
+
+xmin,xmax,dx=0.,4000.,5.
+ymin,ymax=0.,4000.
+
+y0=2000.
+# Amplitude de la sinusoide
+a=50.
+# Nombre de periodes
+n=50
+#
+x = arange(xmin,xmax+dx,dx)
+y = y0 + a*sin(2.*pi*n*(x-xmin)/(xmax-xmin))
+
+x1,y1,x3,y3=xmin,y0,xmax,y0
+x2,y2=(xmin+xmax)/2.,y0+a*3./2.
+
+ma=(y2-y1)/(x2-x1)
+mb=(y3-y2)/(x3-x2)
+
+x1j=complex(x1,y1)
+mx1j=abs(x1j)
+
+xc=(ma*mb*(y1-y3)+mb*(x1+x2)-ma*(x2+x3))/2./(mb-ma)
+yc=-(xc-(x1+x2)/2.)/ma+(y1+y2)/2.
+xcj=complex(xc,yc)
+
+R=abs(x1j-xcj)
+
+ang1=math.atan2(y1-yc,x1-xc)
+ang2=math.atan2(y3-yc,x3-xc)
+angb=msort([ang1,ang2])
+# Nombre de points pour la sinusoide
+nsin=5000.
+da = fabs((ang1-ang2)/nsin)
+ang = arange(angb[0],angb[1]+da,da)
+
+xi = xc + R*cos(ang)
+yi = yc + R*sin(ang)
+
+
+# Interface circulaire modulee
+
+xim = xc + cos(ang)*(R + a*sin(2.*pi*n*(ang-ang[0])/(angb[1]-angb[0])))
+yim = yc + sin(ang)*(R + a*sin(2.*pi*n*(ang-ang[0])/(angb[1]-angb[0])))
+
+#
+# Preparation des donnees en vue ecriture
+#
+DataS=zeros((len(x),2),Float)
+DataC=zeros((len(xi),2),Float)
+DataM=zeros((len(xim),2),Float)
+DataS[:,0]=x
+DataS[:,1]=y
+DataC[:,0]=xi[::-1]
+DataC[:,1]=yi[::-1]
+DataM[:,0]=xim[::-1]
+DataM[:,1]=yim[::-1]
+
+DataM=asarray(DataM)
+DataS=asarray(DataS)
+DataC=asarray(DataC)
+
+plot(x,y,'k',[xc],[yc],'ro')
+plot(x,y,'k',xi,yi,'r')
+plot(xim,yim,'b')
+axis([xmin,xmax,ymin,ymax])
+grid(True)
+title('Interface sinusoidale')
+#
+# Ecriture du fichier interface
+#
+NbElemSpec=150
+#
+# Format d'ecriture en ascii des coordonnées
+fmt='%f'
+#
+SEM_DIR='/home/cristini/SEM/SEM_2D_Dimitri'
+#
+chdir(SEM_DIR+'/DATA')
+#
+FileMod='InterfaceCercleMod'+str(int(a))+'_'+str(n)+'.dat'
+print FileMod
+#
+f = file(FileMod,'w')
+#
+f.write("#\n# Number of interfaces\n#\n")
+f.write("3\n")
+f.write("#\n#\n#\n")
+f.write("#\n# Interface 1 (Bottom of the mesh)\n#\n")
+f.write('2\n%i %i\n%i %i\n'%(xmin,ymin,xmax,ymin))
+f.write("#\n# Interface 2\n#\n")
+f.write('%i\n'%len(xim))
+for row in DataM:
+ f.write(' '.join([fmt%val for val in row]) + '\n')
+f.write("#\n# Interface 3\n#\n")
+f.write('2\n%i %i\n%i %i\n'%(xmin,ymax,xmax,ymax))
+f.write("#\n#Number of spectral elements\n#\n")
+f.write("#\n#\n#\n")
+f.write("%i\n"%NbElemSpec)
+f.write("#\n#\n#\n")
+f.write("%i\n"%NbElemSpec)
+f.close()
+#
+show()
+
Deleted: seismo/2D/SPECFEM2D/trunk/UTILS/create_sinus_Paul_ficsinmod.py
===================================================================
--- seismo/2D/SPECFEM2D/trunk/UTILS/create_sinus_Paul_ficsinmod.py 2011-11-16 11:32:57 UTC (rev 19205)
+++ seismo/2D/SPECFEM2D/trunk/UTILS/create_sinus_Paul_ficsinmod.py 2011-11-16 11:57:37 UTC (rev 19206)
@@ -1,110 +0,0 @@
-#!/usr/bin/env python
-from pylab import *
-from os import *
-from sys import *
-import math
-
-xmin,xmax,dx=0.,4000.,5.
-ymin,ymax=0.,4000.
-
-y0=2000.
-# Amplitude de la sinusoide
-a=50.
-# Nombre de periodes
-n=50
-#
-x = arange(xmin,xmax+dx,dx)
-y = y0 + a*sin(2.*pi*n*(x-xmin)/(xmax-xmin))
-
-x1,y1,x3,y3=xmin,y0,xmax,y0
-x2,y2=(xmin+xmax)/2.,y0+a*3./2.
-
-ma=(y2-y1)/(x2-x1)
-mb=(y3-y2)/(x3-x2)
-
-x1j=complex(x1,y1)
-mx1j=abs(x1j)
-
-xc=(ma*mb*(y1-y3)+mb*(x1+x2)-ma*(x2+x3))/2./(mb-ma)
-yc=-(xc-(x1+x2)/2.)/ma+(y1+y2)/2.
-xcj=complex(xc,yc)
-
-R=abs(x1j-xcj)
-
-ang1=math.atan2(y1-yc,x1-xc)
-ang2=math.atan2(y3-yc,x3-xc)
-angb=msort([ang1,ang2])
-# Nombre de points pour la sinusoide
-nsin=5000.
-da = fabs((ang1-ang2)/nsin)
-ang = arange(angb[0],angb[1]+da,da)
-
-xi = xc + R*cos(ang)
-yi = yc + R*sin(ang)
-
-
-# Interface circulaire modulee
-
-xim = xc + cos(ang)*(R + a*sin(2.*pi*n*(ang-ang[0])/(angb[1]-angb[0])))
-yim = yc + sin(ang)*(R + a*sin(2.*pi*n*(ang-ang[0])/(angb[1]-angb[0])))
-
-#
-# Preparation des donnees en vue ecriture
-#
-DataS=zeros((len(x),2),Float)
-DataC=zeros((len(xi),2),Float)
-DataM=zeros((len(xim),2),Float)
-DataS[:,0]=x
-DataS[:,1]=y
-DataC[:,0]=xi[::-1]
-DataC[:,1]=yi[::-1]
-DataM[:,0]=xim[::-1]
-DataM[:,1]=yim[::-1]
-
-DataM=asarray(DataM)
-DataS=asarray(DataS)
-DataC=asarray(DataC)
-
-plot(x,y,'k',[xc],[yc],'ro')
-plot(x,y,'k',xi,yi,'r')
-plot(xim,yim,'b')
-axis([xmin,xmax,ymin,ymax])
-grid(True)
-title('Interface sinusoidale')
-#
-# Ecriture du fichier interface
-#
-NbElemSpec=150
-#
-# Format d'ecriture en ascii des coordonnées
-fmt='%f'
-#
-SEM_DIR='/home/cristini/SEM/SEM_2D_Dimitri'
-#
-chdir(SEM_DIR+'/DATA')
-#
-FileMod='InterfaceCercleMod'+str(int(a))+'_'+str(n)+'.dat'
-print FileMod
-#
-f = file(FileMod,'w')
-#
-f.write("#\n# Number of interfaces\n#\n")
-f.write("3\n")
-f.write("#\n#\n#\n")
-f.write("#\n# Interface 1 (Bottom of the mesh)\n#\n")
-f.write('2\n%i %i\n%i %i\n'%(xmin,ymin,xmax,ymin))
-f.write("#\n# Interface 2\n#\n")
-f.write('%i\n'%len(xim))
-for row in DataM:
- f.write(' '.join([fmt%val for val in row]) + '\n')
-f.write("#\n# Interface 3\n#\n")
-f.write('2\n%i %i\n%i %i\n'%(xmin,ymax,xmax,ymax))
-f.write("#\n#Number of spectral elements\n#\n")
-f.write("#\n#\n#\n")
-f.write("%i\n"%NbElemSpec)
-f.write("#\n#\n#\n")
-f.write("%i\n"%NbElemSpec)
-f.close()
-#
-show()
-
Added: seismo/2D/SPECFEM2D/trunk/UTILS/process_DATA_Par_files_to_update_their_format_when_new_parameters_are_added.py
===================================================================
--- seismo/2D/SPECFEM2D/trunk/UTILS/process_DATA_Par_files_to_update_their_format_when_new_parameters_are_added.py (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/UTILS/process_DATA_Par_files_to_update_their_format_when_new_parameters_are_added.py 2011-11-16 11:57:37 UTC (rev 19206)
@@ -0,0 +1,75 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Nov 15 09:00:00 2011
+
+Process Par_file to update them to the new format of release 19201
+
+Usage : "python PathTo/SPECFEM2D/UTILS/ProcessParFileTor19201.py PathTo/filename"
+
+ at author: Cristini Paul, Laboratoire de Mecanique et d'Acoustique, CNRS, Marseille, France
+"""
+import sys
+from shutil import move
+from os.path import exists
+
+def ProcessParfileTor19201(fic):
+ # New additions to the Par_file
+ a1='PERFORM_CUTHILL_MCKEE = .true. # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering\n'
+ a2='USER_T0 = 0.0d0 # use this t0 as earliest starting time rather than the automatically calculated one\n'
+ a3='SU_FORMAT = .false. # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)\n'
+ a4='factor_subsample_image = 1 # factor to subsample color images output by the code (useful for very large models)\n'+ \
+ 'POWER_DISPLAY_COLOR = 0.30d0 # non linear display to enhance small amplitudes in color images\n'+ \
+ 'DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true. # display acoustic layers as constant blue in JPEG images, because they likely correspond to water\n'
+ a5='US_LETTER = .false. # US letter paper or European A4\n'+ \
+ 'USE_SNAPSHOT_NUMBER_IN_FILENAME = .false. # use snapshot number in the file name of JPEG color snapshots instead of the time step\n'
+ a6='\n# for horizontal periodic conditions: detect common points between left and right edges\n'+ \
+ 'ADD_PERIODIC_CONDITIONS = .false.\n\n'+ \
+ '# horizontal periodicity distance for periodic conditions\n'+ \
+ 'PERIODIC_horiz_dist = 0.3597d0\n\n'+ \
+ '# grid point detection tolerance for periodic conditions\n'+ \
+ "PERIODIC_DETECT_TOL = 3.3334d-6\n"
+ #
+ f = open(fic,'r')
+ ligs= f.readlines()
+ f.close()
+ #
+ # Ajout des parametres supplementaires
+ # On verifie si le fichier n'a pas deja ete traite
+ if not (ligs[0].endswith('r19201\n')):
+ ligs[0]=ligs[0][:-1]+' r19201\n' # On indique que le fichier est traite pour cette release
+ #
+ Ct=0
+ for ilg, lig in enumerate(ligs):
+ if lig.startswith('partitioning'):
+ ligs.insert(ilg+1,a1)
+ Ct+=1
+ if lig.startswith('deltat'):
+ ligs.insert(ilg+1,a2)
+ Ct+=1
+ if lig.startswith('rec_normal'):
+ ligs.insert(ilg+1,a3)
+ Ct+=1
+ if lig.startswith('subsamp'):
+ ligs[ilg]=string.replace(ligs[ilg],'subsamp ','subsamp_postscript',1)
+ ligs.insert(ilg+1,a4)
+ Ct+=1
+ if lig.startswith('sizemax'):
+ ligs.insert(ilg+1,a5)
+ Ct+=1
+ if lig.startswith('absorbing_conditions'):
+ ligs.insert(ilg+1,a6)
+ Ct+=1
+ print fic, Ct
+# #
+ move(fic,fic+'.old')
+ fm = open(fic,'w')
+ fm.writelines(ligs)
+ fm.close()
+ print 'File : '+fic+' processed'
+ else:
+ print 'File : '+fic+' already processed'
+ return
+
+if __name__=='__main__':
+ ProcessParfileTor19201(sys.argv[1])
+
\ No newline at end of file
Property changes on: seismo/2D/SPECFEM2D/trunk/UTILS/process_DATA_Par_files_to_update_their_format_when_new_parameters_are_added.py
___________________________________________________________________
Name: svn:executable
+ *
Deleted: seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul.python
===================================================================
--- seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul.python 2011-11-16 11:32:57 UTC (rev 19205)
+++ seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul.python 2011-11-16 11:57:37 UTC (rev 19206)
@@ -1,60 +0,0 @@
-#!/usr/bin/env python
-#
-# Python code to generate Seismic Unix files for Ux and Uz and plot them
-#
-from os import *
-import os.path, string
-#
-def createSU():
- SEM=getcwd()
- filename=SEM+'/DATA/Par_file'
- #
- # Variables to be put in SU header
- #
- if path.exists(filename):
- variables=['nt','dt','sismostype']
- else:
- print 'No Par_file found !'
- return
- #
- # Open the file and get the lines
- #
- f = file(filename,'r')
- lignes= f.readlines()
- f.close()
- # Get the title
- for ligne in lignes:
- lsplit=string.split(ligne)
- if lsplit!= []:
- if lsplit[0]=='title':
- title=' '.join(lsplit[2:])
- break
- print '#'*50
- print '# SU file creation for '+title
- print '#'*50
- # Get the variables
- for var in variables:
- for ligne in lignes:
- lsplit=string.split(ligne)
- if lsplit!= []:
- if lsplit[0]==var:
- exec var+'='+string.replace(string.split(''.join(ligne))[2],'d','e')
- break
- #
- print sismostype
- chdir(SEM+'/OUTPUT_FILES')
- labels='label1="Time" label2="Receivers"'
- # Create headers and Su file
- if sismostype==4:
- ordres=['suaddhead < pressure_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > pressure_file.su']
- ordres.append('suxwigb < pressure_file.su perc=96 '+labels+' title=" Pressure : '+title+'"&')
- else:
- ordres=['suaddhead < Ux_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > Ux_file.su']
- ordres.append('suaddhead < Uz_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > Uz_file.su')
- ordres.append('suxwigb < Ux_file.su perc=96 '+labels+' title=" Ux : '+title+'"&')
- ordres.append('suxwigb < Uz_file.su xbox=600 perc=96 '+labels+' title=" Uz : '+title+'"&')
- #
- for i in range(len(ordres)):
- system(ordres[i])
-if __name__=='__main__':
- createSU()
Copied: seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul_Cristini.py (from rev 19205, seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul.python)
===================================================================
--- seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul_Cristini.py (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul_Cristini.py 2011-11-16 11:57:37 UTC (rev 19206)
@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+#
+# Python code to generate Seismic Unix files for Ux and Uz and plot them
+#
+from os import *
+import os.path, string
+#
+def createSU():
+ SEM=getcwd()
+ filename=SEM+'/DATA/Par_file'
+ #
+ # Variables to be put in SU header
+ #
+ if path.exists(filename):
+ variables=['nt','dt','sismostype']
+ else:
+ print 'No Par_file found !'
+ return
+ #
+ # Open the file and get the lines
+ #
+ f = file(filename,'r')
+ lignes= f.readlines()
+ f.close()
+ # Get the title
+ for ligne in lignes:
+ lsplit=string.split(ligne)
+ if lsplit!= []:
+ if lsplit[0]=='title':
+ title=' '.join(lsplit[2:])
+ break
+ print '#'*50
+ print '# SU file creation for '+title
+ print '#'*50
+ # Get the variables
+ for var in variables:
+ for ligne in lignes:
+ lsplit=string.split(ligne)
+ if lsplit!= []:
+ if lsplit[0]==var:
+ exec var+'='+string.replace(string.split(''.join(ligne))[2],'d','e')
+ break
+ #
+ print sismostype
+ chdir(SEM+'/OUTPUT_FILES')
+ labels='label1="Time" label2="Receivers"'
+ # Create headers and Su file
+ if sismostype==4:
+ ordres=['suaddhead < pressure_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > pressure_file.su']
+ ordres.append('suxwigb < pressure_file.su perc=96 '+labels+' title=" Pressure : '+title+'"&')
+ else:
+ ordres=['suaddhead < Ux_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > Ux_file.su']
+ ordres.append('suaddhead < Uz_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > Uz_file.su')
+ ordres.append('suxwigb < Ux_file.su perc=96 '+labels+' title=" Ux : '+title+'"&')
+ ordres.append('suxwigb < Uz_file.su xbox=600 perc=96 '+labels+' title=" Uz : '+title+'"&')
+ #
+ for i in range(len(ordres)):
+ system(ordres[i])
+if __name__=='__main__':
+ createSU()
More information about the CIG-COMMITS
mailing list