Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
data_py.f90
Go to the documentation of this file.
1 ! *********************************************************************
2 ! * Rocstar Simulation Suite *
3 ! * Copyright@2015, Illinois Rocstar LLC. All rights reserved. *
4 ! * *
5 ! * Illinois Rocstar LLC *
6 ! * Champaign, IL *
7 ! * www.illinoisrocstar.com *
8 ! * sales@illinoisrocstar.com *
9 ! * *
10 ! * License: See LICENSE file in top level of distribution package or *
11 ! * http://opensource.org/licenses/NCSA *
12 ! *********************************************************************
13 ! *********************************************************************
14 ! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
15 ! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
16 ! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
17 ! * NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
18 ! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
19 ! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
20 ! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
21 ! * USE OR OTHER DEALINGS WITH THE SOFTWARE. *
22 ! *********************************************************************
23 !**************************************************************************
24 MODULE data_py
25 
26  IMPLICIT NONE
27 
28  INTEGER, PARAMETER :: DBL = SELECTED_REAL_KIND(P=14,R=30)
29 !=========================================================================
30 
31  REAL *8, PARAMETER :: zero = 0.d0, &
32  one = 1.0d0, &
33  two = 2.0d0, &
34  three = 3.0d0, &
35  four = 4.0d0, &
36  eight = 8.0d0, &
37 
38  half = 0.5d0, &
39  quarter = 0.25d0, &
40  eighth = 0.125d0, &
41  third = one/three, &
42  two_thirds = two/three, &
43  root_two = 1.41421356237309515d0, &
44 
45  j2cal = 1.0d0/4.1868d0, &
46  kg2g = 1000.0d0, &
47  m2cm = 100.0d0, &
48  ru = 1.9859d0, &
49  mpa2atm = 9.86923266716d0, &
50  pa2atm = mpa2atm*1d-6, &
51  kgmc2gcc = kg2g/m2cm**3, &
52  j_m2cal_cm = j2cal/m2cm, &
53  j_msq2cal_cmsq = j2cal/m2cm**2, &
54  gcc2kgmc = one/kgmc2gcc, &
55  kcalmc2atm = 4186.8/101325.d0, &
56  j_kg2cal_g = j2cal/kg2g
57 !=======================================================================
58 ! global structure that contains combustion model parameters
59 
61 
62  REAL(DBL) :: a_p,n_p,Pref,Ac,eg_ru,ec_ru,alfac,C,lamg
63  REAL(DBL) :: Tstar0, To, Tignition, Tsurf, film_cons, lamc
64 
65  INTEGER :: comm, igrid, numx, ixsymm
66  INTEGER :: TABUSE
67  REAL(DBL) :: delt, xmax, beta, delz, x_surf_burn
68  REAL(DBL) :: delt_max, delz2inv, delzsqinv, dx1
69  REAL(DBL),POINTER :: x(:), z(:), zx(:), zxx(:)
70  REAL(DBL) :: P_range(2), rb_range(2), Tf_range(2)
71  CHARACTER*50 :: TABNAM
72 
73  TYPE(tabletype),POINTER :: TABLE
74 
75  END TYPE parameter_structure
76 
78 
79  REAL(DBL) :: P,rhoc,Qc,Qcprime,rb,Ts,To,lamc
80  REAL(DBL) :: fs,fsprime,Tstar,Ts0,Tslimit
81  REAL(DBL) :: c1,c2,c3,c4,c5,c6,alfa_eff,fx2
82  LOGICAL :: ignited
83 
84  END TYPE work_structure
85 
86  TYPE, public :: tabletype
87 !A- table dimensions
88  INTEGER :: nx_table,ny_table,nfield
89  LOGICAL :: spline
90 !B- table fields
91  REAL(DBL),POINTER :: tsurf00(:),press00(:), heatflux00(:,:)
92  REAL(DBL),POINTER :: rb00(:,:),fxsq00(:,:),Tgas00(:,:)
93  REAL(DBL),POINTER :: Tstd00(:),rstd00(:),alph00(:,:)
94 !C- Work Arrays
95  REAL(DBL),POINTER :: wrk1(:),wrk2(:),wrk3(:),wrk4(:),wrk5(:)
96  REAL(DBL),POINTER :: wrk6(:),wrk7(:),wrk8(:),wrk9(:),wrk10(:)
97 !D- Table parameters
98  REAL(DBL) :: alpha,chi,small_1,small_2
99  END TYPE tabletype
100 
101 CONTAINS
102 
103 !data interpolation subroutines
104 
105 !*****************************************************************************
106  SUBROUTINE polin2(bp,x1a,x2a,y12a,m,n,x1,x2,y,j,k)
107  implicit NONE
108  TYPE(parameter_structure),POINTER :: bp
109  INTEGER :: i,l,m,n,ixtrap,jsrt,jend
110  INTEGER,INTENT(INOUT) :: j,k
111  REAL(DBL) :: dy,x1,x2,y
112  REAL(DBL) :: x1a(m),x2a(n),y12a(m,n)
113  REAL(DBL) :: prod,difm,difp,a(4),del12
114  REAL(DBL), POINTER :: yntmp(:),ymtmp(:)
115 !-----------------------------------------------
116 
117 ! yntmp => Bp%TABLE%wrk1
118 ! ymtmp => Bp%TABLE%wrk2
119 
120  if(j+k > 0) goto 123
121 
122  j=0
123  ixtrap = 0
124  prod = 1.0
125  do while(prod .gt. 0)
126  j=j+1
127  difm=x1-x1a(j)
128  difp=x1-x1a(j+1)
129  prod = difm*difp
130  if(j .ge. m-1 .AND. prod .gt. 0)then
131  ixtrap = 1
132  prod = -1.0 !break
133  if(abs(x1-x1a(m)) .lt. abs(x1-x1a(1)) ) then
134  j = m-1
135  else
136  j = 2
137  endif
138  endif
139  enddo
140  k=0
141  ixtrap = 0
142  prod = 1.0
143  do while(prod .gt. 0)
144  k=k+1
145  difm=x2-x2a(k)
146  difp=x2-x2a(k+1)
147  prod = difm*difp
148  if(k .ge. n-1 .AND. prod .gt. 0)then
149  ixtrap = 1
150  prod = -1.0 !break
151  if(abs(x2-x2a(n)) .lt. abs(x2-x2a(1)) ) then
152  k = n-1
153  else
154  k = 2
155  endif
156  endif
157  enddo
158 
159 123 CONTINUE
160 
161  del12 = (x2a(k+1) - x2a(k))*(x1a(j+1) - x1a(j))
162 
163  a(1) = (x2a(k+1) - x2) * (x1a(j+1) - x1)
164  a(2) = (x2a(k+1) - x2) * (x1 - x1a(j))
165  a(3) = (x2 - x2a(k)) * (x1a(j+1) - x1)
166  a(4) = (x2 - x2a(k)) * (x1 - x1a(j))
167 
168  a = a/del12
169 
170 !use bilinear interpolation
171  y = a(1)*y12a(j,k) + a(2)*y12a(j+1,k) + a(3) * y12a(j,k+1) + a(4) * y12a(j+1,k+1)
172 
173 !---------------
174  RETURN
175  END SUBROUTINE polin2
176 !*****************************************************************************
177 
178 
179 
180 !*****************************************************************************
181 ! Spline interpoaltion routine.
182  SUBROUTINE polint(bp,xa,ya,n,x,y,dx)
183  implicit none
184  integer i,j,k,l,m,n,ixtrap
185  TYPE(parameter_structure),POINTER :: bp
186  REAL(DBL) :: x,y,prod,difm,difp,yo,der,dx
187  REAL(DBL) :: xa(n),ya(n)
188  REAL(DBL), POINTER :: h(:),alp(:)
189  REAL(DBL), POINTER :: c1(:),c2(:),c3(:)
190  REAL(DBL), POINTER :: v1(:),v2(:),v3(:)
191 !-----------------------------------------------
192 
193 
194  h => bp%TABLE%wrk3
195 
196  prod = 1.0
197  i=0
198  do i =1,n-1
199  h(i) = xa(i+1) - xa(i)
200  enddo
201 
202 !SPLINE
203  if(bp%TABLE%spline) then
204  alp => bp%TABLE%wrk4
205  c1 => bp%TABLE%wrk5
206  c2 => bp%TABLE%wrk6
207  c3 => bp%TABLE%wrk7
208  v1 => bp%TABLE%wrk8
209  v2 => bp%TABLE%wrk9
210  v3 => bp%TABLE%wrk10
211  do i =2,n-1
212  j = i-1
213  alp(i) = 3.0/h(i)*(ya(i+1) -ya(i)) &
214  - 3.0/h(j)*(ya(j+1) -ya(j))
215  enddo
216  v1(1) = 1.d0
217  v2(1) = 0.d0
218  v3(1) = 0.d0
219  do i = 2,n-1
220  v1(i) = 2.0d0*(xa(i+1) -xa(i-1)) - h(i-1)*v2(i-1)
221  v2(i) = h(i)/v1(i)
222  v3(i) = ( alp(i) - h(i-1)*v3(i-1) ) / v1(i)
223  enddo
224 
225  v1(n) = 1.d0
226  v3(n) = 0.d0
227  c2(n) = 0.d0
228  do i = n-1,1,-1
229  c2(i) = v3(i) - v2(i)* c2(i+1)
230  c1(i) = (ya(i+1) -ya(i))/h(i) - &
231  h(i) * ( c2(i+1)+2.0*c2(i) ) / 3.0d0
232  c3(i) = (c2(i+1)-c2(i))/(3.0d0*h(i))
233  enddo
234  endif
235 
236  i=0
237  ixtrap = 0
238  prod = 1.0
239  do while(prod .gt. 0)
240  i=i+1
241  difm=x-xa(i)
242  difp=x-xa(i+1)
243  prod = difm*difp
244  if(i .ge. n-1 .AND. prod .gt. 0)then
245 ! write(*,101)'FAILED SEARCH',x,(xa(m),m=1,n)
246  ixtrap = 1
247  prod = -1.0
248  endif
249  enddo
250 !101 format(2x,a,1p100e10.2)
251 
252 
253  if(bp%TABLE%spline) then
254 ! SPLINE
255  y = ya(i) + c1(i)*(x-xa(i)) + c2(i)*(x-xa(i))**2&
256  + c3(i)*(x-xa(i))**3
257  else
258 ! USE first order interpolation to avoid problem w\ oscillatory solutions
259  y = ya(i) + (ya(i+1)-ya(i))/h(i) * (x-xa(i))
260  endif
261 
262  if(ixtrap .eq. 1) then
263  if(abs(x-xa(n)) .lt. abs(x-xa(1)) ) then
264  yo = ya(n)
265  der = (ya(n)-ya(n-1))/h(n-1)
266  dx = x-xa(n)
267  else
268  yo = ya(1)
269  der = (ya(2)-ya(1))/h(1)
270  dx = x-xa(1)
271  endif
272  y = yo+der*dx
273  endif
274 
275  return
276  END SUBROUTINE polint
277 !*************************************************
278 
279 
280 !=======================================================================
281 
282 END MODULE data_py
283 !************************************************************************
284 
285 
286 
287 
288 
289 
FT m(int i, int j) const
const NT & d
j indices k indices k
Definition: Indexing.h:6
void int int REAL REAL * y
Definition: read.cpp:74
NT dx
subroutine polint(bp, xa, ya, n, x, y, dx)
Definition: data_py.f90:182
blockLoc i
Definition: read.cpp:79
subroutine polin2(bp, x1a, x2a, y12a, m, n, x1, x2, y, j, k)
Definition: data_py.f90:106
void int int REAL * x
Definition: read.cpp:74
const NT & n
j indices j
Definition: Indexing.h:6
NT dy
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE knode jend
RT a() const
Definition: Line_2.h:140