Skip to content

Commit 6a8c2c9

Browse files
illorenzo7Loren
illorenzo7
authored and
Loren
committed
settled finally on term break up for induction--I think I got it right
1 parent 31ac086 commit 6a8c2c9

File tree

2 files changed

+32
-59
lines changed

2 files changed

+32
-59
lines changed

‎src/Diagnostics/Diagnostics_Base.F90‎

Lines changed: 5 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@ Module Diagnostics_Base
8888
Integer, Parameter :: ialt_offset = custom_offset+100 !2300
8989
Integer, Parameter :: ialt_offset2 = ialt_offset+100 !2400
9090
Integer, Parameter :: adv_offset = 2*3
91-
Integer, Parameter :: comp_offset = adv_offset + 5*3
91+
Integer, Parameter :: comp_offset = adv_offset + 3*3
9292
Integer, Parameter :: set_offset4 = comp_offset + 2*3 ! how many terms in a set for ialt2
9393

9494
! actual induction (not multiplied by B)
@@ -225,6 +225,7 @@ Module Diagnostics_Base
225225
Integer, Parameter :: ialtshear_work_r1 = ialt_offset2 + 1
226226
Integer, Parameter :: ialtshear_work_t1 = ialt_offset2 + 2
227227
Integer, Parameter :: ialtshear_work_p1 = ialt_offset2 + 3
228+
228229
Integer, Parameter :: ialtshear_work_r2 = ialt_offset2 + 4
229230
Integer, Parameter :: ialtshear_work_t2 = ialt_offset2 + 5
230231
Integer, Parameter :: ialtshear_work_p2 = ialt_offset2 + 6
@@ -242,17 +243,11 @@ Module Diagnostics_Base
242243
Integer, Parameter :: ialtadvec_work_t3 = ialt_offset2 + 8 + adv_offset
243244
Integer, Parameter :: ialtadvec_work_p3 = ialt_offset2 + 9 + adv_offset
244245

245-
Integer, Parameter :: ialtadvec_work_r4 = ialt_offset2 + 10 + adv_offset
246-
Integer, Parameter :: ialtadvec_work_t4 = ialt_offset2 + 11 + adv_offset
247-
Integer, Parameter :: ialtadvec_work_p4 = ialt_offset2 + 12 + adv_offset
248-
249-
Integer, Parameter :: ialtadvec_work_r5 = ialt_offset2 + 13 + adv_offset
250-
Integer, Parameter :: ialtadvec_work_t5 = ialt_offset2 + 14 + adv_offset
251-
Integer, Parameter :: ialtadvec_work_p5 = ialt_offset2 + 15 + adv_offset
252246

253247
Integer, Parameter :: ialtcomp_work_r1 = ialt_offset2 + 1 + comp_offset
254248
Integer, Parameter :: ialtcomp_work_t1 = ialt_offset2 + 2 + comp_offset
255249
Integer, Parameter :: ialtcomp_work_p1 = ialt_offset2 + 3 + comp_offset
250+
256251
Integer, Parameter :: ialtcomp_work_r2 = ialt_offset2 + 4 + comp_offset
257252
Integer, Parameter :: ialtcomp_work_t2 = ialt_offset2 + 5 + comp_offset
258253
Integer, Parameter :: ialtcomp_work_p2 = ialt_offset2 + 6 + comp_offset
@@ -278,6 +273,7 @@ Module Diagnostics_Base
278273
Integer, Parameter :: ialtshear_r1 = ialt_offset4 + 1
279274
Integer, Parameter :: ialtshear_t1 = ialt_offset4 + 2
280275
Integer, Parameter :: ialtshear_p1 = ialt_offset4 + 3
276+
281277
Integer, Parameter :: ialtshear_r2 = ialt_offset4 + 4
282278
Integer, Parameter :: ialtshear_t2 = ialt_offset4 + 5
283279
Integer, Parameter :: ialtshear_p2 = ialt_offset4 + 6
@@ -295,17 +291,11 @@ Module Diagnostics_Base
295291
Integer, Parameter :: ialtadvec_t3 = ialt_offset4 + 8 + adv_offset
296292
Integer, Parameter :: ialtadvec_p3 = ialt_offset4 + 9 + adv_offset
297293

298-
Integer, Parameter :: ialtadvec_r4 = ialt_offset4 + 10 + adv_offset
299-
Integer, Parameter :: ialtadvec_t4 = ialt_offset4 + 11 + adv_offset
300-
Integer, Parameter :: ialtadvec_p4 = ialt_offset4 + 12 + adv_offset
301-
302-
Integer, Parameter :: ialtadvec_r5 = ialt_offset4 + 13 + adv_offset
303-
Integer, Parameter :: ialtadvec_t5 = ialt_offset4 + 14 + adv_offset
304-
Integer, Parameter :: ialtadvec_p5 = ialt_offset4 + 15 + adv_offset
305294

306295
Integer, Parameter :: ialtcomp_r1 = ialt_offset4 + 1 + comp_offset
307296
Integer, Parameter :: ialtcomp_t1 = ialt_offset4 + 2 + comp_offset
308297
Integer, Parameter :: ialtcomp_p1 = ialt_offset4 + 3 + comp_offset
298+
309299
Integer, Parameter :: ialtcomp_r2 = ialt_offset4 + 4 + comp_offset
310300
Integer, Parameter :: ialtcomp_t2 = ialt_offset4 + 5 + comp_offset
311301
Integer, Parameter :: ialtcomp_p2 = ialt_offset4 + 6 + comp_offset

‎src/Diagnostics/Diagnostics_Custom.F90‎

Lines changed: 27 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -114,9 +114,9 @@ Subroutine Custom_MHD_Diagnostics(buffer)
114114
Real*8, Allocatable :: buff1(:,:,:,:),buff2(:,:,:,:)
115115
Real*8, Allocatable :: inducttmp(:,:,:,:), sheartmp(:,:,:,:), advtmp(:,:,:,:), comptmp(:,:,:,:), bfieldtmp(:,:,:,:)
116116
Real*8, Allocatable :: sheartmp1(:,:,:,:), sheartmp2(:,:,:,:)
117-
Real*8, Allocatable :: advtmp1(:,:,:,:), advtmp2(:,:,:,:), advtmp3(:,:,:,:), advtmp4(:,:,:,:), advtmp5(:,:,:,:)
117+
Real*8, Allocatable :: advtmp1(:,:,:,:), advtmp2(:,:,:,:), advtmp3(:,:,:,:)
118118
Real*8, Allocatable :: comptmp1(:,:,:,:), comptmp2(:,:,:,:)
119-
119+
Real*8, Allocatable :: divv_r(:,:,:), divv_theta(:,:,:), divv_phi(:,:,:)
120120

121121
Real*8 :: del2b
122122
Real*8, Allocatable :: ovstheta(:), ovs2theta(:)
@@ -146,12 +146,14 @@ Subroutine Custom_MHD_Diagnostics(buffer)
146146
Allocate(advtmp1(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max,1:3))
147147
Allocate(advtmp2(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max,1:3))
148148
Allocate(advtmp3(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max,1:3))
149-
Allocate(advtmp4(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max,1:3))
150-
Allocate(advtmp5(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max,1:3))
151149

152150
Allocate(comptmp1(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max,1:3))
153151
Allocate(comptmp2(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max,1:3))
154152

153+
Allocate(divv_r(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max))
154+
Allocate(divv_theta(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max))
155+
Allocate(divv_phi(1:n_phi,my_r%min:my_r%max,my_theta%min:my_theta%max))
156+
155157
ovstheta = 1.0d0/sintheta
156158
ovs2theta = 1.0d0/sin2theta
157159

@@ -1110,36 +1112,39 @@ Subroutine Custom_MHD_Diagnostics(buffer)
11101112
advtmp1(PSI,1) = -buff1(PSI,vr)*buff2(PSI,dbrdr)
11111113
advtmp2(PSI,1) = -buff1(PSI,vtheta)*buff2(PSI,dbrdt)/radius(r)
11121114
advtmp3(PSI,1) = -buff1(PSI,vphi)*buff2(PSI,dbrdp)/radius(r)/sintheta(t)
1113-
advtmp4(PSI,1) = -2.0d0*buff1(PSI,vr)*buff2(PSI,br)/radius(r)
1114-
advtmp5(PSI,1) = 0.0d0
11151115

11161116
advtmp1(PSI,2) = -buff1(PSI,vr)*buff2(PSI,dbtdr)
11171117
advtmp2(PSI,2) = -buff1(PSI,vtheta)*buff2(PSI,dbtdt)/radius(r)
11181118
advtmp3(PSI,2) = -buff1(PSI,vphi)*buff2(PSI,dbtdp)/radius(r)/sintheta(t)
1119-
advtmp4(PSI,2) = buff1(PSI,vr)*buff2(PSI,btheta)/radius(r)
1120-
advtmp5(PSI,2) = -buff1(PSI,vtheta)*buff2(PSI,btheta)*cottheta(t)/radius(r)
11211119

11221120
advtmp1(PSI,3) = -buff1(PSI,vr)*buff2(PSI,dbpdr)
11231121
advtmp2(PSI,3) = -buff1(PSI,vtheta)*buff2(PSI,dbpdt)/radius(r)
11241122
advtmp3(PSI,3) = -buff1(PSI,vphi)*buff2(PSI,dbpdp)/radius(r)/sintheta(t)
1125-
advtmp4(PSI,3) = buff1(PSI,vr)*buff2(PSI,bphi)/radius(r)
1126-
advtmp5(PSI,3) = buff1(PSI,vtheta)*buff2(PSI,bphi)*cottheta(t)/radius(r)
11271123

11281124
Do j = 1,3
1129-
advtmp(PSI,j) = advtmp1(PSI,j) + advtmp2(PSI,j) + advtmp3(PSI,j) +&
1130-
&advtmp4(PSI,j) + advtmp5(PSI,j)
1125+
advtmp(PSI,j) = advtmp1(PSI,j) + advtmp2(PSI,j) + advtmp3(PSI,j)
11311126
Enddo
11321127
END_DO
11331128

11341129
! get compression
11351130
! ALTERNATE compression: -B (div v) (but TRANSVERSE compression only)
1131+
! 3 parts: transverse to r, theta, phi (do this separately, I think it
1132+
! sums to the right thing
11361133
DO_PSI
1137-
comptmp1(PSI,1) = -buff2(PSI,br)*(buff1(PSI,dvtdt)+buff1(PSI,vtheta)*cottheta(t))/radius(r)
1138-
comptmp2(PSI,1) = -buff2(PSI,br)*buff1(PSI,dvpdp)/radius(r)/sintheta(t)
1139-
comptmp1(PSI,2) = -buff2(PSI,btheta)*(buff1(PSI,dvrdr)+2.0d0*buff1(PSI,vr)/radius(r))
1140-
comptmp2(PSI,2) = -buff2(PSI,btheta)*buff1(PSI,dvpdp)/radius(r)/sintheta(t)
1141-
comptmp1(PSI,3) = -buff2(PSI,bphi)*(buff1(PSI,dvrdr)+2.0d0*buff1(PSI,vr)/radius(r))
1142-
comptmp2(PSI,3) = -buff2(PSI,bphi)*(buff1(PSI,dvtdt)+buff1(PSI,vtheta)*cottheta(t))/radius(r)
1134+
divv_r(PSI) = -buff1(PSI, dvrdr)
1135+
divv_theta(PSI) = -(buff1(PSI,dvtdt)+buff1(PSI,vr))/radius(r)
1136+
divv_phi(PSI) = -(buff1(PSI,dvpdp)/sintheta(t) +&
1137+
&buff1(PSI,vr) + cottheta(t)*buff1(PSI,vtheta))/radius(r)
1138+
1139+
comptmp1(PSI,1) = buff2(PSI,br)*divv_theta(PSI)
1140+
comptmp2(PSI,1) = buff2(PSI,br)*divv_phi(PSI)
1141+
1142+
comptmp1(PSI,2) = buff2(PSI,btheta)*divv_r(PSI)
1143+
comptmp2(PSI,2) = buff2(PSI,btheta)*divv_phi(PSI)
1144+
1145+
comptmp1(PSI,3) = buff2(PSI,bphi)*divv_r(PSI)
1146+
comptmp2(PSI,3) = buff2(PSI,bphi)*divv_theta(PSI)
1147+
11431148
Do j = 1,3
11441149
comptmp(PSI,j) = comptmp1(PSI,j) + comptmp2(PSI,j)
11451150
Enddo
@@ -1226,18 +1231,6 @@ Subroutine Custom_MHD_Diagnostics(buffer)
12261231
END_DO
12271232
Call Add_Quantity(qty)
12281233
Endif
1229-
If (compute_quantity(ialtadvec_work_r4+j-1+offset2)) Then
1230-
DO_PSI
1231-
qty(PSI) = advtmp4(PSI,j)*bfieldtmp(PSI,j)
1232-
END_DO
1233-
Call Add_Quantity(qty)
1234-
Endif
1235-
If (compute_quantity(ialtadvec_work_r5+j-1+offset2)) Then
1236-
DO_PSI
1237-
qty(PSI) = advtmp5(PSI,j)*bfieldtmp(PSI,j)
1238-
END_DO
1239-
Call Add_Quantity(qty)
1240-
Endif
12411234

12421235
! Comp. terms
12431236
If (compute_quantity(ialtcomp_work_r1+j-1+offset2)) Then
@@ -1325,18 +1318,6 @@ Subroutine Custom_MHD_Diagnostics(buffer)
13251318
END_DO
13261319
Call Add_Quantity(qty)
13271320
Endif
1328-
If (compute_quantity(ialtadvec_r4+j-1+offset2)) Then
1329-
DO_PSI
1330-
qty(PSI) = advtmp4(PSI,j)
1331-
END_DO
1332-
Call Add_Quantity(qty)
1333-
Endif
1334-
If (compute_quantity(ialtadvec_r5+j-1+offset2)) Then
1335-
DO_PSI
1336-
qty(PSI) = advtmp5(PSI,j)
1337-
END_DO
1338-
Call Add_Quantity(qty)
1339-
Endif
13401321

13411322
! Comp. terms
13421323
If (compute_quantity(ialtcomp_r1+j-1+offset2)) Then
@@ -1376,12 +1357,14 @@ Subroutine Custom_MHD_Diagnostics(buffer)
13761357
DeAllocate(advtmp1)
13771358
DeAllocate(advtmp2)
13781359
DeAllocate(advtmp3)
1379-
DeAllocate(advtmp4)
1380-
DeAllocate(advtmp5)
13811360

13821361
DeAllocate(comptmp1)
13831362
DeAllocate(comptmp2)
13841363

1364+
DeAllocate(divv_r)
1365+
DeAllocate(divv_theta)
1366+
DeAllocate(divv_phi)
1367+
13851368
! Edit Above This Line
13861369
!=============================================
13871370

0 commit comments

Comments
 (0)