Open Shihang-Hu opened 5 months ago
Status-Mirror
commented
5 months ago Hi @Shihang-Hu,
These boundary conditions were implemented before my time as a developer, but are we sure they're wrong? From my understanding, a perfectly conducting boundary has zero electric field perpendicular to the boundary, and zero magnetic field parallel to the boundary.
That is, on x_min and x_max, we have zero Ex, which is what the code is currently doing. Am I missing something?
Cheers, Stuart
Shihang-Hu
commented
5 months ago Hi @Status-Mirror. I have different understanding, a perfectly conducting boundary means that electric and magnetic fields outside boundary are zero (as an example, the internal electromagnetic field of superconductivity is zero).
According to Maxwell's equations, $$\oint\mathbf{E\cdot}d\mathbf{l}=-\int\frac{\partial \mathbf{B}}{\partial t}\cdot d\mathbf{S}$$ taking a infinitely small loop, the right side of the equation is zero. And we will know that the parallel electric fields $\mathbf{E{t}}$ inside and outside boundary are qeual. Because of zero electric fields outside boundary, we will have $\mathbf{E{t}}=0$ or $\mathbf{n\times E}=0$ inside boundary. It is obvious because the superconductivity is a equipotential body and the electric fields is perpendicular to the boundary.
That is, on x_min and x_max, we have $E{y}=0$ and $E{z}=0$.
Meanwhile, according $$\oint\mathbf{E\cdot}d\mathbf{S}=\int\frac{\rho}{\varepsilon{0}}dV$$ taking a infinitely thin cylinder, with zero electric fields outside boundary, we have $\mathbf{n\cdot E}=\rho{S}/\varepsilon{0}$ inside boundary. But the surface density of charge $\rho{S}$ is hard to define in simulation, I think we can assume the volume density of charge $\rho$ on the boundary is zero and we get $\mathbf{\nabla\cdot E}=\rho/\varepsilon{0}=0$, and because of zero parallel electric fields we have $\partial E{n}/\partial n=0$.
That is, on x_min and x_max, we have $\partial E_{x}/\partial x=0$.
Similarly, according $$\oint\mathbf{B\cdot}d\mathbf{S}=0$$ taking a infinitely thin cylinder, we will know that the perpendicular magnetic fields $\mathbf{B{n}}$ inside and outside boundary are qeual. Because of zero magnetic fields outside boundary, we will have $\mathbf{B{n}}=0$ or $\mathbf{n\cdot B}=0$ inside boundary.
That is, on x_min and x_max, we have zero $B_{x}=0$.
Meanwhile, according $$\oint\mathbf{B\cdot}d\mathbf{l}=-\int\mu{0}(\mathbf{J}+\varepsilon{0}\frac{\partial \mathbf{E}}{\partial t})\cdot d\mathbf{S}$$ taking a infinitely small loop, with zero magnetic fields outside boundary, we have $\mathbf{n\times B}=\mu{0}\mathbf{J{l}}$ inside boundary.Similarly, line density of current $\mathbf{J{l}}$ is hard to define in simulation, I think we can use $\mathbf{\nabla\cdot B}=0$ , and because of zero perpendicular magnetic fields we have $\mathbf{\nabla\cdot B{t}}=0$.
That is, on x_min and x_max, we have $\partial B{y}/\partial y=0$ and $\partial B{z}/\partial z=0$.
Perfectly conducting boundary is that:
Perfect Electrical Conductor (PEC) boundary condition is $$\mathbf{n\times E}=0$$ Perfect Magnetic Conductor (PMC) boundary condition is $$\mathbf{n\cdot B}=0$$
But for perpendicular electric fields $\mathbf{E{n}}$ and parallel magnetic fields $\mathbf{B{t}}$ we need more information or make more assumptions.
Shihang-Hu
commented
4 months ago Dear @Status-Mirror . Now I am sure they are wrong. I have used epoch2d to simulation magnetic reconnection. Perfectly conducting boundary condition is applied in the y direction and periodic boundary condition is applied in the x direction, which is typical in magnetic reconnection simulation. And the simulation results are wrong outrageously, but when I modify perfectly conducting boundary as mentioned above, the results are right.
Status-Mirror
commented
4 months ago Hey @Shihang-Hu,
It would be useful to us if you could send the reconnection input deck for our own testing. Could you paste it as a response here inside ``` quotations to skip the GitHub formatting?
Cheers, Stuart
Shihang-Hu
commented
4 months ago Dear @Status-Mirror . input deck:
begin:constant
di = 3.0e-5
wci = 3e8/(40*di)
inv_wci = wci^(-1)
m_ratio = 400
T_ratio = 5.0
bgd_ratio = 0.1
delta_ratio = 0.5
psi1 = 0.1
Lx_ratio = 20
Ly_ratio = 10
dx_ratio = 2.5e-2
dy_ratio = 2.5e-2
dt_ratio = 2.5e-4
npart_Drift = 10
npart_Background = npart_Drift*(bgd_ratio*Ly_ratio)/(2*delta_ratio)
delta = delta_ratio*di
B0 = (m_ratio*me*wci)/qe
n0 = (m_ratio*me)/(mu0*(qe*di)^2)
Te0 = (B0^2)/(2*mu0*n0*kb*(1+T_ratio))
Pde = (2*me*kb*Te0)/(qe*B0*delta)
Pdi = -m_ratio*T_ratio*Pde
end:constant
begin:control
nsteps = 30/dt_ratio
nx = Lx_ratio/dx_ratio
ny = Ly_ratio/dy_ratio
x_min = -Lx_ratio*di/2
x_max = Lx_ratio*di/2
y_min = -Ly_ratio*di/2
y_max = Ly_ratio*di/2
dt_multiplier = dt_ratio*sqrt(dx_ratio^2+dy_ratio^2)/(dx_ratio*dy_ratio)*c/(wci*di)
maxwell_solver = yee # default
stdout_frequency = 1
print_constants = T
end:control
begin:fields
# initial Harris equilibrium
# bx=(1-psi1*di/(delta*cosh(x/delta)*cosh(y/delta)))*B0*tanh(y/delta)
# by=(psi1*B0*di/delta)*tanh(x/delta)/(cosh(x/delta)*cosh(y/delta))
bx=B0*tanh(y/delta)-((pi*psi1*B0)/(Ly_ratio))*cos((2*pi*x)/(Lx_ratio*di))*sin((pi*y)/(Ly_ratio*di))
by=((2*pi*psi1*B0)/(Lx_ratio))*sin((2*pi*x)/(Lx_ratio*di))*cos((pi*y)/(Ly_ratio*di))
# bz=B0
end:fields
begin:boundaries
# periodic conduct reflect open
bc_x_min_field = periodic
bc_x_max_field = periodic
bc_x_min_particle = periodic
bc_x_max_particle = periodic
bc_y_min_field = conduct
bc_y_max_field = conduct
bc_y_min_particle = reflect
bc_y_max_particle = reflect
end:boundaries
begin:species
name = ElectronDrift
charge = -1.0
mass = 1.0
npart = npart_Drift*nx*ny
number_density = n0/(cosh(y/delta))^2
temperature = Te0
drift_z = Pde
end:species
begin:species
name = IonDrift
charge = 1.0
mass = m_ratio
npart = npart_Drift*nx*ny
number_density = n0/(cosh(y/delta))^2
temperature = T_ratio * Te0
drift_z = Pdi
end:species
begin:species
name = ElectronBackground
charge = -1.0
mass = 1.0
npart = npart_Background*nx*ny
number_density = bgd_ratio * n0
temperature = Te0 * 3
end:species
begin:species
name = IonBackground
charge = 1.0
mass = m_ratio
npart = npart_Background*nx*ny
number_density = bgd_ratio * n0
temperature = T_ratio * Te0 * 3
end:species
begin:collisions # collision between particles
use_collisions = F # T -- collisional; F -- collisionless
# use_nanbu = T # use Nanbu-Perez scattering algorithm
# coulomb_log = 30 # Set the fixed Coulomb logarithm ///calculate coulomb logarithm based on local temperature and density
# collide = all
# collide = IonBackground IonBackground off # collisions between background species are omitted
# collide = IonBackground ElectronBackground off # collisions between background species are omitted
# collide = ElectronBackground ElectronBackground off # collisions between background species are omitted
end:collisions
begin:output
name = fields
nstep_average = 0.025/dt_ratio
nstep_snapshot = 0.1/dt_ratio
ex = always + average
ey = always + average
ez = always + average
bx = always
by = always
bz = always
jx = always + species + no_sum
jy = always + species + no_sum
jz = always + species + no_sum
temperature = always + species + no_sum
number_density = always + species + no_sum
end:output
# begin:output
# name = particles
# file_prefix = part
# dump_at_nsteps = 17/dt_ratio,18/dt_ratio,18.4/dt_ratio,18.8/dt_ratio,19/dt_ratio,19.2/dt_ratio,19.4/dt_ratio,19.6/dt_ratio
# id = always
# particle_grid = always
# end:output
begin:output
name = restart
file_prefix = restart
dump_at_nsteps = 15/dt_ratio
restartable = T
dump_first = F
end:output
begin:output_global
force_first_to_be_restartable = F
force_last_to_be_restartable = F
dump_last = F
end:output_global
Before I modify the perfectly conducting boundary, I get magnetic field (0040.sdf):
After I modify the perfectly conducting boundary, I get magnetic field (0060.sdf):
Shihang-Hu
commented
4 months ago Hey @Status-Mirror . Did you get the simulation results with the input deck I sent?
Status-Mirror
commented
4 months ago Hi @Shihang-Hu,
I think you might be right about this - I have another developer looking into it.
We're currently transitioning the code from FORTRAN to C++, and I've passed on the information to those writing up the new boundaries. In the meantime, I've been told to hold off updating the FORTRAN side - feel free to use the fix you've found yourself, and I'll leave this issue up for others who may be confused.
Cheers, Stuart
Shihang-Hu
commented
4 months ago Hey @Status-Mirror . I have got it. Thanks for your time.
There could be some error for perfectly conducting boundaries in SUBROUTINE efield_bcs and SUBROUTINE bfield_bcs.\ As an example, in file epoch2d\src\boundary.F90
I think it should be modified to
Thanks for your time.