Commit 9760baf2 authored by Tilman Steinweg's avatar Tilman Steinweg

rotated coordinate system of the input/output

(z<->y): y is now vertical for input/output like in the whole code.
Note that the output and input of the models stayed the same.
The order of the loops over the model should be the same in every function. (The most efficent order)
The models have as fast (n1) directon the z (horizontal) component.
n1/n2/n3 = z/x/y. So with ximage you will see horizontal (n1(z),n2(x)) slices forom top to bottom.
parent 7d13a5aa
......@@ -9,8 +9,8 @@ nx=nx/outx;ny=ny/outy;nz=nz/outz;
fignum=71;
file_inp1='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/grad/toy.gradvp_160.00Hz_it1001'; %preconditioned gradient
file_inp2='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/grad/toy.gradvp_160.00Hz_it1'; %"raw" gradient
file_inp1='../par/grad/toy_grad.vp_160.00Hz_it1001'; %preconditioned gradient
file_inp2='../par/grad/toy_grad.vp_160.00Hz_it1'; %"raw" gradient
......
......@@ -6,13 +6,13 @@ nx=160; ny=184; nz=160; %ny:vertical
outx=1; outy=1; outz=1;
dh=0.8;
nx=nx/outx;ny=ny/outy;nz=nz/outz;
fignum=764;
fignum=13;
%file_inp1='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/model/toy_real.vp';
%file_inp2='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/model/toy_real.vp';
file_inp1='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/model/toy.vp_it60';
file_inp2='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/model/toy.vp_it60';
file_inp1='../par/model/toy.vp_it0';
file_inp2='../par/model/toy.vp_it60';
phi1=0; %rotation angles to x-z plane of first and second plane
......
......@@ -2,7 +2,7 @@
# PARAMETER FILE FOR IFOS3D
#-----------------------------------------------------------------
#
# Note that z denotes the vertical direction !
# Note that y denotes the vertical direction !
#
#-----------------------------------------------------------------------------------
#------------------------ MODELING PARAMETERS --------------------------------------
......@@ -15,8 +15,8 @@ number_of_processors_in_z-direction_(NPROCZ) = 2
#
#-------------------- 3-D Grid -----------------------------------
number_of_gridpoints_in_x-direction_(NX) = 160
number_of_gridpoints_in_y-direction_(NY) = 160
number_of_gridpoints_in_z-direction_(NZ) = 184
number_of_gridpoints_in_y-direction_(NY) = 184
number_of_gridpoints_in_z-direction_(NZ) = 160
distance_between_gridpoints(in_m)_in_x-direction_(DX) = 0.8
distance_between_gridpoints(in_m)_in_y-direction_(DY) = 0.8
distance_between_gridpoints(in_m)_in_z-direction_(DZ) = 0.8
......@@ -35,7 +35,7 @@ timestep_(in_seconds)_(DT) = 5.0e-5
#--------------------Source---------------------------------------
# Shape_of_source-signal:
(ricker=1;fumue=2;from_SIGNAL_FILE=3;SIN**3=4;deltapulse=5)_(SOURCE_SHAPE) = 4
point_source_(explosive=1;force_in_x=2;in_y=3;in_z=4;custom=5)_(SOURCE_TYPE) = 4
point_source_(explosive=1;force_in_x=2;in_z=3;in_y=4;custom=5)_(SOURCE_TYPE) = 4
# If SOURCE_TYPE <5 the following two lines are ignored
force_angle_between_x_y_(in_degree)_(ALPHA) = 45.0
force_angle_between_x_z_(in_degree)_(BETA) = 45.0
......@@ -81,7 +81,7 @@ output_of_snapshots_(SNAP)(yes>0) = 0
first_snapshot_(in_sec)_(TSNAP1) = 0.01
last_snapshot_(in_sec)_(TSNAP2) = 0.24
increment_(in_sec)_(TSNAPINC) = 0.0075
# Note that z denotes the vertical direction !
# Note that y denotes the vertical direction !
increment_x-direction_(IDX) = 1
increment_y-direction_(IDY) = 1
increment_z-direction_(IDZ) = 1
......@@ -91,9 +91,9 @@ basic_filename_(SNAP_FILE) = ./snap/back
#if SNAP = 1,2 the following line is ignored
(SNAP_PLANE) = 1
#output of snapshots as energy wihout sign SNAP_PLANE=1
#energy with sign true for x-z-plane SNAP_PLANE=2
#energy with sign true for x-y-plane SNAP_PLANE=3
#energy with sign true for y-z-plane SNAP_PLANE=4
#energy with sign true for x-y-plane SNAP_PLANE=2
#energy with sign true for x-z-plane SNAP_PLANE=3
#energy with sign true for z-y-plane SNAP_PLANE=4
#
#----------------------Receiver-----------------------------------
output_of_seismograms_(SEISMO) = 1
......@@ -110,7 +110,7 @@ read_receiver_positions_(READREC) = 2
REC_FILE = ./receiver/receiver.dat
reference_point_for_receiver_coordinate_system_(REFREC) = 0.0 , 0.0 , 0.0
# if READREC=1 the following receiver options are ignored
# Note that z denotes the vertical direction !
# Note that y denotes the vertical direction !
#
# --------------------Receiver line -------------------------------
position_of_first_receiver_(in_m)_(XREC1,YREC1,ZREC1) = 90.0 , 90.0 , 90.0
......@@ -124,7 +124,7 @@ number_of_planes_(no<=0)_(REC_ARRAY) = 1
depth_of_first_(upper)_plane_(in_m)_(REC_ARRAY_DEPTH) = 24.0
vertical_distance_between_planes_(in_m)_(REC_ARRAY_DIST) = 30.0
distance_between_receivers_in_x-direction_(in_gridpoints)_(DRX) = 10
distance_between_receivers_in_y-direction_(in_gridpoints)_(DRY) = 10
distance_between_receivers_in_y-direction_(in_gridpoints)_(DRZ) = 10
#
#-------------------- Seismogram Output-----------------------------
samplingrate_and_timelag_(in_timesteps!)_(NDT,NDTSHIFT) = 1 , 0
......
{
"MODELING PARAMETERS" : "comment",
"Note that y denotes the vertical direction !" : "comment",
"Domain Decomposition" : "comment",
"NPROCX" : "2",
"NPROCY" : "2",
......@@ -8,8 +10,8 @@
"3-D Grid" : "comment",
"NX" : "160",
"NY" : "160",
"NZ" : "184",
"NY" : "184",
"NZ" : "160",
"DX" : "0.8",
"DY" : "0.8",
"DZ" : "0.8",
......@@ -81,7 +83,7 @@
"REC_ARRAY_DEPTH" : "24.0",
"REC_ARRAY_DIST" : "30.0",
"DRX" : "10",
"DRY" : "10",
"DRZ" : "10",
"Seismograms" : "comment",
......
{
"MODELING PARAMETERS" : "comment",
"Note that z denotes the vertical direction !" : "comment",
"Note that y denotes the vertical direction !" : "comment",
"Domain Decomposition" : "comment",
"NPROCX" : "2",
......@@ -10,8 +10,8 @@
"3-D Grid" : "comment",
"NX" : "160",
"NY" : "160",
"NZ" : "184",
"NY" : "184",
"NZ" : "160",
"DX" : "0.8",
"DY" : "0.8",
"DZ" : "0.8",
......@@ -36,11 +36,11 @@
"External signal input instead of SOURCE_SHAPE" : "comment",
"SIGNAL_FILE" : "./STF/stf.su",
"SOURCE_TYPE" : "4",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_y=3;force_in_z=4;custom=5" : "comment",
"SOURCE_TYPE values (point_source): explosive=1;force_in_x=2;force_in_z=3;force_in_y=4;custom=5" : "comment",
"If SOURCE_TYPE <5 the following two lines are ignored" : "comment",
"ALPHA" : "45.0",
"BETA" : "45.0",
"ALPHA=force_angle_between_x_y_(in_degree) (BETA between x_z)" : "comment",
"ALPHA=force_angle_between_x_y(vertical)_(in_degree) (BETA between x_z)" : "comment",
"read_source_positions_from_SOURCE_FILE_(yes=1;Plane wave=2) " : "comment",
"SRCREC" : "1",
"SOURCE_FILE" : "./sources/sources_toy.dat",
......@@ -112,9 +112,9 @@
"if SNAP !=3 the following line is ignored" : "comment",
"SNAP_PLANE" : "1",
"output of snapshots as energy wihout sign SNAP_PLANE=1" : "comment",
"energy with sign true for x-z-plane SNAP_PLANE=2" : "comment",
"energy with sign true for x-y-plane SNAP_PLANE=3" : "comment",
"energy with sign true for y-z-plane SNAP_PLANE=4" : "comment",
"energy with sign true for x-y-plane SNAP_PLANE=2" : "comment",
"energy with sign true for x-z-plane SNAP_PLANE=3" : "comment",
"energy with sign true for z-x-plane SNAP_PLANE=4" : "comment",
"Receiver" : "comment",
"output_of_seismograms" : "comment",
......@@ -131,7 +131,7 @@
"reference_point_for_receiver_coordinate_system" : "comment",
"REFRECX, REFRECY, REFRECZ" : "0.0 , 0.0, 0.0",
"if READREC=1 the following three lines are ignored" : "comment",
"Note that z denotes the vertical direction !" : "comment",
"Note that y denotes the vertical direction !" : "comment",
"position_of_first_receiver_(in_m)" : "comment",
"XREC1, YREC1, ZREC1" : "90.0 , 90.0, 90.0",
"position_of_last_receiver_(in_m)" : "comment",
......@@ -149,7 +149,7 @@
"REC_ARRAY_DIST" : "30.0",
"distance_between_receivers_in_x/y-direction" : "comment",
"DRX" : "10",
"DRY" : "10",
"DRZ" : "10",
"Seismograms" : "comment",
......
......@@ -2,7 +2,7 @@
# PARAMETER FILE FOR IFOS3D
#-----------------------------------------------------------------
#
# Note that z denotes the vertical direction !
# Note that y denotes the vertical direction !
#
#-----------------------------------------------------------------------------------
#------------------------ MODELING PARAMETERS --------------------------------------
......@@ -15,8 +15,8 @@ number_of_processors_in_z-direction_(NPROCZ) = 2
#
#-------------------- 3-D Grid -----------------------------------
number_of_gridpoints_in_x-direction_(NX) = 160
number_of_gridpoints_in_y-direction_(NY) = 160
number_of_gridpoints_in_z-direction_(NZ) = 184
number_of_gridpoints_in_y-direction_(NY) = 184
number_of_gridpoints_in_z-direction_(NZ) = 160
distance_between_gridpoints(in_m)_in_x-direction_(DX) = 0.8
distance_between_gridpoints(in_m)_in_y-direction_(DY) = 0.8
distance_between_gridpoints(in_m)_in_z-direction_(DZ) = 0.8
......@@ -35,7 +35,7 @@ timestep_(in_seconds)_(DT) = 5.0e-5
#--------------------Source---------------------------------------
# Shape_of_source-signal:
(ricker=1;fumue=2;from_SIGNAL_FILE=3;SIN**3=4;deltapulse=5)_(SOURCE_SHAPE) = 4
point_source_(explosive=1;force_in_x=2;in_y=3;in_z=4;custom=5)_(SOURCE_TYPE) = 4
point_source_(explosive=1;force_in_x=2;in_z=3;in_y=4;custom=5)_(SOURCE_TYPE) = 4
# If SOURCE_TYPE <5 the following two lines are ignored
force_angle_between_x_y_(in_degree)_(ALPHA) = 45.0
force_angle_between_x_z_(in_degree)_(BETA) = 45.0
......@@ -81,7 +81,7 @@ output_of_snapshots_(SNAP)(yes>0) = 0
first_snapshot_(in_sec)_(TSNAP1) = 0.01
last_snapshot_(in_sec)_(TSNAP2) = 0.24
increment_(in_sec)_(TSNAPINC) = 0.0075
# Note that z denotes the vertical direction !
# Note that y denotes the vertical direction !
increment_x-direction_(IDX) = 1
increment_y-direction_(IDY) = 1
increment_z-direction_(IDZ) = 1
......@@ -91,9 +91,9 @@ basic_filename_(SNAP_FILE) = ./snap/back
#if SNAP = 1,2 the following line is ignored
(SNAP_PLANE) = 1
#output of snapshots as energy wihout sign SNAP_PLANE=1
#energy with sign true for x-z-plane SNAP_PLANE=2
#energy with sign true for x-y-plane SNAP_PLANE=3
#energy with sign true for y-z-plane SNAP_PLANE=4
#energy with sign true for x-y-plane SNAP_PLANE=2
#energy with sign true for x-z-plane SNAP_PLANE=3
#energy with sign true for z-x-plane SNAP_PLANE=4
#
#----------------------Receiver-----------------------------------
output_of_seismograms_(SEISMO) = 1
......@@ -110,7 +110,7 @@ read_receiver_positions_(READREC) = 2
REC_FILE = ./receiver/receiver.dat
reference_point_for_receiver_coordinate_system_(REFREC) = 0.0 , 0.0 , 0.0
# if READREC=1 the following receiver options are ignored
# Note that z denotes the vertical direction !
# Note that y denotes the vertical direction !
#
# --------------------Receiver line -------------------------------
position_of_first_receiver_(in_m)_(XREC1,YREC1,ZREC1) = 90.0 , 90.0 , 90.0
......@@ -124,7 +124,7 @@ number_of_planes_(no<=0)_(REC_ARRAY) = 1
depth_of_first_(upper)_plane_(in_m)_(REC_ARRAY_DEPTH) = 24.0
vertical_distance_between_planes_(in_m)_(REC_ARRAY_DIST) = 30.0
distance_between_receivers_in_x-direction_(in_gridpoints)_(DRX) = 10
distance_between_receivers_in_y-direction_(in_gridpoints)_(DRY) = 10
distance_between_receivers_in_y-direction_(in_gridpoints)_(DRZ) = 10
#
#-------------------- Seismogram Output-----------------------------
samplingrate_and_timelag_(in_timesteps!)_(NDT,NDTSHIFT) = 1 , 0
......
......@@ -8,8 +8,8 @@
"3-D Grid" : "comment",
"NX" : "160",
"NY" : "160",
"NZ" : "184",
"NY" : "184",
"NZ" : "160",
"DX" : "0.8",
"DY" : "0.8",
"DZ" : "0.8",
......@@ -58,7 +58,7 @@
"REC_ARRAY_DEPTH" : "24.0",
"REC_ARRAY_DIST" : "30.0",
"DRX" : "10",
"DRY" : "10",
"DRZ" : "10",
"Seismograms" : "comment",
"NDT" : "1",
......
{
"MODELING PARAMETERS" : "comment",
"Note that y denotes the vertical direction !" : "comment",
"Domain Decomposition" : "comment",
"NPROCX" : "2",
"NPROCY" : "2",
......@@ -8,8 +10,8 @@
"3-D Grid" : "comment",
"NX" : "160",
"NY" : "160",
"NZ" : "184",
"NY" : "184",
"NZ" : "160",
"DX" : "0.8",
"DY" : "0.8",
"DZ" : "0.8",
......@@ -52,14 +54,14 @@
"Receiver" : "comment",
"SEISMO" : "1",
"READREC" : "2",
"READREC" : "2",
"Receiver array" : "comment",
"REC_ARRAY" : "1",
"REC_ARRAY_DEPTH" : "24.0",
"REC_ARRAY_DIST" : "30.0",
"DRX" : "10",
"DRY" : "10",
"DRZ" : "10",
"Seismograms" : "comment",
......
......@@ -16,5 +16,5 @@
#
# Symbols:
# XSRC= x-coordinate of source point [meter]
# YSRC= y-coordinate of source point [meter]
# ZSRC= z-coordinate of source point [meter] (vertical)
# YSRC= y-coordinate of source point [meter] (vertical)
# ZSRC= z-coordinate of source point [meter]
......@@ -26,8 +26,8 @@
%
% Symbols:
% XSRC= x-coordinate of source point [meter]
% YSRC= y-coordinate of source point [meter]
% ZSRC= z-coordinate of source point [meter] (vertical)
% YSRC= y-coordinate of source point [meter] (vertical)
% 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
......
32.0 32.0 92.0 0.005 300.0 1.0
32.0 52.8 92.0 0.005 300.0 1.0
32.0 74.4 92.0 0.005 300.0 1.0
32.0 96.0 92.0 0.005 300.0 1.0
64.0 32.0 92.0 0.005 300.0 1.0
64.0 52.8 92.0 0.005 300.0 1.0
64.0 74.4 92.0 0.005 300.0 1.0
64.0 96.0 92.0 0.005 300.0 1.0
96.0 32.0 92.0 0.005 300.0 1.0
96.0 52.8 92.0 0.005 300.0 1.0
96.0 74.4 92.0 0.005 300.0 1.0
96.0 96.0 92.0 0.005 300.0 1.0
32.0 92.0 32.0 0.005 300.0 1.0
32.0 92.0 52.8 0.005 300.0 1.0
32.0 92.0 74.4 0.005 300.0 1.0
32.0 92.0 96.0 0.005 300.0 1.0
64.0 92.0 32.0 0.005 300.0 1.0
64.0 92.0 52.8 0.005 300.0 1.0
64.0 92.0 74.4 0.005 300.0 1.0
64.0 92.0 96.0 0.005 300.0 1.0
96.0 92.0 32.0 0.005 300.0 1.0
96.0 92.0 52.8 0.005 300.0 1.0
96.0 92.0 74.4 0.005 300.0 1.0
96.0 92.0 96.0 0.005 300.0 1.0
......@@ -26,8 +26,8 @@
%
% Symbols:
% XSRC= x-coordinate of source point [meter]
% YSRC= y-coordinate of source point [meter]
% ZSRC= z-coordinate of source point [meter] (vertical)
% YSRC= y-coordinate of source point [meter] (vertical)
% 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
......
This diff is collapsed.
......@@ -26,7 +26,7 @@
void exchange_par(void){
/* declaration of extern variables */
extern int NX, NY, NZ, SOURCE_SHAPE, SOURCE_TYPE, SNAP, SNAP_FORMAT, SNAP_PLANE;
extern int DRX, DRY, L, SRCREC, FDORDER;
extern int DRX, DRZ, L, SRCREC, FDORDER;
extern float DX, DY, DZ, TIME, DT, *FL, TS, TAU, PLANE_WAVE_DEPTH, PHI;
extern float XREC1, XREC2, YREC1, YREC2, ZREC1, ZREC2;
extern float ALPHA, BETA, VPPML;
......@@ -108,7 +108,7 @@ void exchange_par(void){
idum[12] = FREE_SURF;
idum[13] = SNAP;
idum[14] = DRX;
idum[15] = DRY;
idum[15] = DRZ;
idum[16] = BOUNDARY;
idum[17] = REC_ARRAY;
idum[18] = SRCREC;
......@@ -236,7 +236,7 @@ void exchange_par(void){
FREE_SURF = idum[12];
SNAP = idum[13];
DRX = idum[14];
DRY = idum[15];
DRZ = idum[15];
BOUNDARY = idum[16];
REC_ARRAY = idum[17];
SRCREC = idum[18];
......
......@@ -31,7 +31,7 @@ float REFREC[4]={0.0, 0.0, 0.0, 0.0}, DAMPING=8.0, VPPML, FPML;
int SEISMO, NDT, NDTSHIFT, NGEOPH, SEIS_FORMAT[6]={0, 0, 0, 0, 0, 0}, FREE_SURF, READMOD, MOD_FORMAT, READREC, REC_ARRAY, LOG, FDORDER, FW=0, ABS_TYPE, BLOCK;
int NX, NY, NZ=1, NT, SOURCE_SHAPE, SOURCE_TYPE, SNAP, SNAP_FORMAT, BOUNDARY, SRCREC, SNAP_PLANE;
float ALPHA, BETA;
int NXG, NYG, NZG, IDX, IDY, IDZ, L=1, NX1, NX2, NY1, NY2, NZ1, NZ2, DRX, DRY, RUN_MULTIPLE_SHOTS, FDCOEFF;
int NXG, NYG, NZG, IDX, IDY, IDZ, L=1, NX1, NX2, NY1, NY2, NZ1, NZ2, DRX, DRZ, RUN_MULTIPLE_SHOTS, FDCOEFF;
char SNAP_FILE[STRING_SIZE], SOURCE_FILE[STRING_SIZE], SIGNAL_FILE[STRING_SIZE], INV_FILE[STRING_SIZE];
char MFILE[STRING_SIZE], REC_FILE[STRING_SIZE];
char SEIS_FILE[STRING_SIZE],MOD_FILE[STRING_SIZE];
......
......@@ -67,20 +67,12 @@ float *** taus, float *** taup, float * eta){
}
/* loop over global grid */
for (k=1;k<=NZG;k++){
for (j=1;j<=NYG;j++){
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
for (k=1;k<=NZG;k++){
Vp=vp1; Vs=vs1; Rho=rho1;
/*for usability reasons, "z" - as commonly used - denotes the depth (vertical direction),
however, internally "y" is used for the vertical coordinate,
we simply switch the "y" and "z" coordinate as read in the input file,
therefore use the variable "j" to, e.g., calculate/address the vertical coordinate:
x=(float)i*DX;
y=(float)k*DZ;
z=(float)j*DY;*/
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
......@@ -105,22 +97,14 @@ float *** taus, float *** taup, float * eta){
}
}
for (k=50;k<=80;k++){
for (j=50;j<=80;j++){
for (i=55;i<=85;i++){
for (j=50;j<=80;j++){
for (k=50;k<=80;k++){
/*for (j=(NYG/2+1);j<=NYG;j++){*/
Vp=0.0; Vs=0.0; Rho=0.0;
Vp=vp2; Vs=vs2; Rho=rho2;
/*for usability reasons, "z" - as commonly used - denotes the depth (vertical direction),
however, internally "y" is used for the vertical coordinate,
we simply switch the "y" and "z" coordinate as read in the input file,
therefore use the variable "j" to, e.g., calculate/address the vertical coordinate:
x=(float)i*DX;
y=(float)k*DZ;
z=(float)j*DY;*/
muv=Vs*Vs*Rho/(1.0+sumu);
piv=Vp*Vp*Rho/(1.0+sumpi);
......
/*
* homogeneous half space
* 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){
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float *eta) {
/*-----------------------------------------------------------------------*/
/* extern variables */
......@@ -19,199 +19,80 @@ float *** taus, float *** taup, float * eta){
float Vp, Vs, Rho;
float *pts=NULL, ts=0.0, tp=0.0, sumu, sumpi;
int i, j, k, l, ii, jj, kk;
int kasten=0;
/* parameters for layer 1 */
const float vp1=6200.0, vs1=3600.0, rho1=2800.0;
/*const float vp2=6200.0, vs2=3600.0, rho2=2800.0;*/
/*const float vp2=7000.0, vs2=3900.0, rho2=2800.0;*/
/*const float vp2=7000.0, vs2=3900.0, rho2=2800.0;*/
/*-----------------------------------------------------------------------*/
sumu=0.0;
sumu=0.0;
sumpi=0.0;
fprintf(FP," start creation of toy-model");
if(L){
if (L) {
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
for (l=1; l<=L; l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
ts=TAU;
ts=TAU;
tp=TAU;
ws=2.0*PI*FL[1];
for (l=1;l<=L;l++){
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=vp1; Vs=vs1; Rho=rho1;
/*for usability reasons, "z" - as commonly used - denotes the depth (vertical direction),
however, internally "y" is used for the vertical coordinate,
we simply switch the "y" and "z" coordinate as read in the input file,
therefore use the variable "j" to, e.g., calculate/address the vertical coordinate:
x=(float)i*DX;
y=(float)k*DZ;
z=(float)j*DY;*/
for (j=1; j<=NYG; j++) { /*vertical*/
for (i=1; i<=NXG; i++) {
for (k=1; k<=NZG; k++) {
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;
if(L){
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;
}
}
}
}
if(kasten){
for (k=50;k<=90;k++){
for (i=50;i<=80;i++){
for (j=50;j<=65;j++){ /*vertical*/
Vp=0.0; Vs=0.0; Rho=0.0;
Vp=vp1+500.0; Vs=vs1+300; Rho=rho1;
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
/* 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))){
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;
if(L){
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;
}
}
}
}
for (k=50;k<=90;k++){
for (i=81;i<=110;i++){
for (j=50;j<=95;j++){
Vp=0.0; Vs=0.0; Rho=0.0;
Vp=vp1-250; Vs=vs1-100; 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;
if(L){
if (L) {
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;
}
}
}
} for (k=50;k<=90;k++){
for (i=50;i<=80;i++){
for (j=66;j<=95;j++){
Vp=0.0; Vs=0.0; Rho=0.0;
Vp=vp1+300; Vs=vs1+200; 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;
if(L){
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;
}
}
}
} for (k=91;k<=110;k++){
for (i=50;i<=110;i++){
for (j=50;j<=95;j++){
}
Vp=0.0; Vs=0.0; Rho=0.0;
Vp=vp1-400; Vs=vs1-150; 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;
if(L){
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;
}
}
}
}
}
MPI_Barrier(MPI_COMM_WORLD);
if(L) free_vector(pts,1,L);
if (L) {
free_vector(pts,1,L);
}
}
......
/*
* homogeneous half space
* 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){
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float *eta) {
/*-----------------------------------------------------------------------*/
/* extern variables */
......@@ -24,61 +24,60 @@ float *** taus, float *** taup, float * eta){
/* parameters for layer 1 */
const float vp1=6200.0, vs1=3600.0, rho1=2800.0;
/*const float vp2=6200.0, vs2=3600.0, rho2=2800.0;*/
/*const float vp2=7000.0, vs2=3900.0, rho2=2800.0;*/
/*const float vp2=7000.0, vs2=3900.0, rho2=2800.0;*/