Skip to content

Commit a2d610f

Browse files
glesurglesur
committed
- use the left domain as the reference domain for EMFs (instead of
 averaging) to be coherent with BXs exchange routine - check that vector potential follow the same boundary logic as the EMF (left domain is the reference) Note: BCs on the vector potential is only apply after boundary conditions, since EMFs boundary conditions ensure that the vector potential will always be consistent later. Todo: - enforce these boundary conditions with vector potential with periodic boundary conditions and without MPI (as is done for EMFs) - check that BXs normal is also consistent when MPI is off with periodic BCs
1 parent dfd52ca commit a2d610f

8 files changed

Lines changed: 77 additions & 86 deletions

File tree

‎src/dataBlock/dataBlock.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -367,6 +367,7 @@ real DataBlock::ComputeTimestep() {
367367
void DataBlock::DeriveVectorPotential() {
368368
if constexpr(DefaultPhysics::mhd) {
369369
#ifdef EVOLVE_VECTOR_POTENTIAL
370+
hydro->emf->EnforceVectorPotentialBoundary(hydro->Ve);
370371
hydro->emf->ComputeMagFieldFromA(hydro->Ve, hydro->Vs);
371372
#endif
372373
}

‎src/fluid/boundary/axis.cpp‎

Lines changed: 9 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -22,10 +22,9 @@ void Axis::ShowConfig() {
2222
}
2323

2424

25-
void Axis::SymmetrizeEx1Side(int jref) {
25+
void Axis::SymmetrizeEx1Side(int jref, IdefixArray3D<real> Ex1) {
2626
#if DIMENSIONS == 3
2727

28-
IdefixArray3D<real> Ex1 = this->ex;
2928
IdefixArray1D<real> Ex1Avg = this->Ex1Avg;
3029

3130
idefix_for("Ex1_ini",0,data->np_tot[IDIR],
@@ -63,9 +62,7 @@ void Axis::SymmetrizeEx1Side(int jref) {
6362
// Hence, we enforce a regularisation of Ex3 for consistancy.
6463

6564

66-
void Axis::RegularizeEx3side(int jref) {
67-
IdefixArray3D<real> Ex3 = this->ez;
68-
65+
void Axis::RegularizeEx3side(int jref, IdefixArray3D<real> Ex3) {
6966
idefix_for("Ex3_Regularise",0,data->np_tot[KDIR],0,data->np_tot[IDIR],
7067
KOKKOS_LAMBDA(int k,int i) {
7168
Ex3(k,jref,i) = 0.0;
@@ -133,18 +130,20 @@ void Axis::RegularizeCurrentSide(int side) {
133130

134131
// Average the Emf component along the axis
135132

136-
void Axis::RegularizeEMFs() {
133+
void Axis::RegularizeEMFs(IdefixArray3D<real> ex,
134+
IdefixArray3D<real> ey,
135+
IdefixArray3D<real> ez) {
137136
idfx::pushRegion("Axis::RegularizeEMFs");
138137

139138
if(this->axisLeft) {
140139
int jref = data->beg[JDIR];
141-
SymmetrizeEx1Side(jref);
142-
RegularizeEx3side(jref);
140+
SymmetrizeEx1Side(jref, ex);
141+
RegularizeEx3side(jref, ez);
143142
}
144143
if(this->axisRight) {
145144
int jref = data->end[JDIR];
146-
SymmetrizeEx1Side(jref);
147-
RegularizeEx3side(jref);
145+
SymmetrizeEx1Side(jref, ex);
146+
RegularizeEx3side(jref, ez);
148147
}
149148

150149
idfx::popRegion();

‎src/fluid/boundary/axis.hpp‎

Lines changed: 5 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -26,15 +26,16 @@ class Axis {
2626
public:
2727
template <typename Phys>
2828
explicit Axis(Boundary<Phys> *); // Initialisation
29-
void RegularizeEMFs(); // Regularize the EMF sitting on the axis
29+
void RegularizeEMFs(IdefixArray3D<real>, IdefixArray3D<real>, IdefixArray3D<real>);
30+
// Regularize the EMF sitting on the axis
3031
void RegularizeCurrent(); // Regularize the currents along the axis
3132
void EnforceAxisBoundary(int side); // Enforce the boundary conditions (along X2)
3233
void RegularizeBX2s(); // Regularize BX2s on the axis
3334
void ShowConfig();
3435

35-
36-
void SymmetrizeEx1Side(int); // Symmetrize on a specific side (internal method)
37-
void RegularizeEx3side(int); // Regularize Ex3 along the axis (internal method)
36+
// Internal methods
37+
void SymmetrizeEx1Side(int, IdefixArray3D<real>); // Symmetrize on a specific side
38+
void RegularizeEx3side(int, IdefixArray3D<real>); // Regularize Ex3 along the axis
3839
void RegularizeCurrentSide(int); // Regularize J along the axis (internal method)
3940
void FixBx2sAxis(int side); // Fix BX2s on the axis using the field around it (internal)
4041
void FixBx2sAxisGhostAverage(int side); //Fix BX2s on the axis using the average of neighbouring
@@ -76,9 +77,6 @@ class Axis {
7677
IdefixArray1D<int> symmetryVc;
7778
IdefixArray1D<int> symmetryVs;
7879

79-
IdefixArray3D<real> ex;
80-
IdefixArray3D<real> ey;
81-
IdefixArray3D<real> ez;
8280
IdefixArray4D<real> J;
8381

8482
IdefixArray4D<real> Vc;
@@ -94,11 +92,6 @@ Axis::Axis(Boundary<Phys> *boundary) {
9492
Vc = boundary->Vc;
9593
Vs = boundary->Vs;
9694
J = boundary->fluid->J;
97-
if constexpr(Phys::mhd) {
98-
ex = boundary->fluid->emf->ex;
99-
ey = boundary->fluid->emf->ey;
100-
ez = boundary->fluid->emf->ez;
101-
}
10295

10396
data = boundary->data;
10497
haveMHD = Phys::mhd;

‎src/fluid/constrainedTransport/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ target_sources(idefix
55
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/constrainedTransport.hpp
66
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/EMFexchange.hpp
77
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/enforceEMFBoundary.hpp
8+
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/EnforceVectorPotentialBoundary.hpp
89
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/evolveMagField.hpp
910
PUBLIC ${CMAKE_CURRENT_LIST_DIR}/evolveVectorPotential.hpp
1011
)

‎src/fluid/constrainedTransport/EMFexchange.hpp‎

Lines changed: 18 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -14,16 +14,18 @@
1414

1515
#ifdef WITH_MPI
1616
template<typename Phys>
17-
void ConstrainedTransport<Phys>::ExchangeAll() {
18-
if(data->mygrid->nproc[IDIR]>1) this->ExchangeX1();
19-
if(data->mygrid->nproc[JDIR]>1) this->ExchangeX2();
20-
if(data->mygrid->nproc[KDIR]>1) this->ExchangeX3();
17+
void ConstrainedTransport<Phys>::ExchangeAll(IdefixArray3D<real> ex,
18+
IdefixArray3D<real> ey,
19+
IdefixArray3D<real> ez) {
20+
if(data->mygrid->nproc[IDIR]>1) this->ExchangeX1(ey,ez);
21+
if(data->mygrid->nproc[JDIR]>1) this->ExchangeX2(ex,ez);
22+
if(data->mygrid->nproc[KDIR]>1) this->ExchangeX3(ex,ey);
2123
}
2224

2325

2426
// Exchange EMFs in X1
2527
template<typename Phys>
26-
void ConstrainedTransport<Phys>::ExchangeX1() {
28+
void ConstrainedTransport<Phys>::ExchangeX1(IdefixArray3D<real> ey, IdefixArray3D<real> ez) {
2729
idfx::pushRegion("Emf::ExchangeX1");
2830

2931

@@ -34,8 +36,6 @@ void ConstrainedTransport<Phys>::ExchangeX1() {
3436

3537
IdefixArray1D<real> BufferLeft=BufferSendX1[faceLeft];
3638
IdefixArray1D<real> BufferRight=BufferSendX1[faceRight];
37-
IdefixArray3D<real> ey=this->ey;
38-
IdefixArray3D<real> ez=this->ez;
3939

4040

4141
// If MPI Persistent, start receiving even before the buffers are filled
@@ -61,7 +61,6 @@ void ConstrainedTransport<Phys>::ExchangeX1() {
6161

6262
idefix_for("LoadBufferX1Emfz",kbeg,kend,jbeg,jend+1,
6363
KOKKOS_LAMBDA (int k, int j) {
64-
BufferLeft( (j-jbeg) + (k-kbeg)*(ny+1) ) = ez(k,j,ileft);
6564
BufferRight( (j-jbeg) + (k-kbeg)*(ny+1) ) = ez(k,j,iright);
6665
}
6766
);
@@ -70,7 +69,6 @@ void ConstrainedTransport<Phys>::ExchangeX1() {
7069

7170
idefix_for("LoadBufferX1Emfy",kbeg,kend+1,jbeg,jend,
7271
KOKKOS_LAMBDA (int k, int j) {
73-
BufferLeft( (j-jbeg) + (k-kbeg)*ny + Vsindex ) = ey(k,j,ileft);
7472
BufferRight( (j-jbeg) + (k-kbeg)*ny + Vsindex ) = ey(k,j,iright);
7573
}
7674
);
@@ -91,29 +89,21 @@ void ConstrainedTransport<Phys>::ExchangeX1() {
9189
BufferLeft=BufferRecvX1[faceLeft];
9290
BufferRight=BufferRecvX1[faceRight];
9391

94-
// We average the edge emfs zones
92+
// Erase the emf with the one coming from the left process
93+
9594
idefix_for("StoreBufferX1Emfz",kbeg,kend,jbeg,jend+1,
9695
KOKKOS_LAMBDA (int k, int j) {
9796
if(lbound == internal || lbound == periodic) {
98-
ez(k,j,ileft) = HALF_F*(
99-
BufferLeft( (j-jbeg) + (k-kbeg)*(ny+1) ) + ez(k,j,ileft) );
100-
}
101-
if(rbound == internal || rbound == periodic) {
102-
ez(k,j,iright) = HALF_F*(
103-
BufferRight( (j-jbeg) + (k-kbeg)*(ny+1) ) + ez(k,j,iright) );
97+
ez(k,j,ileft) = BufferLeft( (j-jbeg) + (k-kbeg)*(ny+1));
10498
}
10599
});
100+
106101
#if DIMENSIONS == 3
107102
Vsindex = (ny+1)*nz;
108103
idefix_for("StoreBufferX1Emfy",kbeg,kend+1,jbeg,jend,
109104
KOKKOS_LAMBDA (int k, int j) {
110105
if(lbound == internal || lbound == periodic) {
111-
ey(k,j,ileft) = HALF_F*(
112-
BufferLeft( (j-jbeg) + (k-kbeg)*ny +Vsindex) + ey(k,j,ileft) );
113-
}
114-
if(rbound == internal || rbound == periodic) {
115-
ey(k,j,iright) = HALF_F*(
116-
BufferRight( (j-jbeg) + (k-kbeg)*ny +Vsindex) + ey(k,j,iright) );
106+
ey(k,j,ileft) = BufferLeft( (j-jbeg) + (k-kbeg)*ny +Vsindex);
117107
}
118108
});
119109
#endif
@@ -124,7 +114,7 @@ void ConstrainedTransport<Phys>::ExchangeX1() {
124114

125115
// Exchange EMFs in X2
126116
template<typename Phys>
127-
void ConstrainedTransport<Phys>::ExchangeX2() {
117+
void ConstrainedTransport<Phys>::ExchangeX2(IdefixArray3D<real> ex, IdefixArray3D<real> ez) {
128118
idfx::pushRegion("Emf::ExchangeX2");
129119

130120
// Load the buffers with data
@@ -133,8 +123,6 @@ void ConstrainedTransport<Phys>::ExchangeX2() {
133123
[[maybe_unused]] int nz;
134124
IdefixArray1D<real> BufferLeft=BufferSendX2[faceLeft];
135125
IdefixArray1D<real> BufferRight=BufferSendX2[faceRight];
136-
IdefixArray3D<real> ex=this->ex;
137-
IdefixArray3D<real> ez=this->ez;
138126

139127
// If MPI Persistent, start receiving even before the buffers are filled
140128
double tStart = MPI_Wtime();
@@ -158,7 +146,6 @@ void ConstrainedTransport<Phys>::ExchangeX2() {
158146

159147
idefix_for("LoadBufferX2Emfz",kbeg,kend,ibeg,iend+1,
160148
KOKKOS_LAMBDA (int k, int i) {
161-
BufferLeft( (i-ibeg) + (k-kbeg)*(nx+1) ) = ez(k,jleft,i);
162149
BufferRight( (i-ibeg) + (k-kbeg)*(nx+1) ) = ez(k,jright,i);
163150
}
164151
);
@@ -167,7 +154,6 @@ void ConstrainedTransport<Phys>::ExchangeX2() {
167154

168155
idefix_for("LoadBufferX1Emfx",kbeg,kend+1,ibeg,iend,
169156
KOKKOS_LAMBDA (int k, int i) {
170-
BufferLeft( (i-ibeg) + (k-kbeg)*nx + Vsindex ) = ex(k,jleft,i);
171157
BufferRight( (i-ibeg) + (k-kbeg)*nx + Vsindex ) = ex(k,jright,i);
172158
}
173159
);
@@ -191,25 +177,15 @@ void ConstrainedTransport<Phys>::ExchangeX2() {
191177
idefix_for("StoreBufferX2Emfz",kbeg,kend,ibeg,iend+1,
192178
KOKKOS_LAMBDA (int k, int i) {
193179
if(lbound == internal || lbound == periodic) {
194-
ez(k,jleft,i) = HALF_F*(
195-
BufferLeft( (i-ibeg) + (k-kbeg)*(nx+1) ) + ez(k,jleft,i) );
196-
}
197-
if(rbound == internal || rbound == periodic) {
198-
ez(k,jright,i) = HALF_F*(
199-
BufferRight( (i-ibeg) + (k-kbeg)*(nx+1) ) + ez(k,jright,i) );
180+
ez(k,jleft,i) = BufferLeft( (i-ibeg) + (k-kbeg)*(nx+1) );
200181
}
201182
});
202183
#if DIMENSIONS == 3
203184
Vsindex = (nx+1)*nz;
204185
idefix_for("StoreBufferX1Emfy",kbeg,kend+1,ibeg,iend,
205186
KOKKOS_LAMBDA (int k, int i) {
206187
if(lbound == internal || lbound == periodic) {
207-
ex(k,jleft,i) = HALF_F*(
208-
BufferLeft( (i-ibeg) + (k-kbeg)*nx +Vsindex) + ex(k,jleft,i) );
209-
}
210-
if(rbound == internal || rbound == periodic) {
211-
ex(k,jright,i) = HALF_F*(
212-
BufferRight( (i-ibeg) + (k-kbeg)*nx +Vsindex) + ex(k,jright,i) );
188+
ex(k,jleft,i) = BufferLeft( (i-ibeg) + (k-kbeg)*nx +Vsindex);
213189
}
214190
});
215191
#endif
@@ -220,7 +196,7 @@ void ConstrainedTransport<Phys>::ExchangeX2() {
220196

221197
// Exchange EMFs in X3
222198
template<typename Phys>
223-
void ConstrainedTransport<Phys>::ExchangeX3() {
199+
void ConstrainedTransport<Phys>::ExchangeX3(IdefixArray3D<real> ex, IdefixArray3D<real> ey) {
224200
idfx::pushRegion("Emf::ExchangeX3");
225201

226202

@@ -229,8 +205,6 @@ void ConstrainedTransport<Phys>::ExchangeX3() {
229205
int nx,ny;
230206
IdefixArray1D<real> BufferLeft=BufferSendX3[faceLeft];
231207
IdefixArray1D<real> BufferRight=BufferSendX3[faceRight];
232-
IdefixArray3D<real> ex=this->ex;
233-
IdefixArray3D<real> ey=this->ey;
234208

235209
int Vsindex = 0;
236210

@@ -259,15 +233,13 @@ void ConstrainedTransport<Phys>::ExchangeX3() {
259233

260234
idefix_for("LoadBufferX3Emfx",jbeg,jend+1,ibeg,iend,
261235
KOKKOS_LAMBDA (int j, int i) {
262-
BufferLeft( (i-ibeg) + (j-jbeg)*nx ) = ex(kleft,j,i);
263236
BufferRight( (i-ibeg) + (j-jbeg)*nx ) = ex(kright,j,i);
264237
}
265238
);
266239
Vsindex = nx*(ny+1);
267240

268241
idefix_for("LoadBufferX3Emfy",jbeg,jend,ibeg,iend+1,
269242
KOKKOS_LAMBDA (int j, int i) {
270-
BufferLeft( (i-ibeg) + (j-jbeg)*(nx+1) + Vsindex ) = ey(kleft,j,i);
271243
BufferRight( (i-ibeg) + (j-jbeg)*(nx+1) + Vsindex ) = ey(kright,j,i);
272244
}
273245
);
@@ -290,29 +262,17 @@ void ConstrainedTransport<Phys>::ExchangeX3() {
290262
idefix_for("StoreBufferX3Emfx",jbeg,jend+1,ibeg,iend,
291263
KOKKOS_LAMBDA (int j, int i) {
292264
if(lbound == internal || lbound == periodic) {
293-
ex(kleft,j,i) = HALF_F*(
294-
BufferLeft( (i-ibeg) + (j-jbeg)*nx ) + ex(kleft,j,i) );
295-
}
296-
if(rbound == internal || rbound == periodic) {
297-
ex(kright,j,i) = HALF_F*(
298-
BufferRight( (i-ibeg) + (j-jbeg)*nx ) + ex(kright,j,i) );
265+
ex(kleft,j,i) = BufferLeft( (i-ibeg) + (j-jbeg)*nx );
299266
}
300267
});
301268

302269
Vsindex = nx*(ny+1);
303270
idefix_for("StoreBufferX3Emfy",jbeg,jend,ibeg,iend+1,
304271
KOKKOS_LAMBDA (int j, int i) {
305272
if(lbound == internal || lbound == periodic) {
306-
ey(kleft,j,i) = HALF_F*(
307-
BufferLeft( (i-ibeg) + (j-jbeg)*(nx+1) + Vsindex ) + ey(kleft,j,i) );
308-
}
309-
if(rbound == internal || rbound == periodic) {
310-
ey(kright,j,i) = HALF_F*(
311-
BufferRight( (i-ibeg) + (j-jbeg)*(nx+1) + Vsindex ) + ey(kright,j,i) );
273+
ey(kleft,j,i) = BufferLeft( (i-ibeg) + (j-jbeg)*(nx+1) + Vsindex );
312274
}
313275
});
314-
315-
316276
idfx::popRegion();
317277
}
318278

‎src/fluid/constrainedTransport/constrainedTransport.hpp‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -101,13 +101,14 @@ class ConstrainedTransport {
101101
// Routines for evolving the magnetic potential (only available when EVOLVE_VECTOR_POTENTIAL)
102102
void EvolveVectorPotential(real, IdefixArray4D<real> &);
103103
void ComputeMagFieldFromA(IdefixArray4D<real> &Vein, IdefixArray4D<real> &Vsout);
104+
void EnforceVectorPotentialBoundary(IdefixArray4D<real> &Vein); // Enforce BCs on A
104105

105106
#ifdef WITH_MPI
106107
// Exchange surface EMFs to remove interprocess round off errors
107-
void ExchangeAll();
108-
void ExchangeX1();
109-
void ExchangeX2();
110-
void ExchangeX3();
108+
void ExchangeAll(IdefixArray3D<real> ex, IdefixArray3D<real> ey, IdefixArray3D<real> ez);
109+
void ExchangeX1(IdefixArray3D<real> ey, IdefixArray3D<real> ez);
110+
void ExchangeX2(IdefixArray3D<real> ex, IdefixArray3D<real> ez);
111+
void ExchangeX3(IdefixArray3D<real> ex, IdefixArray3D<real> ey);
111112
#endif
112113

113114
private:
@@ -449,6 +450,7 @@ void ConstrainedTransport<Phys>::ShowConfig() {
449450
#include "calcRiemannEmf.hpp"
450451
#include "EMFexchange.hpp"
451452
#include "enforceEMFBoundary.hpp"
453+
#include "enforceVectorPotentialBoundary.hpp"
452454
#include "evolveMagField.hpp"
453455
#include "evolveVectorPotential.hpp"
454456

‎src/fluid/constrainedTransport/enforceEMFBoundary.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -31,14 +31,14 @@ void ConstrainedTransport<Phys>::EnforceEMFBoundary() {
3131
this->data->hydro->emfBoundaryFunc(*data, data->t);
3232

3333
if(this->data->hydro->haveAxis) {
34-
this->data->hydro->boundary->axis->RegularizeEMFs();
34+
this->data->hydro->boundary->axis->RegularizeEMFs(this->ex, this->ey, this->ez);
3535
}
3636

3737
#ifdef ENFORCE_EMF_CONSISTENCY
3838
#ifdef WITH_MPI
3939
// This average the EMFs at the domain surface with immediate neighbours
4040
// to ensure the EMFs exactly match
41-
this->ExchangeAll();
41+
this->ExchangeAll(this->ex, this->ey, this->ez);
4242
#endif
4343
#endif
4444

‎src/fluid/constrainedTransport/enforceVectorPotentialBoundary.hpp‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
// ***********************************************************************************
2+
// Idefix MHD astrophysical code
3+
// Copyright(C) Geoffroy R. J. Lesur <geoffroy.lesur@univ-grenoble-alpes.fr>
4+
// and other code contributors
5+
// Licensed under CeCILL 2.1 License, see COPYING for more information
6+
// ***********************************************************************************
7+
8+
#ifndef FLUID_CONSTRAINEDTRANSPORT_ENFORCEVECTORPOTENTIALBOUNDARY_HPP_
9+
#define FLUID_CONSTRAINEDTRANSPORT_ENFORCEVECTORPOTENTIALBOUNDARY_HPP_
10+
#include "constrainedTransport.hpp"
11+
12+
template<typename Phys>
13+
void ConstrainedTransport<Phys>::EnforceVectorPotentialBoundary(IdefixArray4D<real> &Vein) {
14+
idfx::pushRegion("Emf::EnforceVectorPotentialBoundary");
15+
16+
auto Ax1 = Kokkos::subview(Vein, IDIR, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
17+
auto Ax2 = Kokkos::subview(Vein, JDIR, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
18+
auto Ax3 = Kokkos::subview(Vein, KDIR, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
19+
20+
if(this->hydro->haveAxis) {
21+
this->hydro->boundary->axis->RegularizeEMFs(Ax1, Ax2, Ax3);
22+
}
23+
24+
#ifdef ENFORCE_EMF_CONSISTENCY
25+
#ifdef WITH_MPI
26+
// This average the vector potential at the domain surface with immediate neighbours
27+
// to ensure the vector potentials exactly match
28+
29+
this->ExchangeAll(Ax1, Ax2, Ax3);
30+
#endif
31+
#endif
32+
33+
idfx::popRegion();
34+
}
35+
#endif // FLUID_CONSTRAINEDTRANSPORT_ENFORCEVECTORPOTENTIALBOUNDARY_HPP_

0 commit comments

Comments
 (0)