error for perfectly conducting boundaries

[Shihang-Hu]

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

   SUBROUTINE efield_bcs

      INTEGER :: i

      ! These are the MPI boundaries
      CALL field_bc(ex, ng)
      CALL field_bc(ey, ng)
      CALL field_bc(ez, ng)

      ! Perfectly conducting boundaries
      DO i = c_bd_x_min, c_bd_x_max, c_bd_x_max - c_bd_x_min
         IF (bc_field(i) == c_bc_conduct) THEN
            CALL field_clamp_zero(ex, ng, c_stagger_ex, i)
            CALL field_zero_gradient(ey, c_stagger_ey, i)
            CALL field_zero_gradient(ez, c_stagger_ez, i)
         END IF
      END DO

      DO i = c_bd_y_min, c_bd_y_max, c_bd_y_max - c_bd_y_min
         IF (bc_field(i) == c_bc_conduct) THEN
            CALL field_zero_gradient(ex, c_stagger_ex, i)
            CALL field_clamp_zero(ey, ng, c_stagger_ey, i)
            CALL field_zero_gradient(ez, c_stagger_ez, i)
         END IF
      END DO

      DO i = 1, 2*c_ndims
         ! These apply zero field boundary conditions on the edges
         IF (bc_field(i) == c_bc_clamp &
            .OR. bc_field(i) == c_bc_simple_laser &
            .OR. bc_field(i) == c_bc_simple_outflow) THEN
            CALL field_clamp_zero(ex, ng, c_stagger_ex, i)
            CALL field_clamp_zero(ey, ng, c_stagger_ey, i)
            CALL field_clamp_zero(ez, ng, c_stagger_ez, i)
         END IF

         ! These apply zero gradient boundary conditions on the edges
         IF (bc_field(i) == c_bc_zero_gradient &
            .OR. bc_field(i) == c_bc_cpml_laser &
            .OR. bc_field(i) == c_bc_cpml_outflow) THEN
            CALL field_zero_gradient(ex, c_stagger_ex, i)
            CALL field_zero_gradient(ey, c_stagger_ey, i)
            CALL field_zero_gradient(ez, c_stagger_ez, i)
         END IF
      END DO

   END SUBROUTINE efield_bcs



   SUBROUTINE bfield_bcs(mpi_only)

      LOGICAL, INTENT(IN) :: mpi_only
      INTEGER :: i

      ! These are the MPI boundaries
      CALL field_bc(bx, ng)
      CALL field_bc(by, ng)
      CALL field_bc(bz, ng)

      IF (mpi_only) RETURN

      ! Perfectly conducting boundaries
      DO i = c_bd_x_min, c_bd_x_max, c_bd_x_max - c_bd_x_min
         IF (bc_field(i) == c_bc_conduct) THEN
            CALL field_zero_gradient(bx, c_stagger_bx, i)
            CALL field_clamp_zero(by, ng, c_stagger_by, i)
            CALL field_clamp_zero(bz, ng, c_stagger_bz, i)
         END IF
      END DO

      DO i = c_bd_y_min, c_bd_y_max, c_bd_y_max - c_bd_y_min
         IF (bc_field(i) == c_bc_conduct) THEN
            CALL field_clamp_zero(bx, ng, c_stagger_bx, i)
            CALL field_zero_gradient(by, c_stagger_by, i)
            CALL field_clamp_zero(bz, ng, c_stagger_bz, i)
         END IF
      END DO

      DO i = 1, 2*c_ndims
         ! These apply zero field boundary conditions on the edges
         IF (bc_field(i) == c_bc_clamp &
            .OR. bc_field(i) == c_bc_simple_laser &
            .OR. bc_field(i) == c_bc_simple_outflow) THEN
            CALL field_clamp_zero(bx, ng, c_stagger_bx, i)
            CALL field_clamp_zero(by, ng, c_stagger_by, i)
            CALL field_clamp_zero(bz, ng, c_stagger_bz, i)
         END IF

         ! These apply zero gradient boundary conditions on the edges
         IF (bc_field(i) == c_bc_zero_gradient &
            .OR. bc_field(i) == c_bc_cpml_laser &
            .OR. bc_field(i) == c_bc_cpml_outflow) THEN
            CALL field_zero_gradient(bx, c_stagger_bx, i)
            CALL field_zero_gradient(by, c_stagger_by, i)
            CALL field_zero_gradient(bz, c_stagger_bz, i)
         END IF
      END DO

   END SUBROUTINE bfield_bcs

I think it should be modified to

   SUBROUTINE efield_bcs

      INTEGER :: i

      ! These are the MPI boundaries
      CALL field_bc(ex, ng)
      CALL field_bc(ey, ng)
      CALL field_bc(ez, ng)

      ! Perfectly conducting boundaries
      DO i = c_bd_x_min, c_bd_x_max, c_bd_x_max - c_bd_x_min
         IF (bc_field(i) == c_bc_conduct) THEN
            CALL field_zero_gradient(ex, c_stagger_ex, i)
            CALL field_clamp_zero(ey, ng, c_stagger_ey, i)
            CALL field_clamp_zero(ez, ng, c_stagger_ez, i)
         END IF
      END DO

      DO i = c_bd_y_min, c_bd_y_max, c_bd_y_max - c_bd_y_min
         IF (bc_field(i) == c_bc_conduct) THEN
            CALL field_clamp_zero(ex, ng, c_stagger_ex, i)
            CALL field_zero_gradient(ey, c_stagger_ey, i)
            CALL field_clamp_zero(ez, ng, c_stagger_ez, i)
         END IF
      END DO

      DO i = 1, 2*c_ndims
         ! These apply zero field boundary conditions on the edges
         IF (bc_field(i) == c_bc_clamp &
            .OR. bc_field(i) == c_bc_simple_laser &
            .OR. bc_field(i) == c_bc_simple_outflow) THEN
            CALL field_clamp_zero(ex, ng, c_stagger_ex, i)
            CALL field_clamp_zero(ey, ng, c_stagger_ey, i)
            CALL field_clamp_zero(ez, ng, c_stagger_ez, i)
         END IF

         ! These apply zero gradient boundary conditions on the edges
         IF (bc_field(i) == c_bc_zero_gradient &
            .OR. bc_field(i) == c_bc_cpml_laser &
            .OR. bc_field(i) == c_bc_cpml_outflow) THEN
            CALL field_zero_gradient(ex, c_stagger_ex, i)
            CALL field_zero_gradient(ey, c_stagger_ey, i)
            CALL field_zero_gradient(ez, c_stagger_ez, i)
         END IF
      END DO

   END SUBROUTINE efield_bcs



   SUBROUTINE bfield_bcs(mpi_only)

      LOGICAL, INTENT(IN) :: mpi_only
      INTEGER :: i

      ! These are the MPI boundaries
      CALL field_bc(bx, ng)
      CALL field_bc(by, ng)
      CALL field_bc(bz, ng)

      IF (mpi_only) RETURN

      ! Perfectly conducting boundaries
      DO i = c_bd_x_min, c_bd_x_max, c_bd_x_max - c_bd_x_min
         IF (bc_field(i) == c_bc_conduct) THEN
            CALL field_clamp_zero(bx, ng, c_stagger_bx, i)
            CALL field_zero_gradient(by, c_stagger_by, i)
            CALL field_zero_gradient(bz, c_stagger_bz, i)
         END IF
      END DO

      DO i = c_bd_y_min, c_bd_y_max, c_bd_y_max - c_bd_y_min
         IF (bc_field(i) == c_bc_conduct) THEN
            CALL field_zero_gradient(bx, c_stagger_bx, i)
            CALL field_clamp_zero(by, ng, c_stagger_by, i)
            CALL field_zero_gradient(bz, c_stagger_bz, i)
         END IF
      END DO

      DO i = 1, 2*c_ndims
         ! These apply zero field boundary conditions on the edges
         IF (bc_field(i) == c_bc_clamp &
            .OR. bc_field(i) == c_bc_simple_laser &
            .OR. bc_field(i) == c_bc_simple_outflow) THEN
            CALL field_clamp_zero(bx, ng, c_stagger_bx, i)
            CALL field_clamp_zero(by, ng, c_stagger_by, i)
            CALL field_clamp_zero(bz, ng, c_stagger_bz, i)
         END IF

         ! These apply zero gradient boundary conditions on the edges
         IF (bc_field(i) == c_bc_zero_gradient &
            .OR. bc_field(i) == c_bc_cpml_laser &
            .OR. bc_field(i) == c_bc_cpml_outflow) THEN
            CALL field_zero_gradient(bx, c_stagger_bx, i)
            CALL field_zero_gradient(by, c_stagger_by, i)
            CALL field_zero_gradient(bz, c_stagger_bz, i)
         END IF
      END DO

   END SUBROUTINE bfield_bcs

Thanks for your time.

[Status-Mirror]

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]

Perfectly conducting boundary

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).

Electric fields

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$.

https://eng.libretexts.org/Bookshelves/Electrical_Engineering/Electro-Optics/Book%3A_Electromagnetics_I_(Ellingson)/05%3A_Electrostatics/5.17%3A_Boundary_Conditions_on_the_Electric_Field_Intensity_(E)

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$.

Magnetic fields

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$.

https://eng.libretexts.org/Bookshelves/Electrical_Engineering/Electro-Optics/Book%3A_Electromagnetics_I_(Ellingson)/07%3A_Magnetostatics/7.10%3A_Boundary_Conditions_on_the_Magnetic_Flux_Density_(B)

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$.

Conclusion

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]

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]

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]

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]

Hey @Status-Mirror .
Did you get the simulation results with the input deck I sent?

[Status-Mirror]

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]

Hey @Status-Mirror .
I have got it. Thanks for your time.