-
Notifications
You must be signed in to change notification settings - Fork 26
/
potential_root.f90
631 lines (523 loc) · 26.5 KB
/
potential_root.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
submodule (potential_comm) potential_root
implicit none (type, external)
contains
module procedure potential_root_mpi_curv
!! ROOT MPI COMM./SOLVE ROUTINE FOR POTENTIAL. THIS VERSION
!! INCLUDES THE POLARIZATION CURRENT TIME DERIVATIVE PART
!! AND CONVECTIVE PARTS IN MATRIX SOLUTION.
!! STATE VARIABLES VS2,3 INCLUDE GHOST CELLS. FOR NOW THE
!! POLARIZATION TERMS ARE PASSED BACK TO MAIN FN, EVEN THOUGH
!! THEY ARE NOT USED (THEY MAY BE IN THE FUTURE)
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: v2,v3
real(wp), dimension(1:size(Phiall,1),1:size(Phiall,2),1:size(Phiall,3)) :: srctermall
real(wp), dimension(1:size(Phiall,2),1:size(Phiall,3)), target :: Vminx1,Vmaxx1 !allow pointer aliases for these vars.
real(wp), dimension(1:size(Phiall,2),1:size(Phiall,3)) :: Vminx1buf,Vmaxx1buf
real(wp), dimension(1:size(Phiall,1),1:size(Phiall,3)) :: Vminx2,Vmaxx2
real(wp), dimension(1:size(Phiall,1),1:size(Phiall,2)) :: Vminx3,Vmaxx3
integer :: flagdirich
real(wp), dimension(1:size(Phiall,2),1:size(Phiall,3)) :: v2slaball,v3slaball !stores drift velocs. for pol. current
! real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: paramtrim !to hold trimmed magnetic field
real(wp), dimension(1:size(Phiall,1),1:size(Phiall,2),1:size(Phiall,3)) :: E01all,E02all,E03all !background fields
! real(wp), dimension(1:size(Phiall,1),1:size(Phiall,2),1:size(Phiall,3)) :: divJperpall2,divJperpall3,divJperpall
!! more work arrays
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: integrand,sigintegral !general work array for doing integrals
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: grad2E,grad3E !more work arrays for pol. curr.
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: DE2Dt,DE3Dt !pol. drift
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: J1pol,J2pol,J3pol
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: E01,E02,E03 !distributed background fields
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: srcterm,divJperp
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: E1prev,E2prev,E3prev
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: Phi
real(wp), dimension(1:size(E1,2),1:size(E1,3)) :: SigPint2,SigPint3,SigHint,incapint,srctermint
real(wp), dimension(1:size(Phiall,2),1:size(Phiall,3)) :: SigPint2all,SigPint3all,SigHintall,incapintall,srctermintall
real(wp), dimension(1:size(Phiall,2),1:size(Phiall,3)) :: Phislab,Phislab0
real(wp), dimension(0:size(E1,1)+1,0:size(E1,2)+1,0:size(E1,3)+1) :: divtmp
!! one extra grid point on either end to facilitate derivatives
real(wp), dimension(-1:size(E1,1)+2,-1:size(E1,2)+2,-1:size(E1,3)+2) :: J1halo,J2halo,J3halo
!! haloing assumes existence of two ghost cells
real(wp), dimension(1:size(Phiall,1),1:size(Phiall,2),1:size(Phiall,3)) :: sig0scaledall,sigPscaledall,sigHscaledall
real(wp), dimension(1:size(E1,1),1:size(E1,2),1:size(E1,3)) :: sig0scaled,sigPscaled,sigHscaled
logical :: perflag !MUMPS stuff
real(wp), dimension(1:size(Phiall,3)) :: Vminx2slice,Vmaxx2slice
real(wp), dimension(1:size(Phiall,2)) :: Vminx3slice,Vmaxx3slice
real(wp), dimension(1:size(E1,2),1:size(E1,3)) :: Vminx1slab,Vmaxx1slab
real(wp), dimension(1:size(E1,2),1:size(E1,3)) :: v2slab,v3slab
integer :: iid, ierr
integer :: ix1,ix2,ix3,lx1,lx2,lx3,lx3all
integer :: idleft,idright,iddown,idup
real(wp) :: tstart,tfin
!SIZES - PERHAPS SHOULD BE TAKEN FROM GRID MODULE INSTEAD OF RECOMPUTED?
lx1=size(sig0,1)
lx2=size(sig0,2)
lx3=size(sig0,3)
lx3all=size(Phiall,3)
!USE PREVIOUS MUMPS PERMUTATION (OLD CODE? BUT MIGHT BE WORTH REINSTATING?)
! perflag=.false.
perflag=.true.
!R-------
!! POPULATE BACKGROUND AND BOUNDARY CONDITION ARRAYS
!! - IDEALLY ROOT ONLY SINCE IT INVOLVES FILE INPUT, ALTHOUGH THE INTERPOLATION MAY BE SLOW...
call cpu_time(tstart)
if (flagE0file==1) then
call potentialBCs2D_fileinput(dt,dtE0,t,ymd,UTsec,E0dir,x,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3,Vmaxx3, &
E01all,E02all,E03all,flagdirich)
else
call potentialBCs2D(t,x,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3,Vmaxx3, &
E01all,E02all,E03all,flagdirich) !user needs to manually swap x2 and x3 in this function.
end if
call cpu_time(tfin)
if (debug) print *, 'Root has computed BCs in time: ',tfin-tstart
!R-------
!R--------
ierr=0
do iid=1,lid-1 !communicate intent for solve to workers so they know whether or not to call mumps fn.
call mpi_send(flagdirich,1,MPI_INTEGER,iid,tag%flagdirich,MPI_COMM_WORLD,ierr)
end do
if (ierr/=0) error stop 'mpi_send failed to send solve intent'
if (debug) print *, 'Root has communicated type of solve to workers: ',flagdirich
!R--------
!Need to broadcast background fields from root
!Need to also broadcast x1 boundary conditions for source term calculations.
call bcast_send(E01all,tag%E01,E01)
call bcast_send(E02all,tag%E02,E02)
call bcast_send(E03all,tag%E03,E03)
!These are pointer targets so don't assume contiguous in memory - pack them into a buffer to be safe
Vmaxx1buf=Vmaxx1; Vminx1buf=Vminx1;
call bcast_send(Vminx1buf,tag%Vminx1,Vminx1slab)
call bcast_send(Vmaxx1buf,tag%Vmaxx1,Vmaxx1slab)
!-------
!CONDUCTION CURRENT BACKGROUND SOURCE TERMS FOR POTENTIAL EQUATION. MUST COME AFTER CALL TO BC CODE.
J1=0d0 !so this div is only perp components
if (flagswap==1) then
J2=sigP*E02+sigH*E03 !BG x2 current
J3=-1*sigH*E02+sigP*E03 !BG x3 current
else
J2=sigP*E02-sigH*E03 !BG x2 current
J3=sigH*E02+sigP*E03 !BG x3 current
end if
J1halo(1:lx1,1:lx2,1:lx3)=J1
J2halo(1:lx1,1:lx2,1:lx3)=J2
J3halo(1:lx1,1:lx2,1:lx3)=J3
call halo_pot(J1halo,tag%J1,x%flagper,.false.)
call halo_pot(J2halo,tag%J2,x%flagper,.false.)
call halo_pot(J3halo,tag%J3,x%flagper,.false.)
divtmp=div3D(J1halo(0:lx1+1,0:lx2+1,0:lx3+1),J2halo(0:lx1+1,0:lx2+1,0:lx3+1), &
J3halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
srcterm=divtmp(1:lx1,1:lx2,1:lx3)
if (debug) print *, 'Root has computed background field source terms...',minval(srcterm), maxval(srcterm)
!print*, myid, any(ieee_is_nan(J1halo(0:lx1+1,1:lx2,1:lx3))), &
! any(ieee_is_nan(J2halo(1:lx1,0:lx2+1,1:lx3))), &
! any(ieee_is_nan(J3halo(1:lx1,1:lx2,0:lx3+1))), &
! any(ieee_is_nan(divtmp(1:lx1,1:lx2,1:lx3)))
!-------
!-------
!NEUTRAL WIND SOURCE TERMS FOR POTENTIAL EQUATION, SIMILAR TO ABOVE BLOCK OF CODE
J1=0d0 !so this div is only perp components
if (flagswap==1) then
J2=-1*sigP*vn3*B1(1:lx1,1:lx2,1:lx3)+sigH*vn2*B1(1:lx1,1:lx2,1:lx3)
!! wind x2 current, note that all workers already have a copy of this.
J3=sigH*vn3*B1(1:lx1,1:lx2,1:lx3)+sigP*vn2*B1(1:lx1,1:lx2,1:lx3)
!! wind x3 current
else
J2=sigP*vn3*B1(1:lx1,1:lx2,1:lx3)+sigH*vn2*B1(1:lx1,1:lx2,1:lx3)
!! wind x2 current
J3=sigH*vn3*B1(1:lx1,1:lx2,1:lx3)-sigP*vn2*B1(1:lx1,1:lx2,1:lx3)
!! wind x3 current
end if
J1halo(1:lx1,1:lx2,1:lx3)=J1
J2halo(1:lx1,1:lx2,1:lx3)=J2
J3halo(1:lx1,1:lx2,1:lx3)=J3
call halo_pot(J1halo,tag%J1,x%flagper,.false.)
call halo_pot(J2halo,tag%J2,x%flagper,.false.)
call halo_pot(J3halo,tag%J3,x%flagper,.false.)
divtmp=div3D(J1halo(0:lx1+1,0:lx2+1,0:lx3+1),J2halo(0:lx1+1,0:lx2+1,0:lx3+1), &
J3halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
srcterm=srcterm+divtmp(1:lx1,1:lx2,1:lx3)
if (debug) print *, 'Root has computed wind source terms...',minval(srcterm), maxval(srcterm)
!-------
!!!!!!!!
!-----AT THIS POINT WE MUST DECIDE WHETHER TO DO AN INTEGRATED SOLVE OR A 2D FIELD-RESOLVED SOLVED
!-----DECIDE BASED ON THE SIZE OF THE X2 DIMENSION
if (lx2/=1) then !either field-resolved 3D or integrated 2D solve for 3D domain
if (potsolve == 1) then !2D, field-integrated solve
if (debug) print *, 'Beginning field-integrated solve...'
!> INTEGRATE CONDUCTANCES AND CAPACITANCES FOR SOLVER COEFFICIENTS
integrand=sigP*x%h1(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)/x%h2(1:lx1,1:lx2,1:lx3)
sigintegral=integral3D1(integrand,x,1,lx1) !no haloing required for a field-line integration
SigPint2=sigintegral(lx1,:,:)
integrand=sigP*x%h1(1:lx1,1:lx2,1:lx3)*x%h2(1:lx1,1:lx2,1:lx3)/x%h3(1:lx1,1:lx2,1:lx3)
sigintegral=integral3D1(integrand,x,1,lx1)
SigPint3=sigintegral(lx1,:,:)
integrand=x%h1(1:lx1,1:lx2,1:lx3)*sigH
sigintegral=integral3D1(integrand,x,1,lx1)
SigHint=sigintegral(lx1,:,:)
sigintegral=integral3D1(incap,x,1,lx1)
incapint=sigintegral(lx1,:,:)
!-------
!> PRODUCE A FIELD-INTEGRATED SOURCE TERM
if (flagdirich /= 1) then !Neumann conditions; incorporate a source term and execute the solve
if (debug) print *, 'Using FAC boundary condition...'
!-------
integrand=x%h1(1:lx1,1:lx2,1:lx3)*x%h2(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)*srcterm
sigintegral=integral3D1(integrand,x,1,lx1)
srctermint=sigintegral(lx1,:,:)
srctermint=srctermint+x%h2(lx1,1:lx2,1:lx3)*x%h3(lx1,1:lx2,1:lx3)*Vmaxx1slab- &
x%h2(1,1:lx2,1:lx3)*x%h3(1,1:lx2,1:lx3)*Vminx1slab
!! workers don't have access to boundary conditions, unless root sends
!-------
v2=vs2(1:lx1,1:lx2,1:lx3,1); v3=vs3(1:lx1,1:lx2,1:lx3,1);
! must be set since used later by the polarization current calculation
v2slab=v2(lx1,1:lx2,1:lx3); v3slab=v3(lx1,1:lx2,1:lx3);
!! need to pick out the ExB drift here (i.e. the drifts from highest altitudes); but this is only valid for Cartesian,
!! so it's okay for the foreseeable future
!RADD--- ROOT NEEDS TO PICK UP *INTEGRATED* SOURCE TERMS AND COEFFICIENTS FROM WORKERS
call gather_recv(srctermint,tag%src,srctermintall)
call gather_recv(incapint,tag%incapint,incapintall)
call gather_recv(SigPint2,tag%SigPint2,SigPint2all)
call gather_recv(SigPint3,tag%SigPint3,SigPint3all)
call gather_recv(SigHint,tag%SigHint,SigHintall)
call gather_recv(v2slab,tag%v2electro,v2slaball)
call gather_recv(v3slab,tag%v3electro,v3slaball)
! if (t>480) then
! open(newunit=utrace, form='unformatted', access='stream', file='scrtermintall.raw8', status='replace', action='write')
! write(utrace) srctermintall
! close(utrace)
! error stop 'DEBUG'
! end if
!R------
!EXECUTE FIELD-INTEGRATED SOLVE
Vminx2slice=Vminx2(lx1,:) !slice the boundaries into expected shape
Vmaxx2slice=Vmaxx2(lx1,:)
Vminx3slice=Vminx3(lx1,:)
Vmaxx3slice=Vmaxx3(lx1,:)
Phislab0=Phiall(lx1,:,:) !root already possess the fullgrid potential from prior solves...
if (debug) print *, 'Root is calling MUMPS...'
!R-------
!R------ EXECUTE THE MUMPS SOLVE FOR FIELD-INT
call cpu_time(tstart)
if (.not. x%flagper) then !nonperiodic mesh
if (debug) print *, '!!!User selected aperiodic solve...'
Phislab=potential2D_polarization(srctermintall,SigPint2all,SigPint3all,SigHintall,incapintall,v2slaball,v3slaball, &
Vminx2slice,Vmaxx2slice,Vminx3slice,Vmaxx3slice, &
dt,x,Phislab0,perflag,it)
!! note tha this solver is only valid for cartesian meshes, unless the inertial capacitance is set to zero
else
if (debug) print *, '!!!User selected periodic solve...'
Phislab = potential2D_polarization_periodic(srctermintall,SigPint2all,SigHintall,incapintall,v2slaball,v3slaball, &
Vminx2slice,Vmaxx2slice,Vminx3slice,Vmaxx3slice, &
dt,x,Phislab0,perflag,it)
!! !note that either sigPint2 or 3 will work since this must be cartesian...
end if
call cpu_time(tfin)
if (debug) print *, 'Root received results from MUMPS which took time: ',tfin-tstart
!R-------
else
!! Dirichlet conditions - since this is field integrated we just copy BCs specified by user
!! to other locations along field line
!R------
Phislab=Vmaxx1
!! potential is whatever user specifies, since we assume equipotential field lines,
!! it doesn't really matter whether we use Vmaxx1 or Vminx1.
!! Note however, tha thte boundary conditions subroutines should explicitly
!! set these to be equal with Dirichlet conditions, for consistency.
if (debug) print *, 'Dirichlet conditions selected with field-integrated solve. Copying BCs along x1-direction...'
!R------
end if
!R------ AFTER ANY TYPE OF FIELD-INT SOLVE COPY THE BCS ACROSS X1 DIMENSION
do ix1=1,lx1
Phiall(ix1,:,:)=Phislab(:,:)
!! copy the potential across the ix1 direction; past this point there is no real difference with 3D,
!! note that this is still valid in curvilinear form
end do
!R------
else
!! resolved 3D solve
!! ZZZ - conductivities need to be properly scaled here...
!! So does the source term... Maybe leave as broken for now since I don't really plan to use this code
if (debug) print *, 'Beginning field-resolved 3D solve... Type; ',flagdirich
!-------
!PRODUCE SCALED CONDUCTIVITIES TO PASS TO SOLVER, ALSO SCALED SOURCE TERM,
!need to adopt for curvilinear case...
sig0scaled=x%h2(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)/x%h1(1:lx1,1:lx2,1:lx3)*sig0
if (flagswap==1) then
sigPscaled=x%h1(1:lx1,1:lx2,1:lx3)*x%h2(1:lx1,1:lx2,1:lx3)/x%h3(1:lx1,1:lx2,1:lx3)*sigP !remember to swap 2-->3
else
sigPscaled=x%h1(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)/x%h2(1:lx1,1:lx2,1:lx3)*sigP
end if
srcterm=srcterm*x%h1(1:lx1,1:lx2,1:lx3)*x%h2(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)
sigHscaled=x%h1(1:lx1,1:lx2,1:lx3)*sigH
!RADD--- ROOT NEEDS TO PICK UP FIELD-RESOLVED SOURCE TERM AND COEFFICIENTS FROM WORKERS
call gather_recv(sigPscaled,tag%sigP,sigPscaledall)
call gather_recv(sigHscaled,tag%sigH,sigHscaledall)
call gather_recv(sig0scaled,tag%sig0,sig0scaledall)
call gather_recv(srcterm,tag%src,srctermall)
!R------
if (debug) print *, '!Beginning field-resolved 3D solve (could take a very long time)...'
! Phiall=elliptic3D_curv(srctermall,sig0scaledall,sigPscaledall,sigHscaledall,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3,Vmaxx3, &
! x,flagdirich,perflag,it)
if( maxval(abs(Vminx1))>1e-12_wp .or. maxval(abs(Vmaxx1))>1e-12_wp ) then
do iid=1,lid-1
call mpi_send(1,1,MPI_INTEGER,iid,tag%flagdirich,MPI_COMM_WORLD,ierr)
end do
Phiall=potential3D_fieldresolved_decimate(srctermall,sig0scaledall,sigPscaledall,sigHscaledall, &
Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3,Vmaxx3, &
x,flagdirich,perflag,it)
else
do iid=1,lid-1
call mpi_send(0,1,MPI_INTEGER,iid,tag%flagdirich,MPI_COMM_WORLD,ierr)
end do
if (debug) print*, 'Boundary conditions too small to require solve, setting everything to zero...'
Phiall=0e0_wp
end if
!R------
end if
else !lx1=1 so do a field-resolved 2D solve over x1,x3
if (debug) print *, 'Beginning field-resolved 2D solve... Type; ',flagdirich
!-------
!PRODUCE SCALED CONDUCTIVITIES TO PASS TO SOLVER, ALSO SCALED SOURCE TERM
sig0scaled=x%h2(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)/x%h1(1:lx1,1:lx2,1:lx3)*sig0
if (flagswap==1) then
sigPscaled=x%h1(1:lx1,1:lx2,1:lx3)*x%h2(1:lx1,1:lx2,1:lx3)/x%h3(1:lx1,1:lx2,1:lx3)*sigP !remember to swap 2-->3
else
sigPscaled=x%h1(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)/x%h2(1:lx1,1:lx2,1:lx3)*sigP
end if
srcterm=srcterm*x%h1(1:lx1,1:lx2,1:lx3)*x%h2(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)
! srcterm=-1d0*srcterm
!! in a 2D solve negate this due to it being a cross produce and the fact that we've permuted the 2 and 3 dimensions.
!! ZZZ - NOT JUST THIS WORKS WITH BACKGROUND FIELDS???
!-------
!RADD--- NEED TO GET THE RESOLVED SOURCE TERMS AND COEFFICIENTS FROM WORKERS
call gather_recv(sigPscaled,tag%sigP,sigPscaledall)
call gather_recv(sig0scaled,tag%sig0,sig0scaledall)
call gather_recv(srcterm,tag%src,srctermall)
!! EXECUTE THE SOLVE WITH MUMPS AND SCALED TERMS
!! NOTE THE LACK OF A SPECIAL CASE HERE TO CHANGE THE POTENTIAL PROBLEM
!! - ONLY THE HALL TERM CHANGES (SINCE RELATED TO EXB) BUT THAT DOESN'T APPEAR IN THIS EQN!
Phiall=potential2D_fieldresolved(srctermall,sig0scaledall,sigPscaledall,Vminx1,Vmaxx1,Vminx3,Vmaxx3, &
x,flagdirich,perflag,it)
end if
if (debug) print *, 'MUMPS time: ',tfin-tstart
!!!!!!!!!
!RADD--- ROOT NEEDS TO PUSH THE POTENTIAL BACK TO ALL WORKERS FOR FURTHER PROCESSING (BELOW)
call bcast_send(Phiall,tag%Phi,Phi)
!-------
!! STORE PREVIOUS TIME TOTAL FIELDS BEFORE UPDATING THE ELECTRIC FIELDS WITH NEW POTENTIAL
!! (OLD FIELDS USED TO CALCULATE POLARIZATION CURRENT)
E1prev=E1
E2prev=E2
E3prev=E3
!-------
!-------
!CALCULATE PERP FIELDS FROM POTENTIAL
! E20all=grad3D2(-1d0*Phi0all,dx2(1:lx2))
!! causes major memory leak. maybe from arithmetic statement argument?
!! Left here as a 'lesson learned' (or is it a gfortran bug...)
! E30all=grad3D3(-1d0*Phi0all,dx3all(1:lx3all))
Phi=-1d0*Phi
!COMPUTE THE 2 COMPONENT OF THE ELECTRIC FIELD
J1halo(1:lx1,1:lx2,1:lx3)=Phi
call halo_pot(J1halo,tag%J1,x%flagper,.true.)
divtmp=grad3D2(J1halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
E2=divtmp(1:lx1,1:lx2,1:lx3)
!COMPUTE THE 3 COMPONENT OF THE ELECTRIC FIELD
J1halo(1:lx1,1:lx2,1:lx3)=Phi
call halo_pot(J1halo,tag%J1,x%flagper,.false.)
divtmp=grad3D3(J1halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
E3=divtmp(1:lx1,1:lx2,1:lx3)
Phi=-1d0*Phi !put things back for later use
!--------
!R-------
!JUST TO JUDGE THE IMPACT OF MI COUPLING
if (debug) then
print *, 'Max integrated inertial capacitance: ',maxval(incapintall)
!print *, 'Max integrated Pedersen conductance (includes metric factors): ',maxval(SigPint2all)
!print *, 'Max integrated Hall conductance (includes metric factors): ',minval(SigHintall), maxval(SigHintall)
print *, 'Max E2,3 BG and response values are: ',maxval(E02), maxval(E03),maxval(E2),maxval(E3)
print *, 'Min E2,3 BG and response values are: ',minval(E02), minval(E03),minval(E2),minval(E3)
print *, 'Min/Max values of potential: ',minval(Phi),maxval(Phi)
print *, 'Min/Max values of full grid potential: ',minval(Phiall),maxval(Phiall)
endif
!R-------
!--------
!ADD IN BACKGROUND FIELDS BEFORE DISTRIBUTING TO WORKERS
E2=E2+E02
E3=E3+E03
!--------
!--------
!COMPUTE TIME DERIVATIVE NEEDED FOR POLARIZATION CURRENT. ONLY DO THIS IF WE HAVE SPECIFIC NONZERO INERTIAL CAPACITANCE
!if (maxval(incap) > 0._wp) then
!! ZZZ this is really bad needs to be a global test rather than having each worker test since there
!! is message passing embedded in here and everyone needs to do the same thing!!!
if (flagcap/=0) then
if (debug) print*, 'Working on polarization currents...'
!differentiate E2 in x2 (needs haloing)
J1halo(1:lx1,1:lx2,1:lx3)=E2
call halo_pot(J1halo,tag%J1,x%flagper,.false.)
divtmp=grad3D2(J1halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
grad2E=divtmp(1:lx1,1:lx2,1:lx3)
!Now differentiate E2 in the x3 direction
J1halo(1:lx1,1:lx2,1:lx3)=E2
call halo_pot(J1halo,tag%J1,x%flagper,.false.)
divtmp=grad3D3(J1halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
grad3E=divtmp(1:lx1,1:lx2,1:lx3)
!total derivative needed to compute the x2 component of the polarization current
DE2Dt=(E2-E2prev)/dt+v2*grad2E+v3*grad3E
!differentiate E3 in x2
J1halo(1:lx1,1:lx2,1:lx3)=E3
call halo_pot(J1halo,tag%J1,x%flagper,.false.)
divtmp=grad3D2(J1halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
grad2E=divtmp(1:lx1,1:lx2,1:lx3)
!differentiate E3 in x3
J1halo(1:lx1,1:lx2,1:lx3)=E3
call halo_pot(J1halo,tag%J1,x%flagper,.false.)
divtmp=grad3D3(J1halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
grad3E=divtmp(1:lx1,1:lx2,1:lx3)
!total derivative needed for x3 component of pol. current
DE3Dt=(E3-E3prev)/dt+v2*grad2E+v3*grad3E
!convert derivative to current density
J1pol=0d0
J2pol=incap*DE2Dt
J3pol=incap*DE3Dt
else !pure electrostatic solve was done
if (debug) print*, 'Electrostatics used, skipping polarization current...'
DE2Dt=0d0
DE3Dt=0d0
J1pol=0d0
J2pol=0d0
J3pol=0d0
end if
!--------
!--------
if (debug) print *, 'Working on perp. conduction currents...'
if (flagswap==1) then
J2=sigP*E2+sigH*E3 !BG field already added to E above
J3=-1*sigH*E2+sigP*E3
else
J2=sigP*E2-sigH*E3 !BG field already added to E above
J3=sigH*E2+sigP*E3
end if
!WHAT I THINK THE NEUTRAL WIND CURRENTS SHOULD BE IN 2D
if (flagswap==1) then
J2=J2-sigP*vn3*B1(1:lx1,1:lx2,1:lx3)+sigH*vn2*B1(1:lx1,1:lx2,1:lx3)
J3=J3+sigH*vn3*B1(1:lx1,1:lx2,1:lx3)+sigP*vn2*B1(1:lx1,1:lx2,1:lx3)
else
J2=J2+sigP*vn3*B1(1:lx1,1:lx2,1:lx3)+sigH*vn2*B1(1:lx1,1:lx2,1:lx3)
J3=J3+sigH*vn3*B1(1:lx1,1:lx2,1:lx3)-sigP*vn2*B1(1:lx1,1:lx2,1:lx3)
end if
!-------
!!!!!!!!
!NOW DEAL WITH THE PARALLEL FIELDS AND ALL CURRENTS
if (lx2/=1 .and. potsolve ==1) then !we did a field-integrated solve above
if (debug) print*, 'Appear to need to differentiate to get J1...'
!-------
!NOTE THAT A DIRECT E1ALL CALCULATION WILL GIVE ZERO, SO USE INDIRECT METHOD, AS FOLLOWS
J1=0d0 !a placeholder so that only the perp divergence is calculated - will get overwritten later.
! divJperp=div3D(J1,J2,J3,x,1,lx1,1,lx2,1,lx3)
J1halo(1:lx1,1:lx2,1:lx3)=J1
J2halo(1:lx1,1:lx2,1:lx3)=J2
J3halo(1:lx1,1:lx2,1:lx3)=J3
call halo_pot(J1halo,tag%J1,x%flagper,.false.)
call halo_pot(J2halo,tag%J2,x%flagper,.false.)
call halo_pot(J3halo,tag%J3,x%flagper,.false.)
divtmp=div3D(J1halo(0:lx1+1,0:lx2+1,0:lx3+1),J2halo(0:lx1+1,0:lx2+1,0:lx3+1), &
J3halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
divJperp=x%h1(1:lx1,1:lx2,1:lx3)*x%h2(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)*divtmp(1:lx1,1:lx2,1:lx3)
if (flagdirich /= 1) then
!! Neumann conditions, this is boundary location-agnostic since both bottom and top FACs are known
!! - they have to be loaded into VVmaxx1 and Vminx1
if (debug) print *, 'Nuemann boundaries, integrated from highest altitude down to preserve accuracy...'
if (gridflag==0) then !closed dipole grid, really would be best off integrating from the source hemisphere
if (debug) print *, 'Closed dipole grid; integration starting in source hemisphere (if applicable)...', &
minval(Vmaxx1slab), &
maxval(Vmaxx1slab)
if (sourcemlat>=0d0) then !integrate from northern hemisphere
if (debug) print *, 'Source is in northern hemisphere (or there is no source)...'
J1=integral3D1_curv_alt(divJperp,x,1,lx1) !int divperp of BG current, go from maxval(x1) to location of interest
do ix1=1,lx1
J1(ix1,:,:)=1d0/x%h2(ix1,1:lx2,1:lx3)/x%h3(ix1,1:lx2,1:lx3)* &
(x%h2(1,1:lx2,1:lx3)*x%h3(1,1:lx2,1:lx3)*Vmaxx1slab+J1(ix1,:,:))
end do
else
if (debug) print *, 'Source in southern hemisphere...'
J1=integral3D1(divJperp,x,1,lx1) !int divperp of BG current starting from minx1
do ix1=1,lx1
J1(ix1,:,:)=1d0/x%h2(ix1,1:lx2,1:lx3)/x%h3(ix1,1:lx2,1:lx3)* &
(x%h2(1,1:lx2,1:lx3)*x%h3(1,1:lx2,1:lx3)*Vminx1slab-J1(ix1,:,:))
end do
end if
elseif (gridflag==1) then !this would be an inverted grid, this max altitude corresponds to the min value of x1
if (debug) print *, 'Inverted grid; integration starting at min x1 (highest alt. or southern hemisphere)...', &
minval(Vminx1slab), &
maxval(Vminx1slab)
J1=integral3D1(divJperp,x,1,lx1) !int divperp of BG current
do ix1=1,lx1
J1(ix1,:,:)=1d0/x%h2(ix1,1:lx2,1:lx3)/x%h3(ix1,1:lx2,1:lx3)* &
(x%h2(1,1:lx2,1:lx3)*x%h3(1,1:lx2,1:lx3)*Vminx1slab-J1(ix1,:,:))
end do
else !minx1 is at teh bottom of the grid to integrate from max x1
if (debug) print *, 'Non-inverted grid; integration starting at max x1...', minval(Vmaxx1slab), maxval(Vmaxx1slab)
J1=integral3D1_curv_alt(divJperp,x,1,lx1) !int divperp of BG current, go from maxval(x1) to location of interest
do ix1=1,lx1
J1(ix1,:,:)=1d0/x%h2(ix1,1:lx2,1:lx3)/x%h3(ix1,1:lx2,1:lx3)* &
(x%h2(1,1:lx2,1:lx3)*x%h3(1,1:lx2,1:lx3)*Vmaxx1slab+J1(ix1,:,:))
end do
end if
else
!! Dirichlet conditions - we need to integrate from the ***lowest altitude***
!! (where FAC is known to be zero, note this is not necessarilty the logical bottom of the grid), upwards (to where it isn't)
if (gridflag/=2) then !inverted grid (logical top is the lowest altitude)
if (debug) print *, 'Inverted grid detected - integrating logical top downward to compute FAC...'
J1=integral3D1_curv_alt(divJperp,x,1,lx1) !int divperp of BG current
do ix1=1,lx1
J1(ix1,:,:)=1d0/x%h2(ix1,1:lx2,1:lx3)/x%h3(ix1,1:lx2,1:lx3)* &
(J1(ix1,:,:)) !FAC AT TOP ASSUMED TO BE ZERO
end do
else !non-inverted grid (logical bottom is the lowest altitude - so integrate normy)
if (debug) print *, 'Non-inverted grid detected - integrating logical bottom to top to compute FAC...'
J1=integral3D1(divJperp,x,1,lx1) !int divperp of BG current
do ix1=1,lx1
J1(ix1,:,:)=1d0/x%h2(ix1,1:lx2,1:lx3)/x%h3(ix1,1:lx2,1:lx3)* &
(-1d0*J1(ix1,:,:)) !FAC AT THE BOTTOM ASSUMED TO BE ZERO
end do
end if
end if
E1=J1/sig0
!-------
else !we resolved the field line (either 2D solve or full 3D) so just differentiate normally
!-------
Phi=-1d0*Phi
E1=grad3D1(Phi,x,1,lx1,1,lx2,1,lx3) !no haloing required since x1-derivative
Phi=-1d0*Phi
J1=sig0*E1
!-------
end if
!!!!!!!!!
!R-------
if (debug) then
print *, 'Max topside FAC (abs. val.) computed to be: ',maxval(abs(J1(1,:,:))) !ZZZ - this rey needsz to be current at the "top"
print *, 'Max polarization J2,3 (abs. val.) computed to be: ',maxval(abs(J2pol)), &
maxval(abs(J3pol))
! print *, 'Max conduction J2,3 (abs. val.) computed to be: ',maxval(abs(J2)), &
! maxval(abs(J3))
print *, 'Max conduction J2,3 computed to be: ',maxval(J2), &
maxval(J3)
print *, 'Min conduction J2,3 computed to be: ',minval(J2), &
minval(J3)
print *, 'Max conduction J1 (abs. val.) computed to be: ',maxval(abs(J1))
print *, 'flagswap: ',flagswap
endif
!R-------
!-------
!GRAND TOTAL FOR THE CURRENT DENSITY: TOSS IN POLARIZATION CURRENT SO THAT OUTPUT FILES ARE CONSISTENT
J1=J1+J1pol
J2=J2+J2pol
J3=J3+J3pol
!-------
! if (t>11) then
! open(newunit=utrace, form='unformatted', access='stream',file='Phiall.raw8', status='replace', action='write')
! write(utrace) Phiall
! close(utrace)
! error stop 'DEBUG'
! end if
end procedure potential_root_mpi_curv
end submodule potential_root