hh_toy_start.c 1.89 KB
Newer Older
Simone Butzer's avatar
Simone Butzer committed
1
/*
2
 *   homogeneous half space
Simone Butzer's avatar
Simone Butzer committed
3 4 5 6 7
 *   last update 02.11.02, T. Bohlen
 */

#include "fd.h"

Tilman Steinweg's avatar
Tilman Steinweg committed
8
void model(st_model *mod) {
Simone Butzer's avatar
Simone Butzer committed
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

	/*-----------------------------------------------------------------------*/
	/* extern variables */

	extern float DT, *FL, TAU;
	extern int NX, NY, NZ, NXG, NYG, NZG, POS[4], L;
	extern FILE *FP;
	/* local variables */
	float muv, piv, ws;
	float Vp, Vs, Rho;
	float *pts=NULL, ts=0.0, tp=0.0, sumu, sumpi;
	int i, j, k, l, ii, jj, kk;

	/* 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;*/
25
	/*const float vp2=7000.0, vs2=3900.0, rho2=2800.0;*/
Simone Butzer's avatar
Simone Butzer committed
26 27

	/*-----------------------------------------------------------------------*/
28
	sumu=0.0;
Simone Butzer's avatar
Simone Butzer committed
29 30 31
	sumpi=0.0;

	fprintf(FP," start creation of toy-model");
32 33

	if (L) {
Simone Butzer's avatar
Simone Butzer committed
34 35
		/* vector for maxwellbodies */
		pts=vector(1,L);
36 37

		for (l=1; l<=L; l++) {
Simone Butzer's avatar
Simone Butzer committed
38
			pts[l]=1.0/(2.0*PI*FL[l]);
Tilman Steinweg's avatar
Tilman Steinweg committed
39
			mod->eta[l]=DT/pts[l];
Simone Butzer's avatar
Simone Butzer committed
40
		}
41 42

		ts=TAU;
Simone Butzer's avatar
Simone Butzer committed
43 44
		tp=TAU;
		ws=2.0*PI*FL[1];
45 46

		for (l=1; l<=L; l++) {
Simone Butzer's avatar
Simone Butzer committed
47 48 49
			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]));
		}
50
	}
Simone Butzer's avatar
Simone Butzer committed
51 52

	/* loop over global grid */
53 54 55
	for (j=1; j<=NYG; j++) { /*vertical*/
		for (i=1; i<=NXG; i++) {
			for (k=1; k<=NZG; k++) {
Simone Butzer's avatar
Simone Butzer committed
56 57


58 59 60
				Vp=vp1;
				Vs=vs1;
				Rho=rho1;
Simone Butzer's avatar
Simone Butzer committed
61 62 63 64

				muv=Vs*Vs*Rho/(1.0+sumu);
				piv=Vp*Vp*Rho/(1.0+sumpi);

65
				/* only the PE which belongs to the current global gridpoint
Simone Butzer's avatar
Simone Butzer committed
66
				  is saving model parameters in his local arrays */
67 68 69
				if ((POS[1]==((i-1)/NX)) &&
				        (POS[2]==((j-1)/NY)) &&
				        (POS[3]==((k-1)/NZ))) {
Simone Butzer's avatar
Simone Butzer committed
70 71 72 73
					ii=i-POS[1]*NX;
					jj=j-POS[2]*NY;
					kk=k-POS[3]*NZ;

74
					if (L) {
Tilman Steinweg's avatar
Tilman Steinweg committed
75 76
						mod->taus[jj][ii][kk]=ts;
						mod->taup[jj][ii][kk]=tp;
Simone Butzer's avatar
Simone Butzer committed
77 78
					}

Tilman Steinweg's avatar
Tilman Steinweg committed
79 80 81
					mod->u[jj][ii][kk]=muv;
					mod->rho[jj][ii][kk]=Rho;
					mod->pi[jj][ii][kk]=piv;
Simone Butzer's avatar
Simone Butzer committed
82 83 84
				}
			}
		}
85
	}
Simone Butzer's avatar
Simone Butzer committed
86

87
	
Simone Butzer's avatar
Simone Butzer committed
88 89 90


	MPI_Barrier(MPI_COMM_WORLD);
91 92 93 94

	if (L)	{
		free_vector(pts,1,L);
	}
Simone Butzer's avatar
Simone Butzer committed
95 96 97 98
}