...
 
Commits (15)
To compile the program do:
cd sofi3D/src
then
"make sofi3D" to compile the standard staggered grid version (SSG)
(you probably need to change the compiler options in sofi3D/src/Makefile)
To run the program on 8 CPUs do
cd sofi3D/par
mpirun -np 8 ../bin/sofi3D ./in_and_out/sofi3D.json > ./in_and_out/sofi3D.jout
you may also use the shell script startSOFI3D.sh
The file sofi3D/par/in_and_out/sofi3D.jout shows the obtained screen output.
The modelling parameters are specified in sofi3D.json. I hope that the
parameters are more or less self-explanatory. The synthetic seismograms
and snapshots are written to the files specified in :
sofi3D/par/in_and_out/sofi3D.json.
In the current distribution the model is generated on the fly by the
function sofi3D/src/hh_elastic.c . This function generates a homogeneous
medium with Vp=3500 m/s, Vs=2000 m/s, and rho=2000 kg/(m*m*m). The function
readmod.c can be used to read model info from external grid files.
See readmod.c for the format of the external files. You can change the
function which generates the model grid at the beginning of sofi3D/src/Makefile.
For more information, have a look at the SOFI3D manual which you can compile in ./sofi3D/doc/guide_sofi3D.
Additional you can find the manual at https://git.scc.kit.edu/GPIAG-Software/SOFI3D/wikis/home
See README.md in par/
# What is SOFI3D?
This is a minimalistic version of SOFI3D for the _WAVE_-Project.
**Do not use this version for production!**
SOFI3D stands for Seismic mOdeling with FInite differences and denotes our 3-D viscoelastic time domain massive parallel modeling code.
The manual and a reference paper is included in the download archive or can be downloaded [here](https://git.scc.kit.edu/GPIAG-Software/SOFI3D/wikis/home)
The manual is included in the download archive or can be downloaded [here](https://git.scc.kit.edu/GPIAG-Software/SOFI3D/wikis/home).
# Download and Newsletter
## Content
You can Download the [latest stable Release](https://git.scc.kit.edu/GPIAG-Software/SOFI3D/tree/Release) or the current [beta-version](https://git.scc.kit.edu/GPIAG-Software/SOFI3D/tree/master).
Also a SOFI3D branch with additional benchmarks is available: [overnightbuilt](https://git.scc.kit.edu/GPIAG-Software/SOFI3D/tree/overnightbuilt)
- Source code: `src/`
- Run scripts etc.: `par/`
- Executable: `bin/`
- Documentation: `doc/` (for full SOFI3D version)
- _WAVE_-Documentation: `par/WAVE_doc/`
To receive news and updates please [register](https://www.gpi.kit.edu/Software-WS.php) on the email list sofi@lists.kit.edu.
Please use this list also to ask questions on using the software or to report problems or bugs.
\ No newline at end of file
See file `par/README.md` to get started.
/*
* homogeneous half space
* last update 02.11.02, T. Bohlen
*/
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern float DT, *FL, TAU, TS, DX, DY;
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L, MYID;
extern char MFILE[STRING_SIZE];
/* local variables */
float muv, piv, ws;
float Vp, Vs, Rho, y,x , ampund, aund, shiftx, undf;
float *pts, ts, tp, sumu, sumpi;
int i, j, k, l, ii, jj, kk;
char modfile[STRING_SIZE];
/*----------------------SPECIFY HERE MATERIAL PROPERTIES -------------- */
ampund=400; aund=0.01; shiftx=1200;
/*-----------------------------------------------------------------------*/
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
ws=2.0*PI*FL[1];
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
/* loop over global grid */
for (k=1;k<=NZG;k++){
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
Vp=3500.0; Vs=2000.0; Rho=2000.0;
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
}
}
}
}
writemod(MFILE,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(MFILE,3);
free_vector(pts,1,L);
}
1
2610.0 2610.0 2100.0 0.0 5.0 1.0
%
% Definition of (distributed) source positions:
% position, time-delay, centre frequency, and amplitude
%
% (comment line is indicated by # as first character)
%
% Symbols:
% XSRC= x-coordinate of source point [meter]
% YSRC= y-coordinate of source point [meter] (y denotes vertical direction in FDMPI !)
% ZSRC= z-coordinate of source point [meter]
% TD= excitation time (time-delay) for source node [s]
% FC= centre frequency of source signal [Hz]
% AMP= maximum amplitude of source signal
% NSRC= number of source positions
%
% Parameters for each source node (one source node per line):
%
% NSRC (one line)
% XSRC ZSRC YSRC TD FC AMP (NSRC lines)
/*
* homogeneous space with a tunnel
* last update 05.02.99, T. Bohlen
*/
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern float DT, DH, *FL, TAU, TS;
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L, MYID;
extern char MFILE[STRING_SIZE];
/* local variables */
float muv, piv, ws;
float Vp, Vs, Rho;
float x, y, z;
float *pts, ts, tp, sumu, sumpi;
int i, j, k, l, ii, jj, kk;
char modfile[STRING_SIZE];
/* background */
const float vp0=6000.0, vs0=3000.0, rho0=2000.0;
/* position of chrack */
const float vpt=0.0, vst=0.0, rhot=1e-3;
const float xc=5.0, yc=5.0, zc=5.0, xd=15.0, yd=10.0, zd=10.0;
/*-----------------------------------------------------------------------*/
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
ws=2.0*PI*FL[1];
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
/* loop over global grid */
for (k=1;k<=NZG;k++){
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
Vp=vp0; Vs=vs0; Rho=rho0;
x=(float)i*DH;
y=(float)j*DH;
z=(float)k*DH;
if ((x>xc)&&(x<(xc+xd))&&(y>yc)&&(y<(yc+yd))&&(z>zc)&&(z<(zc+zd))){
Vp=vpt; Vs=vst; Rho=rhot;}
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
}
}
}
}
sprintf(modfile,"model/model.bin");
writemod(modfile,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
free_vector(pts,1,L);
}
/*
* Dike-Model-3D
* D.Koehn
* Kiel, 8.11.2003
*
* update: 16.11.2003 - added "groundwater-line-polynom"
* 20.12.2003 - fixed polynomial coefficient bug
* 18.03.2004 - NR don't run on Altix - calculate Polynomial Coefficients with MATLAB !
* 28.03.2004 - included dike-shape-function f3
*/
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern float DT, DH, *FL, TAU, TS;
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L, MYID;
extern char MFILE[STRING_SIZE];
/* local variables */
float muv, piv, ws;
float Vp, Vs, Rho, xi, yi, zi;
float *pts, ts, tp, sumu, sumpi;
float f1, f2;
int i, j, k, l, ii, jj, kk;
char modfile[STRING_SIZE];
/* some mathematical constants*/
const float pim = 4.0*atan(1.0);
const float pic = pim/180.0;
/* geometrical parameters*/
const float alpha=36.32*pic, delta=24.52*pic;
const float l1=1.48, l2=1.0, l3=2.485, h=1.19;
const float dp=0.05, da=0.04, T=0.3;
/* parameters for air */
/*const float Vp0=0.0, Vs0=0.0001, Rho0=1.25;*/
const float Vp0=0.0, Vs0=0.0001, Rho0=1.25;
/* parameters for water*/
//const float Vp1=1500.0, Vs1=1.0e-4, Rho1=1000.0;
/* parameters for Lehmkern */
const float Vp2=500.0, Vs2=250.0, Rho2=1800.0;
/* parameters for Plexiglas (PMMA)*/
/*const float Vp2w=3150.0, Vs2w=1570.0, Rho2w=1180.0;*/
/*const float Vp2w=0.0, Vs2w=0.0001, Rho2w=1.25;*/
const float Vp2w=0.0, Vs2w=0.0001, Rho2w=1.25;
/*-----------------------------------------------------------------------*/
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
ws=2.0*PI*FL[1];
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
/* loop over global grid */
for (k=1;k<=NZG;k++){
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
xi=(float)i*DH;
yi=(float)j*DH;
zi=(float)k*DH;
/* fill model with clay*/
Vp=Vp2; Vs=Vs2; Rho=Rho2;
/* area lambda 1 */
/* ------------ */
/* calculate "dike-shape-function - f1" */
/* to describe the boundary of the dike on landside. */
f1=(-tan(alpha)*xi)+h+(da);
/* air */
if((yi<=f1)){
Vp=Vp0; Vs=Vs0; Rho=Rho0;
}
/* area lambda 2 */
/* ------------ */
/* calculate "dike-shape-functions-f2" */
/* to describe the boundary of the dike on seaside. */
f2 = (tan(delta) * xi) - (tan(delta)*((2*dp)+l1+l2+l3+da)) + h;
/* air */
if(yi<f2){
Vp=Vp0; Vs=Vs0; Rho=Rho0;
}
/* boundaries*/
/* fill front with air */
if(zi<(da+dp)){
Vp=Vp0; Vs=Vs0; Rho=Rho0;
}
/* front boundary */
if((zi<(da+dp))&&(zi>da)){
Vp=Vp2w; Vs=Vs2w; Rho=Rho2w;
}
/* background boundary */
/* fill back with plexiglas*/
if(zi>(da+dp+T)){
Vp=Vp2w; Vs=Vs2w; Rho=Rho2w;
}
/* lower boundary */
if((yi>(da+h-dp))){
Vp=Vp2w; Vs=Vs2w; Rho=Rho2w;
}
/* left boundary */
if((xi<=da+dp)){
Vp=Vp2w; Vs=Vs2w; Rho=Rho2w;
}
/* right boundary */
if((xi>=da+dp+l1+l2+l3)){
Vp=Vp2w; Vs=Vs2w; Rho=Rho2w;
}
/* upper boundary*/
if((yi<=da)){
Vp=Vp0; Vs=Vs0; Rho=Rho0;
}
/* lower boundary */
if((yi>=(h+da))){
Vp=Vp0; Vs=Vs0; Rho=Rho0;
}
/* left boundary */
if((xi<=da)){
Vp=Vp0; Vs=Vs0; Rho=Rho0;
}
/* right boundary */
if((xi>da+dp+l1+l2+l3+dp)){
Vp=Vp0; Vs=Vs0; Rho=Rho0;
}
/* fill back with air */
if(zi>(da+dp+T+dp)){
Vp=Vp0; Vs=Vs0; Rho=Rho0;
}
/* fill back with air */
if(k<4){
Vp=Vp0; Vs=Vs0; Rho=Rho0;
}
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
}
}
}
}
writemod(MFILE,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(MFILE,3);
free_vector(pts,1,L);
}
/*
* homogeneous space with a tunnel
* last update 05.02.99, T. Bohlen
*/
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern float DT, DH, *FL, TS, TAU, REFREC[4];
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L, MYID;
extern char MFILE[STRING_SIZE];
/* local variables */
float muv, piv, ws;
float Vp, Vs, Rho;
float x, y, z, r, yhill;
float *pts, ts, tp, sumu, sumpi;
int i, j, k, l, ii, jj, kk;
/* air */
const float vp1=0.0, vs1=0.0, rho1=1.25, Qp1=10000.0, Qs1=10000.0;
/* const float vp1=3000.0, vs1=1500.0, rho1=1200.0, Qp1=10000.0, Qs1=10000.0;*/
/* hill */
const float h=1000.0, a=1000.0, x0=REFREC[1], z0=REFREC[3], y0=REFREC[2],
vp2=3000.0, vs2=1500.0, rho2=1200.0;
/*-----------------------------------------------------------------------*/
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
ws=2.0*PI*FL[1];
/* loop over global grid */
for (k=1;k<=NZG;k++){
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
x=(float)i*DH-x0;
y=(float)j*DH;
z=(float)k*DH-z0;
r=sqrt(x*x+z*z);
yhill=h*(1-exp(-(r*r)/(a*a)))+y0;
if (y<yhill){
Vp=vp1; Vs=vs1; Rho=rho1; tp=2.0/Qp1; ts=2.0/Qs1;}
else {
Vp=vp2; Vs=vs2; Rho=rho2; tp=TAU; ts=TAU;}
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
}
}
}
}
/* writemod(MFILE,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(MFILE,3);
*/
free_vector(pts,1,L);
}
/*
* homogeneous space with a tunnel
* last update 05.02.99, T. Bohlen
*/
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern float DT, DH, *FL, TS;
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L, MYID;
extern char MFILE[STRING_SIZE];
/* local variables */
float muv, piv, ws;
float Vp, Vs, Rho;
float y, depth;
float *pts, ts, tp, sumu, sumpi;
int i, j, k, l, ii, jj, kk;
FILE *fp;
char modfile[STRING_SIZE];
/* mantle */
const float vp2=8000.0, vs2=4000.0, rho2=3190.0, Qp2=1000.0, Qs2=1000.0;
/* crust */
const float vp1=6200.0, vs1=3000.0, rho1=2900.0, Qp1=500.0, Qs1=500.0;
/*-----------------------------------------------------------------------*/
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
ws=2.0*PI*FL[1];
printf(" PE %d: opening file %s ...",MYID,MFILE);
fp=fopen(MFILE,"r");
if (fp==NULL) err(" could not open file MFILE ");
/* loop over global grid */
for (k=1;k<=NZG;k++){
for (i=1;i<=NXG;i++){
fscanf(fp,"%e", &depth);
depth*=1000.0;
for (j=1;j<=NYG;j++){
y=(float)j*DH;
if (y<=depth){
Vp=vp1; Vs=vs1; Rho=rho1; tp=2.0/Qp1; ts=2.0/Qs1;}
else{
Vp=vp2; Vs=vs2; Rho=rho2; tp=2.0/Qp2; ts=2.0/Qs2;}
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
}
}
}
}
fclose(fp);
/* sprintf(modfile,"model/iceland.bin");
writemod(modfile,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
*/
free_vector(pts,1,L);
}
/* $Id: jastram_acoustic.c,v 1.1 2007/10/28 18:02:15 koehn Exp $ */
/*
* jastram model
*/
#include "fd.h"
void model_acoustic(float *** rho, float *** pi){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], MYID;
extern char MFILE[STRING_SIZE];
extern float DX, DY;
/* local variables */
float piv;
float Vp, Rho, x, y, undf, ampund,shiftx, aund;
int i, j, k, ii, jj, kk;
char filename[STRING_SIZE];
/*-----------------------------------------------------------------------*/
aund = 0.01;
ampund = 400.0;
shiftx = 1200.0;
for (j=1;j<=NYG;j++){
for (i=1;i<=NXG;i++){
for (k=1;k<=NZG;k++){
x=i*DX;
y=j*DY;
Vp=1900.0; Rho=1000.0;
/* calculate undulation function */
undf=ampund*exp(-(aund*(x-shiftx))*(aund*(x-shiftx)));
if(y<=200.0){
Vp=1500; Rho=1000;
}
if(y>200.0){
Vp=1000; Rho=1200;
}
if(y>250.0){
Vp=2000; Rho=1500;
}
if(y>(460.0+undf)){
Vp=3000; Rho=2000;
}
piv=Vp*Vp*Rho;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
}
}
}
}
sprintf(filename,"%s.fdmpi.pi",MFILE);
writemod(filename,pi,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(filename,3);
}
/*
* homogeneous half space
* last update 02.11.02, T. Bohlen
*/
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern float DT, *FL, TAU, TS, DX, DY;
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L, MYID;
extern char MFILE[STRING_SIZE];
/* local variables */
float muv, piv, ws;
float Vp, Vs, Rho, y,x , ampund, aund, shiftx, undf;
float *pts, ts, tp, sumu, sumpi;
int i, j, k, l, ii, jj, kk;
char modfile[STRING_SIZE];
/*----------------------SPECIFY HERE MATERIAL PROPERTIES -------------- */
ampund=400; aund=0.01; shiftx=1200;
/*-----------------------------------------------------------------------*/
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
ws=2.0*PI*FL[1];
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
/* loop over global grid */
for (k=1;k<=NZG;k++){
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
/* calculate cartesian coordinates*/
x=i*DX;
y=j*DY;
/* calculate undulation function */
undf=ampund*exp(-(aund*(x-shiftx))*(aund*(x-shiftx)));
if(y<=200.0){
Vp=1500; Vs=1e-6; Rho=1000;
}
if(y>200.0){
Vp=1000; Vs=400; Rho=1200;
}
if(y>250.0){
Vp=2000; Vs=1155; Rho=1500;
}
if(y>(460.0+undf)){
Vp=3000; Vs=1732; Rho=2000;
}
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
}
}
}
}
writemod(MFILE,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(MFILE,3);
free_vector(pts,1,L);
}
1
1200.0 500.0 75.0 0.0 50.0 1.0
%
% Definition of (distributed) source positions:
% position, time-delay, centre frequency, and amplitude
%
% (comment line is indicated by # as first character)
%
% Symbols:
% XSRC= x-coordinate of source point [meter]
% YSRC= y-coordinate of source point [meter] (y denotes vertical direction in FDMPI !)
% ZSRC= z-coordinate of source point [meter]
% TD= excitation time (time-delay) for source node [s]
% FC= centre frequency of source signal [Hz]
% AMP= maximum amplitude of source signal
% NSRC= number of source positions
%
% Parameters for each source node (one source node per line):
%
% NSRC (one line)
% XSRC ZSRC YSRC TD FC AMP (NSRC lines)
/* $Id: benchmark_acoustic.c,v 1.1 2007/11/03 01:07:37 koehn Exp $ */
/*
* jastram model
*/
#include "fd.h"
void model_acoustic(float *** rho, float *** pi){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], MYID;
extern char MFILE[STRING_SIZE];
extern float DX, DY;
/* local variables */
float piv;
float Vp, Rho;
int i, j, k, ii, jj, kk;
char filename[STRING_SIZE];
/*-----------------------------------------------------------------------*/
for (j=1;j<=NYG;j++){
for (i=1;i<=NXG;i++){
for (k=1;k<=NZG;k++){
Vp=3500.0; Rho=2000.0;
piv=Vp*Vp*Rho;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
}
}
}
}
sprintf(filename,"%s.fdmpi.pi",MFILE);
writemod(filename,pi,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(filename,3);
}
/*
* homogeneous half space
* last update 02.11.02, T. Bohlen
*/
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern float DT, *FL, TAU, TS, DX, DY;
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L, MYID;
extern char MFILE[STRING_SIZE];
/* local variables */
float muv, piv, ws;
float Vp, Vs, Rho, y,x , ampund, aund, shiftx, undf;
float *pts, ts, tp, sumu, sumpi;
int i, j, k, l, ii, jj, kk;
char modfile[STRING_SIZE];
/*----------------------SPECIFY HERE MATERIAL PROPERTIES -------------- */
ampund=400; aund=0.01; shiftx=1200;
/*-----------------------------------------------------------------------*/
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
ws=2.0*PI*FL[1];
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
/* loop over global grid */
for (k=1;k<=NZG;k++){
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
Vp=3500.0; Vs=2000.0; Rho=2000.0;
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
}
}
}
}
writemod(MFILE,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(MFILE,3);
free_vector(pts,1,L);
}
4
2610.0 2610.0 1050.0 0.0 5.0 1.0
2610.0 2610.0 4200.0 0.0 5.0 1.0
2610.0 4610.0 2100.0 0.0 5.0 1.0
4610.0 2610.0 1050.0 0.0 5.0 1.0
%
% Definition of (distributed) source positions:
% position, time-delay, centre frequency, and amplitude
%
% (comment line is indicated by # as first character)
%
% Symbols:
% XSRC= x-coordinate of source point [meter]
% YSRC= y-coordinate of source point [meter] (y denotes vertical direction in FDMPI !)
% ZSRC= z-coordinate of source point [meter]
% TD= excitation time (time-delay) for source node [s]
% FC= centre frequency of source signal [Hz]
% AMP= maximum amplitude of source signal
% NSRC= number of source positions
%
% Parameters for each source node (one source node per line):
%
% NSRC (one line)
% XSRC ZSRC YSRC TD FC AMP (NSRC lines)
/*
* homogeneous half space
* last update 02.11.02, T. Bohlen
*/
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern float DT, DH, *FL, TAU, TS;
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L, MYID;
extern char MFILE[STRING_SIZE];
/* local variables */
float muv, piv, ws;
float Vp, Vs, Rho, y;
float *pts, ts, tp, sumu, sumpi;
int i, j, k, l, ii, jj, kk;
char modfile[STRING_SIZE];
/* parameters for layer 1 */
const float vp1=0.0, vs1=0.0001, rho1=1.25, iw=5.0;
/* parameters for layer 2 */
const float vp2=500.0, vs2=250.0, rho2=1800.0;
/*-----------------------------------------------------------------------*/
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
ws=2.0*PI*FL[1];
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
/* loop over global grid */
for (k=1;k<=NZG;k++){
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
Vp=vp2; Vs=vs2; Rho=rho2;
if ((i<=iw)||(j<=iw)||(k<=iw)||(k>=NZG-iw)||(i>=NXG-iw)||(j>=NYG-iw)){
Vp=vp1; Vs=vs1; Rho=rho1; }
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;