SPEED
Functions/Subroutines
GET_INITIAL_CONDITIONS.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine get_initial_conditions (nnod_loc, u0, u1, v1, ne_loc, cs_loc, cs_nnz_loc, nm, prop_mat, sdeg_mat, loc_n_num, xs_loc, ys_loc, zs_loc, mpi_id, dt)
 Set initial conditions.
 

Function/Subroutine Documentation

◆ get_initial_conditions()

subroutine get_initial_conditions ( integer*4  nnod_loc,
real*8, dimension(3*nnod_loc)  u0,
real*8, dimension(3*nnod_loc)  u1,
real*8, dimension(3*nnod_loc)  v1,
integer*4  ne_loc,
integer*4, dimension(0:cs_nnz_loc)  cs_loc,
integer*4  cs_nnz_loc,
integer*4  nm,
real*8, dimension(nm,4)  prop_mat,
integer*4, dimension(nm)  sdeg_mat,
integer*4, dimension(nnod_loc)  loc_n_num,
real*8, dimension(nnod_loc)  xs_loc,
real*8, dimension(nnod_loc)  ys_loc,
real*8, dimension(nnod_loc)  zs_loc,
integer*4  mpi_id,
real*8  dt 
)

Set initial conditions.

Author
Ilario Mazzieri
Date
September, 2013
Version
1.0
Parameters
[in]nnod_locnumber of local nodes
[in]ne_locnumber of local elements
[in]lengthof cs_loc
[in]cs_loclocal connectivity vector
[in]nmnumber of materials
[in]prop_matmaterial properties (rho,lambda,mu,gamma)
[in]sdeg_matpolinomial degree vector
[in]local_n_num
[in]xs_locx-coordinate of spectral nodes
[in]ys_locy-coordinate of spectral nodes
[in]zs_locz-coordinate of spectral nodes
[in]mpi_idMPI process id
[in]dttime step
[out]u0displacement at time -dt
[out]u1displacement at time 0
[out]v1velocity at time 0

Definition at line 44 of file GET_INITIAL_CONDITIONS.f90.

48
49
50 implicit none
51
52 character*100000 :: input_line
53 character*4 :: keyword
54
55 integer*4 :: nnod_loc, mpi_id, cs_nnz_loc
56 integer*4 :: nm, ne_loc, status
57 integer*4 :: ie, im, nn, fn,is,in,iaz,i,j,k,ileft,iright,nval
58 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
59
60 integer*4, dimension(nnod_loc) :: loc_n_num
61 integer*4, dimension(nm) :: sdeg_mat
62
63
64 real*8 :: time, term1,term2,term3, dt, pi, vp, omega, k1, k2, k3
65 real*8 :: a1,a2,a3,b1,b2,b3,c1,c2,c3,w1, rho,lambda,mu, cost1,cost2,cost3
66 real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j, u1_ex,u2_ex,u3_ex
67
68 real*8, dimension(nm,4) :: prop_mat
69 real*8, dimension(ne_loc) :: alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,alfa31,alfa32,alfa33
70 real*8, dimension(ne_loc) :: beta11,beta12,beta13,beta21,beta22,beta23,beta31,beta32,beta33
71 real*8, dimension(ne_loc) :: gamma1,gamma2,gamma3,delta1,delta2,delta3
72 real*8, dimension(3*nnod_loc) :: u0, u1, v1
73 real*8, dimension(nnod_loc) :: xs_loc, ys_loc, zs_loc
74 real*8, dimension(:), allocatable :: val_data
75
76
77
78 pi = 4.d0*datan(1.d0)
79
80! open(file='INITCOND.input',unit=40+mpi_id)
81!
82! do
83! read(40+mpi_id,'(A)',IOSTAT = status) input_line
84!
85! if (status.ne.0) exit
86!
87! keyword = input_line(1:4)
88!
89! ileft = 0
90! iright = len(input_line)
91! do i = 1,iright
92! if (input_line(i:i).eq.' ') exit
93! enddo
94! ileft = i
95!
96! if (keyword .eq. 'DATA') then
97! read(input_line(ileft:iright),*) nval
98! endif
99! enddo
100!
101! close(40+mpi_id)
102! allocate(val_data(nval))
103!
104! open(file='INITCOND.input',unit=40+mpi_id)
105!
106! do
107! read(40+mpi_id,'(A)',IOSTAT = status) input_line
108!
109! if (status.ne.0) exit
110!
111! keyword = input_line(1:4)
112!
113! ileft = 0
114! iright = len(input_line)
115! do i = 1,iright
116! if (input_line(i:i).eq.' ') exit
117! enddo
118! ileft = i
119!
120! if (keyword .eq. 'DATA') then
121! read(input_line(ileft:iright),*) nval, (val_data(j), j = 1, nval)
122! endif
123! enddo
124!
125! close(40+mpi_id)
126
127
128 do ie = 1,ne_loc
129
130 im = cs_loc(cs_loc(ie -1))
131 nn = sdeg_mat(im) +1
132 rho = prop_mat(im,1); lambda = prop_mat(im,2); mu = prop_mat(im,3)
133
134 k1 = 0;
135 k2 = 0;
136 k3 = 1;
137
138 omega = 120.d0*2*pi
139 vp = dsqrt((lambda+2*mu)/rho)
140 k3 = k3*omega/vp;
141
142! a1 = val_data(1); b1 = val_data(2); c1 = val_data(3)
143! a2 = val_data(4); b2 = val_data(5); c2 = val_data(6)
144! a3 = val_data(7); b3 = val_data(8); c3 = val_data(9)
145! w1 = val_data(10);
146
147 do k = 1,nn
148 do j = 1,nn
149 do i = 1,nn
150
151 is = nn*nn*(k -1) +nn*(j -1) +i
152 in = cs_loc(cs_loc(ie -1) +is)
153 iaz = 3*(in -1) + 1
154
155 !PAPER WITH BLANCA INITIAL CONDITIONS
156 u0(iaz+0) = + dsin(3.d0*pi*dt)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
157 u0(iaz+1) = - dsin(3.d0*pi*dt)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
158 u0(iaz+2) = - dsin(3.d0*pi*dt)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));
159
160 u1(iaz+0) = - 0.d0*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
161 u1(iaz+1) = + 0.d0*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
162 u1(iaz+2) = + 0.d0*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));
163
164 v1(iaz+0) = - 3.d0*pi*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
165 v1(iaz+1) = + 3.d0*pi*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
166 v1(iaz+2) = + 3.d0*pi*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));
167
168
169 !PLANE-WAVE SOLUTION INITIAL CONDITIONS
170! u0(iaz+0) = 0;
171! u0(iaz+1) = 0;
172
173! u1(iaz+0) = 0;
174! u1(iaz+1) = 0;
175
176! v1(iaz+0) = 0;
177! v1(iaz+1) = 0;
178
179! if (zs_loc(in) .ge. 50) then
180! u0(iaz+2) = -dcos(k1*xs_loc(in)+k2*ys_loc(in)+k3*zs_loc(in) + omega*dt);;
181! u1(iaz+2) = -dcos(k1*xs_loc(in)+k2*ys_loc(in)+k3*zs_loc(in));
182! v1(iaz+2) = dsin(k1*xs_loc(in)+k2*ys_loc(in)+k3*zs_loc(in));
183
184! else
185! u0(iaz+2) = 0;
186! u1(iaz+2) = 0;
187! v1(iaz+2) = 0;
188! endif
189
190
191
192
193 enddo
194 enddo
195 enddo
196
197 enddo
198
199
200 ! deallocate(val_data)
201