MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
mod_DFT_CHAIN_d.F90
Go to the documentation of this file.
1 
5 
6 
7 
8 ! Generated by TAPENADE (INRIA, Tropics team)
9 ! Tapenade 3.10 (r5498) - 20 Jan 2015 09:48
10 !
11 MODULE mod_dft_chain_d
12  IMPLICIT NONE
13  PRIVATE
14  PUBLIC chain_aux_d
15  PUBLIC chain_dfdrho_d
16 
17 CONTAINS
18 ! Differentiation of chain_aux in forward (tangent) mode:
19 ! variations of useful results: rhobar lambda
20 ! with respect to varying inputs: rhop
21  SUBROUTINE chain_aux_d(rhop, rhopd, rhobar, rhobard, lambda, lambdad, &
22 & user)
23  USE basic_variables, ONLY : ncomp
24  USE eos_variables, Only: dhs,rho
25  USE mod_dft, ONLY : zp, dzp, fa
26 
27  !PETSc module
28  Use petscmanagement
29 
30  IMPLICIT NONE
31 
32 #include <finclude/petscsys.h>
33 
34 
35 !passed
36 type(userctx) :: user
37 petscscalar :: rhop(ncomp,user%gxs:user%gxe)
38 petscscalar :: rhopd(ncomp,user%gxs:user%gxe)
39 REAL,INTENT(OUT) :: rhobar(user%gxs:user%gxe,ncomp)
40 REAL,INTENT(OUT) :: rhobard(user%gxs:user%gxe,ncomp)
41 REAL,INTENT(OUT) :: lambda(user%gxs:user%gxe,ncomp)
42 REAL,INTENT(OUT) :: lambdad(user%gxs:user%gxe,ncomp)
43 
44 !local
45  INTEGER :: i, j, k
46  REAL :: dhsk
47  INTEGER :: fak, n
48  REAL :: zz, dz, xlo, xhi, integral_lamb, integral_rb
49  REAL :: integral_lambd, integral_rbd
50  INTEGER, parameter :: nmax = 800
51  REAL, DIMENSION(NMAX) :: x_int, lamb_int, rb_int
52  REAL, DIMENSION(NMAX) :: lamb_intd, rb_intd
53  REAL, DIMENSION(NMAX) :: y2_lamb, y2_rb
54  REAL, DIMENSION(NMAX) :: y2_lambd, y2_rbd
55  REAL :: rhopjk, rhopjp1k
56  REAL :: rhopjkd, rhopjp1kd
57  INTRINSIC epsilon
58  REAL :: result1
59 
60 
61 
62  rhobard = 0.0
63  lambdad = 0.0
64  y2_rbd = 0.0
65  y2_lambd = 0.0
66 !fak = maxval(fa(1:ncomp))
67  DO k=1,ncomp
68  dhsk = dhs(k)
69  fak = fa(k) + 1
70 
71  Do i = user%xs-fak,user%xe+fak !lambda und rhobar werden bis i+-sig gebraucht -> Schleife bis +- fa
72  n = 1
73  x_int = 0.0
74  lamb_int = 0.0
75  rb_int = 0.0
76  rb_intd = 0.0
77  lamb_intd = 0.0
78 !um lambda bei i zu berechnen, muss bis +- sig um i integriert werden -> Schleife bis +- fa
79  DO j=i-fak,i+fak
80  rhopjkd = rhopd(k, j)
81  rhopjk = rhop(k, j)
82  rhopjp1kd = rhopd(k, j+1)
83  rhopjp1k = rhop(k, j+1)
84  IF (zp(i) - zp(j+1) .LT. dhsk .AND. zp(i) - zp(j) .GE. dhsk) &
85 & THEN
86 !the position of j+1 is already within i-d while j is still outside this range in this case, the integration steplength (dz) is j
87 !ust the distance, which j+1 overlaps with i-d and what is integrated is the interpolated value of the integrand
88 !here always n=1!
89  IF (n .NE. 1) THEN
90  GOTO 100
91  ELSE
92 !distance between grid points j and i
93  zz = zp(j) - zp(i)
94 !the part of the intervall between zp(j) and zp(j+1) which is already within i-d
95  dz = zp(j+1) - (zp(i)-dhsk)
96 !if(dz < epsilon(dz)) dz = epsilon(dz) !bei unguenstiger Kombination von sig und ngrid kann dz unter Machinengenauigkeit epsilon
97 !liegen, dann ist x(2) = x(1) + dz = x(1) -> das fuehrt zu Abbruch in Spline Interpolation
98 !array containing x-values for spline integration
99  x_int(n) = 0.0
100 !lineare interpolation analog zum FMT Teil
101  lamb_intd(n) = rhopjkd + (dzp-dz)*(rhopjp1kd-rhopjkd)/dzp
102  lamb_int(n) = rhopjk + (rhopjp1k-rhopjk)/dzp*(dzp-dz)
103 !erklärung analog wie bei n3_int in FMT Teil
104  rb_intd(n) = 0.0
105  rb_int(n) = 0.0
106  END IF
107  ELSE IF (zp(j) .GT. zp(i) - dhsk .AND. zp(j) .LE. zp(i) + dhsk&
108 & ) THEN
109 !grid point j within i+-d
110  n = n + 1
111 !first time in this If condition, dz is stil the old value from above!
112  x_int(n) = x_int(n-1) + dz
113  zz = zp(j) - zp(i)
114  dz = dzp
115  lamb_intd(n) = rhopjkd
116  lamb_int(n) = rhopjk
117  rb_intd(n) = (dhsk**2-zz**2)*rhopjkd
118  rb_int(n) = rhopjk*(dhsk**2-zz**2)
119  IF (zp(j) .LT. zp(i) + dhsk .AND. zp(j+1) .GE. zp(i) + dhsk&
120 & ) THEN
121 !zp(j) is still within zp(i)+d but zp(j+1) is already outside zp(i)+d
122  dz = zp(i) + dhsk - zp(j)
123 !If(dz <= epsilon(dz)) exit !wie oben, kann auch hier bei ungluecklicher Wahl von sig und ngrid dz < eps werden und somit x(n)
124 != x(n-1) -> Abbruch in Spline interpolation. Dann einfach ngrid aendern!
125  zz = zp(j) - zp(i)
126  n = n + 1
127  x_int(n) = x_int(n-1) + dz
128  rb_intd(n) = 0.0
129  rb_int(n) = 0.0
130  lamb_intd(n) = rhopjkd + dz*(rhopjp1kd-rhopjkd)/dzp
131  lamb_int(n) = rhopjk + (rhopjp1k-rhopjk)/dzp*dz
132 ! If(x_int(n) == x_int(n-1)) Then
133 ! n = n - 1
134 ! exit
135 ! End If
136  END IF
137  END IF
138  END DO
139 !spline integration
140  xlo = x_int(1)
141  xhi = x_int(n)
142  CALL spline_d(x_int, lamb_int, lamb_intd, n, 1.e30, 1.e30, &
143 & y2_lamb, y2_lambd)
144  CALL spline_d(x_int, rb_int, rb_intd, n, 1.e30, 1.e30, y2_rb, &
145 & y2_rbd)
146  CALL splint_integral_d(x_int, lamb_int, lamb_intd, y2_lamb, &
147 & y2_lambd, n, xlo, xhi, integral_lamb, &
148 & integral_lambd)
149  CALL splint_integral_d(x_int, rb_int, rb_intd, y2_rb, y2_rbd, n&
150 & , xlo, xhi, integral_rb, integral_rbd)
151  lambdad(i, k) = 0.5*integral_lambd/dhsk
152  lambda(i, k) = 0.5*integral_lamb/dhsk
153  rhobard(i, k) = 0.75*integral_rbd/dhsk**3
154  rhobar(i, k) = 0.75*integral_rb/dhsk**3
155  result1 = epsilon(dz)
156  IF (lambda(i, k) .LT. result1) lambda(i, k) = epsilon(dz)
157  END DO
158  END DO
159  GOTO 110
160  100 stop
161  110 CONTINUE
162  END SUBROUTINE chain_aux_d
163 
164 
165 
166 
167 ! Differentiation of chain_dfdrho in forward (tangent) mode:
168 ! variations of useful results: df_drho_chain rhop
169 ! with respect to varying inputs: rhobar df_drho_chain lambda
170 ! rhop
171  SUBROUTINE chain_dfdrho_d(i, rhop, rhopd, lambda, lambdad, rhobar, &
172 & rhobard, df_drho_chain, df_drho_chaind, user)
173  USE basic_variables, ONLY : ncomp,parame
174  USE eos_variables, Only: dhs
175  USE mod_dft, ONLY : zp, dzp, fa
176 
177  !PETSc module
178  Use petscmanagement
179  IMPLICIT NONE
180 
181 #include <finclude/petscsys.h>
182 
183 
184 !passed
185 INTEGER, INTENT(IN) :: i !the grid point to calculate dFdrho at
186 type(userctx) :: user
187 petscscalar :: rhop(ncomp,user%gxs:user%gxe)
188 petscscalar :: rhopd(ncomp,user%gxs:user%gxe)
189 REAL,INTENT(IN) :: rhobar(user%gxs:user%gxe,ncomp)
190 REAL,INTENT(IN) :: rhobard(user%gxs:user%gxe,ncomp)
191 REAL,INTENT(IN) :: lambda(user%gxs:user%gxe,ncomp)
192 REAL,INTENT(IN) :: lambdad(user%gxs:user%gxe,ncomp)
193 
194 
195 
196 REAL,INTENT(OUT) :: df_drho_chain(user%xs:user%xe,ncomp)
197 REAL,INTENT(OUT) :: df_drho_chaind(user%xs:user%xe,ncomp)
198 
199 
200 !local
201  INTEGER :: j, k, n
202  REAL :: dhsk
203  INTEGER :: fak
204  REAL :: rhopjk, rhopjp1k, logrho, xlo, xhi
205  REAL :: rhopjkd, rhopjp1kd
206  REAL :: ycorr(ncomp), dlny(ncomp, ncomp)
207  REAL :: ycorrd(ncomp), dlnyd(ncomp, ncomp)
208  INTEGER, parameter :: nmax = 800
209  REAL, DIMENSION(NMAX) :: x_int, int_1, int_2
210  REAL, DIMENSION(NMAX) :: int_1d, int_2d
211  REAL, DIMENSION(NMAX) :: y2_1, y2_2
212  REAL, DIMENSION(NMAX) :: y2_1d, y2_2d
213  REAL :: dz, zz, integral_1, integral_2
214  REAL :: integral_1d, integral_2d
215  REAL :: lamy
216  INTRINSIC sum
217  INTRINSIC epsilon
218  INTRINSIC log
219  REAL, DIMENSION(ncomp) :: arg1
220  REAL, DIMENSION(ncomp) :: arg1d
221  REAL :: result1
222  REAL :: arg10
223  REAL :: arg10d
224 
225 
226 
227 
228 
229  CALL cavity_mix_d(rhobar(i, 1:ncomp), rhobard(i, 1:ncomp), ycorr, &
230 & ycorrd, dlny, dlnyd)
231  y2_1d = 0.0
232  y2_2d = 0.0
233  DO k=1,ncomp
234  fak = fa(k)
235  dhsk = dhs(k)
236  n = 1
237  x_int = 0.0
238  int_1 = 0.0
239  int_2 = 0.0
240  int_1d = 0.0
241  int_2d = 0.0
242 !es muss bis +- sig um i integriert werden
243  DO j=i-fak,i+fak
244  rhopjkd = rhopd(k, j)
245  rhopjk = rhop(k, j)
246  rhopjp1kd = rhopd(k, j+1)
247  rhopjp1k = rhop(k, j+1)
248  IF (zp(i) - zp(j+1) .LT. dhsk .AND. zp(i) - zp(j) .GE. dhsk) &
249 & THEN
250 !the position of j+1 is already within i-d while j is still outside this range in this case, the integration steplength (dz) is j
251 !ust the distance, which j+1 overlaps with i-d and what is integrated is the interpolated value of the integrand
252 !here always n=1!
253  IF (n .NE. 1) THEN
254  GOTO 100
255  ELSE
256 !distance between grid points j and i
257  zz = zp(j) - zp(i)
258 !the part of the intervall between zp(j) and zp(j+1) which is already within i-d
259  dz = zp(j+1) - (zp(i)-dhsk)
260 !if(dz < epsilon(dz)) dz = epsilon(dz) !bei unguenstiger Kombination von sig und ngrid kann dz unter Machinengenauigkeit epsilon
261 !liegen, dann ist x(2) = x(1) + dz = x(1) -> das fuehrt zu Abbruch in Spline Interpolation
262 !array containing x-values for spline integration
263  x_int(n) = 0.0
264 !erklärung analog wie bei n3_int in FMT Teil
265  int_1d(n) = 0.0
266  int_1(n) = 0.0
267 !lineare interpolation analog zum FMT Teil
268  int_2d(n) = rhopjkd*lambda(j, k) + rhopjk*lambdad(j, k) + (&
269 & dzp-dz)*(rhopjp1kd*lambda(j+1, k)+rhopjp1k*lambdad(j+1, k)&
270 & -rhopjkd*lambda(j, k)-rhopjk*lambdad(j, k))/dzp
271  int_2(n) = rhopjk*lambda(j, k) + (rhopjp1k*lambda(j+1, k)-&
272 & rhopjk*lambda(j, k))/dzp*(dzp-dz)
273  END IF
274  ELSE IF (zp(j) .GT. zp(i) - dhsk .AND. zp(j) .LE. zp(i) + dhsk) &
275 & THEN
276 !grid point j within i+-d
277  n = n + 1
278 !first time in this If condition, dz is stil the old value from above!
279  x_int(n) = x_int(n-1) + dz
280  zz = zp(j) - zp(i)
281  dz = dzp
282  arg1d(:) = (parame(1:ncomp,1)-1.0)*(rhopd(1:ncomp, j)*dlny(1:ncomp&
283 & , k)+rhop(1:ncomp, j)*dlnyd(1:ncomp, k))
284  arg1(:) = (parame(1:ncomp,1)-1.0)*rhop(1:ncomp, j)*dlny(1:ncomp, k&
285 & )
286  int_1d(n) = (dhsk**2-zz**2)*0.75*sum(arg1d(:))/dhsk**3
287  int_1(n) = sum(arg1(:))*0.75/dhsk**3*(dhsk**2-zz**2)
288  int_2d(n) = (rhopjkd*lambda(j, k)-rhopjk*lambdad(j, k))/lambda&
289 & (j, k)**2
290  int_2(n) = rhopjk/lambda(j, k)
291  IF (zp(j) .LT. zp(i) + dhsk .AND. zp(j+1) .GE. zp(i) + dhsk) &
292 & THEN
293 !zp(j) is still within zp(i)+d but zp(j+1) is already outside zp(i)+d
294  dz = zp(i) + dhsk - zp(j)
295 !If(dz <= epsilon(dz)) exit !wie oben, kann auch hier bei ungluecklicher Wahl von sig und ngrid dz < eps werden und somit x(n)
296 != x(n-1) -> Abbruch in Spline interpolation. Dann einfach ngrid aendern!
297  zz = zp(j) - zp(i)
298  n = n + 1
299  x_int(n) = x_int(n-1) + dz
300  int_1d(n) = 0.0
301  int_1(n) = 0.0
302  int_2d(n) = rhopjkd*lambda(j, k) + rhopjk*lambdad(j, k) + dz&
303 & *(rhopjp1kd*lambda(j+1, k)+rhopjp1k*lambdad(j+1, k)-&
304 & rhopjkd*lambda(j, k)-rhopjk*lambdad(j, k))/dzp
305  int_2(n) = rhopjk*lambda(j, k) + (rhopjp1k*lambda(j+1, k)-&
306 & rhopjk*lambda(j, k))/dzp*dz
307 ! If(x_int(n) == x_int(n-1)) Then
308 ! n = n - 1
309 ! exit
310 ! End If
311  END IF
312  END IF
313  END DO
314 !Spline Integration
315  xlo = x_int(1)
316  xhi = x_int(n)
317  CALL spline_d(x_int, int_1, int_1d, n, 1.e30, 1.e30, y2_1, y2_1d)
318  CALL splint_integral_d(x_int, int_1, int_1d, y2_1, y2_1d, n, xlo, &
319 & xhi, integral_1, integral_1d)
320  CALL spline_d(x_int, int_2, int_2d, n, 1.e30, 1.e30, y2_2, y2_2d)
321  CALL splint_integral_d(x_int, int_2, int_2d, y2_2, y2_2d, n, xlo, &
322 & xhi, integral_2, integral_2d)
323  result1 = epsilon(dz)
324  IF (rhop(k, i) .LT. result1) THEN
325  rhop = epsilon(dz)
326  rhopd = 0.0
327  END IF
328  arg10d = ycorrd(k)*lambda(i, k) + ycorr(k)*lambdad(i, k)
329  arg10 = ycorr(k)*lambda(i, k)
330  df_drho_chaind(i, k) = (parame(k,1)-1.0)*rhopd(k, i)/rhop(k, i) - (&
331 & parame(k,1)-1.0)*(arg10d/arg10+0.5*integral_2d/dhsk) - integral_1d
332  df_drho_chain(i, k) = (parame(k,1)-1.0)*log(rhop(k, i)) - (parame(k,1)-1.0&
333 & )*(log(arg10)-1.0+0.5*integral_2/dhsk) - integral_1
334  END DO
335  GOTO 110
336  100 stop
337  110 CONTINUE
338  END SUBROUTINE chain_dfdrho_d
339 
340 
341 
342 ! Differentiation of cavity_mix in forward (tangent) mode:
343 ! variations of useful results: ycorr dlnydr
344 ! with respect to varying inputs: rhoi
345 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
346 !
347 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
348 !
349  SUBROUTINE cavity_mix_d(rhoi, rhoid, ycorr, ycorrd, dlnydr, dlnydrd)
350 !
351  USE parameters, ONLY : pi
352  USE basic_variables, ONLY : ncomp,parame
353  USE eos_variables, Only: dhs
354  IMPLICIT NONE
355 !
356 ! ----------------------------------------------------------------------
357  REAL, INTENT(IN) :: rhoi(ncomp)
358  REAL, INTENT(IN) :: rhoid(ncomp)
359  REAL, INTENT(OUT) :: ycorr(ncomp)
360  REAL, INTENT(OUT) :: ycorrd(ncomp)
361 ! this is: d( ln( yij ) ) / d( rho(k) ) used with i=j
362  REAL, INTENT(OUT) :: dlnydr(ncomp, ncomp)
363  REAL, INTENT(OUT) :: dlnydrd(ncomp, ncomp)
364 !
365 ! ----------------------------------------------------------------------
366  INTEGER :: i, j, k
367  REAL :: z0, z1, z2, z3, zms, z1_rk, z2_rk, z3_rk
368  REAL :: z2d, z3d, zmsd
369  REAL, DIMENSION(ncomp, ncomp) :: dij_ab, gij, gij_rk
370  REAL, DIMENSION(ncomp, ncomp) :: gijd, gij_rkd
371  INTRINSIC sum
372  REAL, DIMENSION(ncomp) :: arg1
373  REAL, DIMENSION(ncomp) :: arg1d
374 ! ----------------------------------------------------------------------
375  arg1(:) = rhoi(1:ncomp)*parame(1:ncomp,1)
376  z0 = pi/6.0*sum(arg1(:))
377  arg1(:) = rhoi(1:ncomp)*parame(1:ncomp,1)*dhs(1:ncomp)
378  z1 = pi/6.0*sum(arg1(:))
379  arg1d(:) = parame(1:ncomp,1)*dhs(1:ncomp)**2*rhoid(1:ncomp)
380  arg1(:) = rhoi(1:ncomp)*parame(1:ncomp,1)*dhs(1:ncomp)**2
381  z2d = pi*sum(arg1d(:))/6.0
382  z2 = pi/6.0*sum(arg1(:))
383  arg1d(:) = parame(1:ncomp,1)*dhs(1:ncomp)**3*rhoid(1:ncomp)
384  arg1(:) = rhoi(1:ncomp)*parame(1:ncomp,1)*dhs(1:ncomp)**3
385  z3d = pi*sum(arg1d(:))/6.0
386  z3 = pi/6.0*sum(arg1(:))
387  zmsd = -z3d
388  zms = 1.0 - z3
389  DO i=1,ncomp
390  DO j=1,ncomp
391  dij_ab(i, j) = dhs(i)*dhs(j)/(dhs(i)+dhs(j))
392  END DO
393  END DO
394  ycorrd = 0.0
395  dlnydrd = 0.0
396  gij_rkd = 0.0
397  gijd = 0.0
398  DO k=1,ncomp
399  DO i=1,ncomp
400  z1_rk = pi/6.0*parame(k,1)*dhs(k)
401  z2_rk = pi/6.0*parame(k,1)*dhs(k)*dhs(k)
402  z3_rk = pi/6.0*parame(k,1)*dhs(k)**3
403 !DO j = 1, ncomp
404  j = i
405  gijd(i, j) = ((3.0*dij_ab(i, j)*z2d*zms-3.0*dij_ab(i, j)*z2*zmsd&
406 & )/zms-3.0*dij_ab(i, j)*z2*zmsd/zms)/zms**2 - zmsd/zms**2 + (&
407 & 2.0*2*dij_ab(i, j)**2*z2*z2d*zms**3-2.0*dij_ab(i, j)**2*z2**2*&
408 & 3*zms**2*zmsd)/(zms**3)**2
409  gij(i, j) = 1.0/zms + 3.0*dij_ab(i, j)*z2/zms/zms + 2.0*(dij_ab(&
410 & i, j)*z2)**2/zms**3
411 !dgijdz(i,j)= 1.0/zms/zms + 3.0*dij_ab(i,j)*z2*(1.0+z3)/z3/zms**3 &
412 ! + (dij_ab(i,j)*z2/zms/zms)**2 *(4.0+2.0*z3)/z3
413  gij_rkd(i, j) = (-2)*(z3_rk*zmsd/zms)/zms**2 + ((3.0*dij_ab(i, j&
414 & )*(2.0*z3_rk*z2d*zms-2.0*z2*z3_rk*zmsd)/zms-3.0*dij_ab(i, j)*(&
415 & z2_rk+2.0*z2*z3_rk/zms)*zmsd)/zms-3.0*dij_ab(i, j)*(z2_rk+2.0*&
416 & z2*z3_rk/zms)*zmsd/zms)/zms**2 + (dij_ab(i, j)**2*z2d*zms**3-&
417 & dij_ab(i, j)**2*z2*3*zms**2*zmsd)*(4.0*z2_rk+6.0*z2*z3_rk/zms)&
418 & /zms**6 + dij_ab(i, j)**2*z2*(6.0*z3_rk*z2d*zms-6.0*z2*z3_rk*&
419 & zmsd)/zms**5
420  gij_rk(i, j) = z3_rk/zms/zms + 3.0*dij_ab(i, j)*(z2_rk+2.0*z2*&
421 & z3_rk/zms)/zms/zms + dij_ab(i, j)**2*z2/zms**3*(4.0*z2_rk+6.0*&
422 & z2*z3_rk/zms)
423 !END DO
424  ycorrd(i) = gijd(i, i)
425  ycorr(i) = gij(i, i)
426  dlnydrd(i, k) = (gij_rkd(i, i)*gij(i, i)-gij_rk(i, i)*gijd(i, i)&
427 & )/gij(i, i)**2
428  dlnydr(i, k) = gij_rk(i, i)/gij(i, i)
429  END DO
430  END DO
431  END SUBROUTINE cavity_mix_d
432 
433 END MODULE mod_dft_chain_d
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
lambda
double lambda
Latent heat of blowing agent, J/kg.
Definition: experimentalInputs.h:118
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
petscmanagement
In this module, the application context is defined.
Definition: mod_PETSc.F90:7
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29