Context Navigation

Heat.F90

main

Last change on this file was aad342c, checked in by Stephen Siegel <siegel@…>, 3 years ago

Performing huge refactor to incorporate ABC, GMC, and SARL into CIVL repo and use Java modules.

git-svn-id: svn://vsl.cis.udel.edu/civl/trunk@5664 fb995dde-84ed-4084-dfe6-e5aef3e2452c

Property mode set to 100644

File size: 7.3 KB

Line
1	!!***if source/physics/sourceTerms/Heat/HeatMain/Neutrino/Heat
2	!!
3	!! NAME
4	!!
5	!! Heat
6	!!
7	!!
8	!! SYNOPSIS
9	!!
10	!! call Heat (integer(IN) :: blockCount,
11	!! integer(IN) :: blockList(blockCount),
12	!! real(IN) :: dt,
13	!! real(IN) :: time)
14	!!
15	!!
16	!!
17	!! DESCRIPTION
18	!!
19	!! Calculates local heating and cooling due to neutrinos
20	!! according to the approach of Murphy & Burrows (2008, ApJ, 688, 1159).
21	!!
22	!! NOTES
23	!! The FLASH implementation of this method is described in
24	!! S.M. Couch (2013, ApJ, 765, 29). Citation of this latter
25	!! reference is appreciated if this unit is used in published work.
26	!!
27	!! ARGUMENTS
28	!!
29	!! blockCount : number of blocks to operate on
30	!! blockList : list of blocks to operate on
31	!! dt : current timestep
32	!! time : current time
33	!!
34	!!***
35	!
36	!==============================================================================
37	!
38	subroutine Heat(blockCount, blockList, dt, time)
39
40	#include "Flash.h"
41	#include "constants.h"
42
43
44	#ifdef _CIVL
45	#define solnData(v,i,j,k) data(v,i,j,k,bId)
46	use Grid_interface, ONLY : &
47	Grid_getBlkIndexLimits, Grid_getCellCoords
48
49	!$ civl $input
50	real,dimension(:,:,:,:,:) :: data_in
51	!$ civl $output
52	real,dimension(:,:,:,:,:) :: data_out
53	!$ civl $input
54	logical :: useHeat
55	!$ civl $input
56	real :: ht_Lneut, ht_Tneut, ht_bounceTime, ht_postBounce
57
58	#else
59	use Heat_data, ONLY : useHeat, ht_Lneut, ht_Tneut, ht_bounceTime, ht_postBounce, &
60	ht_useEntr
61	use Grid_interface, ONLY : Grid_getBlkPtr, Grid_releaseBlkPtr, &
62	Grid_getBlkIndexLimits, Grid_getCellCoords
63	use Deleptonize_interface, ONLY : Deleptonize_getBounce
64	use Eos_interface, ONLY : Eos_wrapped, Eos_nucOneZone
65	#endif
66
67	implicit none
68
69	integer,intent(IN) :: blockCount
70	integer,dimension(blockCount),intent(IN)::blockList
71	real,intent(IN) :: dt,time
72
73	real, pointer, dimension(:,:,:,:) :: solnData
74	integer, dimension(LOW:HIGH,MDIM) :: blkLimits, blkLimitsGC
75
76	integer :: blockID
77	integer :: i,j,k,n
78	integer,dimension(MDIM) :: dimSize
79	real,allocatable, dimension(:) :: xCenter, yCenter, zCenter
80	real :: xx, yy, zz
81	logical :: gcell = .true.
82	real :: radius, dEneut, ek, tauNu
83	real, parameter :: MeVtoK = 1.16045221d10
84
85	real :: bounceTime, temp
86
87	real :: xDens, xTemp, xEner, xEntr, xYe, outYe
88	real :: del_ye, del_entr
89	real :: abar, zbar, sumY, Ye0, Ye, dXneut, abarInv
90	real :: xPres, mu_nu, xCs2
91	real :: xXp, xXn, xXa, xXh,xdedt,xdpderho, tmp
92	logical :: threadBlockList
93	logical :: eosCall2
94
95	#ifndef _CIVL
96
97	#ifdef ST_THREAD_BLOCK_LIST
98	threadBlockList = .true.
99
100	#ifdef ST_THREAD_WITHIN_BLOCK
101	call Driver_abortFlash("Cannot include both threading strategies")
102	#endif
103
104	#else
105	threadBlockList = .false.
106	#endif
107
108	if (.NOT. useHeat) return
109
110	if (.NOT. ht_postBounce) &
111	call Deleptonize_getBounce(ht_postBounce, ht_bounceTime)
112
113	if (.NOT. ht_postBounce) return
114
115	!$omp parallel if(threadBlockList) &
116	!$omp private(n,blockID,k,j,i,solnData,dimSize,zCenter,yCenter, &
117	!$omp xCenter,radius,dEneut,ek,blkLimits,blkLimitsGC,tmp,xDens,xTemp, &
118	!$omp xYe,xEner,xPres,xEntr,xdedt,xdpderho,mu_nu,xXp,xXn,xXa,xXh,xCs2,temp) &
119	!$omp shared(blockCount,blockList,ht_Lneut,ht_Tneut,dt,ht_postBounce, &
120	!$omp gcell,eosCall2)
121
122	#else
123
124	!$omp parallel &
125	!$omp private(n,blockID,k,j,i,dimSize,zCenter,yCenter, &
126	!$omp xCenter,radius,dEneut,ek,blkLimits,blkLimitsGC,tmp,xDens,xTemp, &
127	!$omp xYe,xEner,xPres,xEntr,xdedt,xdpderho,mu_nu,xXp,xXn,xXa,xXh,xCs2,temp) &
128	!$omp shared(blockCount,blockList,ht_Lneut,ht_Tneut,dt,ht_postBounce, &
129	!$omp data, &
130	!$omp gcell,eosCall2)
131
132	! memcpy data_in to data
133
134	#endif
135
136	!$omp do schedule(static)
137	do n = 1, blockCount
138	blockID = blockList(n)
139
140	call Grid_getBlkIndexLimits(blockID,blkLimits,blkLimitsGC)
141
142	#ifndef _CIVL
143	call Grid_getBlkPtr(blockID,solnData)
144	#endif
145
146
147	dimSize(:)=blkLimitsGC(HIGH,:)-blkLimitsGC(LOW,:)+1
148	if (NDIM > 2)then
149	allocate(zCenter(dimSize(KAXIS)))
150	call Grid_getCellCoords(KAXIS,blockID,&
151	CENTER,gcell,zCenter,dimSize(KAXIS))
152	end if
153	if (NDIM > 1)then
154	allocate(yCenter(dimSize(JAXIS)))
155	call Grid_getCellCoords(JAXIS,blockID,&
156	CENTER,gcell,yCenter,dimSize(JAXIS))
157	end if
158
159	allocate(xCenter(dimSize(IAXIS)))
160	call Grid_getCellCoords(IAXIS,blockID,&
161	CENTER,gcell,xCenter,dimSize(IAXIS))
162
163	#ifdef DELE_VAR
164	solnData(DELE_VAR,:,:,:) = 0.0
165	#endif
166
167	do k = blkLimits(LOW,KAXIS), blkLimits(HIGH,KAXIS)
168	do j = blkLimits(LOW,JAXIS), blkLimits(HIGH,JAXIS)
169	do i = blkLimits(LOW,IAXIS), blkLimits(HIGH,IAXIS)
170
171	#ifndef _CIVL
172	dEneut = 0.
173
174	xDens = solnData(DENS_VAR,i,j,k)
175	xTemp = solnData(TEMP_VAR,i,j,k)
176	xYe = solnData(YE_MSCALAR,i,j,k)
177	xEner = solnData(EINT_VAR,i,j,k)
178
179	radius = (xCenter(i))**2
180	if (NDIM > 1) then
181	radius = radius + (yCenter(j))**2
182	endif
183	if (NDIM > 2) then
184	radius = radius + zCenter(k)**2
185	endif
186	radius = sqrt(radius)
187
188	! Calculate heating
189	dEneut = 1.544e20 * (ht_Lneut/1.e52) * (1.e7 / radius)*2 (ht_Tneut / 4.)**2
190
191	! Now subtract cooling
192	dEneut = dEneut - 1.399e20 * (xTemp / (2.MeVtoK))*6
193
194	dEneut = dEneut * exp(-tauNu(xDens))
195
196	! Now call Eos to get Yp and Yn
197	call Eos_nucOneZone(xDens,xTemp,xYe,xEner,xPres,xEntr,&
198	tmp,tmp,tmp,tmp,tmp,tmp,xXp,16,MODE_DENS_TEMP)
199	call Eos_nucOneZone(xDens,xTemp,xYe,xEner,xPres,xEntr,&
200	tmp,tmp,tmp,tmp,tmp,tmp,xXn,15,MODE_DENS_TEMP)
201
202	dEneut = dEneut * (xXp + xXn)
203
204	solnData(EINT_VAR,i,j,k) = solnData(EINT_VAR,i,j,k) + dEneut*dt
205
206	ek = 0.5e0(solnData(VELX_VAR,i,j,k)*2 + &
207	solnData(VELY_VAR,i,j,k)**2 + &
208	solnData(VELZ_VAR,i,j,k)**2)
209
210	solnData(ENER_VAR,i,j,k) = solnData(EINT_VAR,i,j,k) + ek
211
212	! Now store any auxilliary variables
213	#ifdef DELE_VAR
214	solnData(DELE_VAR,i,j,k) = dEneut
215	#endif
216	#ifdef TAUN_VAR
217	solnData(TAUN_VAR,i,j,k) = tauNu(xDens)
218	#endif
219	#ifdef YP_VAR
220	solnData(YP_VAR,i,j,k) = xXp
221	#endif
222	#ifdef YN_VAR
223	solnData(YN_VAR,i,j,k) = xXn
224	#endif
225	#ifdef YA_VAR
226	solnData(YA_VAR,i,j,k) = xXa
227	#endif
228	#ifdef YH_VAR
229	solnData(YH_VAR,i,j,k) = xXh
230	#endif
231
232	#else
233	! CIVL's update subroutine
234
235	#endif
236
237	enddo
238	enddo
239	enddo
240
241	#ifndef _CIVL
242	call Grid_releaseBlkPtr(blockID,solndata)
243	call Eos_wrapped(MODE_DENS_EI,blkLimits,blockID)
244	#else
245	call Eos_wrapped(MODE_DENS_EI,blkLimits,data(:,:,:,:,blockID))
246	#endif
247
248	deallocate(xCenter)
249	if(NDIM>1)deallocate(yCenter)
250	if(NDIM>2)deallocate(zCenter)
251
252	enddo
253	!$omp enddo
254	!$omp end parallel
255
256	#ifdef _CIVL
257	! memcpy data to data_out
258	#endif
259
260	return
261	end subroutine Heat
262
263
264	function tauNu(dens)
265
266	implicit none
267
268	real, intent(IN) :: dens
269	real :: tauNu
270
271	! Other optical depth approximations may be used here.
272	tauNu = dens * 1.0e-11
273
274	return
275	end function tauNu

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format