Unexpected behavior with 2D M4 solver

[tj9726]

Hi, I found a possible bug in the 2D M4 solver.
The script I ran was

import math
# -------------------
# MY PYTHON VARIABLES
# -------------------

micrometer  = 7.75701890
femtosecond = 2.32549576

a0        = 9.79333075
pduration = 160.*femtosecond/math.sqrt(2.*math.log(2.))
pwidth    = 40.*femtosecond
pcenter   = 80.*femtosecond/math.sqrt(2.*math.log(2.))

nx = 1280    # number of cells
dx = 60.*micrometer/nx   # spatial resolution
nt = 10800  # 108000 # number of time steps
dt = 540.*femtosecond/nt # 5400.*femtosecond/nt
npatch = 128 # number of patches
output = nt/18  #nt/180 # output every 30 femtosecond

# --------------------------------------
# SMILEI's VARIABLES (DEFINED IN BLOCKS)
# --------------------------------------

Main(
    geometry = "2Dcartesian",
    
    interpolation_order = 4,
    interpolator = "wt",
    
    number_of_timesteps = nt,
    timestep = dt,
    
    cell_length = [dx, dx],
    number_of_cells  = [nx, nx],
    
    number_of_patches = [npatch, npatch],
    patch_arrangement = "hilbertian",
    
    EM_boundary_conditions = [["silver-muller"]],
    maxwell_solver = "M4"
)

LoadBalancing(
    initial_balance = True,
    every = 150,
    cell_load = 1.,
    frozen_particle_load = 0.1
)

LaserGaussian2D(
    box_side        = "xmin",
    a0              = a0, # normalized amplitude
    omega           = 1.,
    focus           = [50.*micrometer, 50.*micrometer], # coordinates of laser focus
    waist           = 10.*micrometer,
    incidence_angle = 0.,
    time_envelope   = tgaussian(duration=pduration, fwhm=pwidth, center=pcenter)
)
LaserGaussian2D(
    box_side        = "ymax",
    a0              = a0, # normalized amplitude
    omega           = 1.,
    focus           = [10.*micrometer, 10.*micrometer], # coordinates of laser focus
    waist           = 10.*micrometer,
    incidence_angle = 0.,
    time_envelope   = tgaussian(duration=pduration, fwhm=pwidth, center=pcenter)
)

### DIAGNOSTICS

DiagScalar(
    every = output
)

DiagFields(
    every = output,
    fields = ["Bx_m","By_m","Bz_m","Ex","Ey","Ez"] 
)

The Bz plot I got is

This is different than the expected behavior because there should be two identical laser pulses entering from Xmin and Ymax.
By changing the solver to Yee maxwell_solver = "Yee", I got the following result

This is the expected result.

I looked at the code src/ElectroMagnSolver/MF_Solver2D_M4.cpp and found that special treatment in Ymin and Ymax is absent.
In MF_Solver2D_Bouchard.cpp, there are lines after comment //at Ymin+dy - treat using simple discretization of the curl. (and for Ymax too).
Also, both Y and Z boundaries seem to be treated specially in MF_Solver3D_M4.

It is just a simple guess, but could it be one of the reasons for the unexpected behavior?

[Z10Frank]

Hello, the M4 solver tries to mitigate the lowest order numerical Cherenkov radiation for a plasma drifting along x (https://www.sciencedirect.com/science/article/abs/pii/S0021999120301625?via%3Dihub) so I don't know if it will work well with a laser propagating perpendicularly.

[tj9726]

Thanks,
I will use other solvers when with laser.

Is the reason for this limitation known?
I also tested with a point source electric field with M4, and it propagated isotropically.

Is there anything specific to the laser module?

[mccoys]

The M4 solver was added in #486
Maybe @yingchaolu will have an answer to your questions

[Z10Frank]

Hello, I don't understand your point source was static? In that case the field was computed with the Poisson solver.

[tj9726]

Hi. This is what I tried. There seems to be no problem with light propagating in the y-direction, at least for the internal region.
Could the problem be related to something with laser boundary conditions?

import math
# -------------------
# MY PYTHON VARIABLES
# -------------------

nx = 1280
dx = 1.
nt = 1000
dt = 0.5
npatch = 32
output = nt/10

# --------------------------------------
# SMILEI's VARIABLES (DEFINED IN BLOCKS)
# --------------------------------------

Main(
    geometry = "2Dcartesian",
    
    number_of_timesteps = nt,
    timestep = dt,
    
    cell_length = [dx, dx],
    number_of_cells  = [nx, nx],
    
    number_of_patches = [npatch, npatch],
    patch_arrangement = "hilbertian",
    
    EM_boundary_conditions = [["silver-muller"]],
    maxwell_solver = "M4"
)

ExternalField(
    field = "Bz",
    profile = gaussian(1., xfwhm = 10., yfwhm = 10.)
)

DiagFields(
    every = output,
    fields = ["Bx_m","By_m","Bz_m","Ex","Ey","Ez"]
)

[Z10Frank]

Interesting, the isotropic propagation seems to work as expected.
For the case with two laser pulses, to understand if the boundary conditions are the problem, first I would try to dump a more frequent output and see with slide() the evolution of the simulation, to check if something unexpected happens at the border.

[tj9726]

This is the time evolution I got (the color bar range is not the same as the figure above).

1. The downward propagating laser is weak from the beginning.
2. The propagation in the internal region seems correct.
3. The downward propagating laser reflects at the Ymin boundary, whereas the rightward propagating laser exits the simulation box without reflection (expected behavior with SM boundary)

[Z10Frank]

Thank you, my theory at the moment is that the boundary conditions are not adapted to the M4 solver.
Can you try these tests?

• rotate the laser polarisation by 90 degrees, e.g. y polarisation to z polarisation and vice versa (your namelist has y polarised lasers it seems?)
• try using only one laser in the simulation and see if something changes in its propagation
• using another boundary condition, e.g. the Perfectly Matched Layers

[tj9726]

Hi,

I tried different interpolation orders (2 and 4), different boundary conditions (SM and PML), and different polarization (polarization_phi = 0 and pi/2) with 1 downward propagating laser (name list at the bottom).
Results were almost the same as the cases above:
Small amplitude from the entrance and spurious reflection at ymin, ymax boundaries.
I also think there is something wrong with the y boundary conditions of the M4 solver.

In MF_Solver2D_Yee.cpp, the y index runs from 1 to ny_d-2:

        for( unsigned int y = 1; y < ny_d - 1; ++y ) {
            Bz2D[x * ny_d + y] += dt_ov_dy * ( Ex2D[x * ny_p + y] - Ex2D[x * ny_p + y - 1] ) -
                                  dt_ov_dx * ( Ey2D[x * ny_d + y] - Ey2D[( x - 1 ) * ny_d + y] );
        }

This is the same for the Cowan solver, and the Bouchard solver updates y=1 and ny_d-2 as special cases.

Whereas in MF_Solver2D_M4.cpp, the y index is from 2 to ny_d-3:

        for (unsigned int j=2 ; j<ny_d-2 ; j++) {
            (*Bz2D)(i,j) += Ay * ((*Ex2D)(i,j)-(*Ex2D)(i,j-1))
                          + By * ((*Ex2D)(i+1,j)-(*Ex2D)(i+1,j-1) + (*Ex2D)(i-1,j)-(*Ex2D)(i-1,j-1))
                          + Dy * ((*Ex2D)(i,j+1)-(*Ex2D)(i,j-2))
                          + Ax * ((*Ey2D)(i-1,j)-(*Ey2D)(i,j))
                          + Bx * ((*Ey2D)(i-1,j+1)-(*Ey2D)(i,j+1) + (*Ey2D)(i-1,j-1)-(*Ey2D)(i,j-1))
                          + Dx * ((*Ey2D)(i-2,j)-(*Ey2D)(i+1,j));
        }

j =1 and ny_d-2 are not updated in the MF solver for M4 (and Grassi solver as well?).
Could this be causing the problem at the boundary (only at ymin and ymax and not at borders of patches)?
However, I do not know whether these indices are designed to be treated here or in a boundary condition function.

This is the name list I used

import math

micrometer  = 7.75701890
femtosecond = 2.32549576

a0        = 9.79333075
pduration = 160.*femtosecond/math.sqrt(2.*math.log(2.))
pwidth    = 40.*femtosecond
pcenter   = 80.*femtosecond/math.sqrt(2.*math.log(2.))

nx = 1280
dx = 60.*micrometer/nx
nt = 10800
dt = 540.*femtosecond/nt
npatch = 128
output = nt/18

Main(
    geometry = "2Dcartesian",
    
    interpolation_order = 4,
    #interpolation_order = 2,
    
    number_of_timesteps = nt,
    timestep = dt,
    
    cell_length = [dx, dx],
    number_of_cells  = [nx, nx],
    
    number_of_patches = [npatch, npatch],
    patch_arrangement = "hilbertian",
    
    EM_boundary_conditions = [["silver-muller"]],
    #EM_boundary_conditions = [["PML"]],

    maxwell_solver = "M4"
    #maxwell_solver = "Yee"
)

LaserGaussian2D(
    box_side         = "ymax",
    a0               = a0,
    omega            = 1.,
    focus            = [10.*micrometer, 10.*micrometer],
    waist            = 10.*micrometer,
    incidence_angle  = 0.,
    polarization_phi = 0.,
    #polarization_phi = math.pi/2.,
    ellipticity      = 0.,
    time_envelope    = tgaussian(duration=pduration, fwhm=pwidth, center=pcenter)
)

DiagScalar(
    every = output
)

DiagFields(
    every = output,
    fields = ["Bx_m","By_m","Bz_m","Ex","Ey","Ez"]
)

[Z10Frank]

I don't know why those indices are like this in the M4 solver.
Have you tried to make the y index run as Yee, recompile and launch again the simulation?

[tj9726]

The problems with a laser can be solved by adding these lines (which I copied from the Bouchard solver) to the M4 solver.
Also, I found that the problem is not unique to the M4.
Other solvers with the same y index (such as Grassi) also have the same problem and can be fixed by adding the same lines.
I think all of these solvers (and possibly their 3D counterparts) need fixes.
However, I do not know the code well enough to say anything conclusive.

    // at Ymin+dy - treat using simple discretization of the curl (will be overwritten if not at the ymin-border)
    for (unsigned int i=0 ; i<nx_p ; i++) {
        (*Bx2D)(i,1) += dt_ov_dy * ( (*Ez2D)(i,0) - (*Ez2D)(i,1) );
    }
    // at Ymax-dy - treat using simple discretization of the curl (will be overwritten if not at the ymax-border)
    for (unsigned int i=0 ; i<nx_p ; i++) {
        (*Bx2D)(i,ny_d-2) += dt_ov_dy * ( (*Ez2D)(i,ny_d-3) - (*Ez2D)(i,ny_d-2) );
    }
    // at Ymin+dy - treat using simple discretization of the curl (will be overwritten if not at the ymin-border)
    for (unsigned int i=2 ; i<nx_d-2 ; i++) {
        (*Bz2D)(i,1) += dt_ov_dx * ( (*Ey2D)(i-1,1) - (*Ey2D)(i,1)   )
        +               dt_ov_dy * ( (*Ex2D)(i,1) - (*Ex2D)(i,0) );
    }
    // at Ymax-dy - treat using simple discretization of the curl (will be overwritten if not at the ymax-border)
    for (unsigned int i=2 ; i<nx_d-2 ; i++) {
        (*Bz2D)(i,ny_d-2) += dt_ov_dx * ( (*Ey2D)(i-1,ny_d-2) - (*Ey2D)(i,ny_d-2)   )
        +                    dt_ov_dy * ( (*Ex2D)(i,ny_d-2) - (*Ex2D)(i,ny_d-3) );
    }
  

[Z10Frank]

Thank you! We'll test these changes and in case add them to the code.