SPEED
GET_INITIAL_CONDITIONS.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
19
40!
41!**************************************************************************************************
42
43
44 subroutine get_initial_conditions(nnod_loc, u0, u1, v1, ne_loc, cs_loc, cs_nnz_loc,&
45 nm,prop_mat, sdeg_mat, loc_n_num, &
46 xs_loc,ys_loc,zs_loc,&
47 mpi_id,dt)
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
202 end subroutine get_initial_conditions
get_initial_conditions
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.
Definition GET_INITIAL_CONDITIONS.f90:48