search

SmileiPIC / Smilei

Particle-in-cell code for plasma simulation
https://smileipic.github.io/Smilei
337 stars 118 forks source link

Moving window in 2D #309

Closed egelfer closed 4 years ago

egelfer commented 4 years ago

Description

Dear SMILEI developers,

I have a problem with running of 2D simulation with moving window. When the window starts to move, the length of the shift is not equal to velocity_x c dt, but to the length of the simulation box. It leads to the loosing of all the particles and fields.

For example, when I use the command

for tt in timesteps: gxm+=[Diag.getXmoved(tt)/l0] print gxm

it provides [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, 140.0, 140.0, 140.0, 140.0, 140.0, 140.0, 140.0, 140.0, 140.0, 140.0, 140.0]

and in the input file (see below) grid_length=70*l0.

In the corresponding 1D simulation this problem does not appear. 

# ----------------------------------------------------------------------------------------
#                     SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI
# ----------------------------------------------------------------------------------------

import math                 
import numpy as np
from numpy import s_

l0 = 2.0*math.pi       # laser wavelength
t0 = l0                # optical cycle
Lsim = [70.*l0,64.*l0]  # length of the simulation
resx = 50.            # nb of cells in on laser wavelength
rest = math.sqrt(2.)*resx/0.95
dt   = t0/rest

tau = 6*t0

ww=4.*l0
cut=4.*t0
dd=0.*l0

edge1=30.*l0
edge2=31*l0

aa0=400.0                        

nn=56
plasma_width=4.0*l0
focusd=0.0*l0

## LASER PROPERTIES
theta  = 0.                             # angle of incidence
waist  = ww/2.                          # laser waist

time_delay    = 10.*tau           # time delay to ensure full width at 8w0 to enter the box

xsim=200.*l0

Tsim         = time_delay + xsim

ss=5

A=1.
Z=1.

c = 299792458
electron_mass = 9.10938356e-31
electron_charge = 1.60217662e-19
lambdar = 1e-6                # Normalization wavelength
wr = 2*math.pi*c/lambdar      # Normalization frequency

def time_envelope(t):
    tloc = t-time_delay
    return math.exp(-(tloc)**2/(tau/2)**2)

def carrier(t):
    return np.cos(t)

def space_envelope(y):
    return np.exp(-(y-Lsim[1]/2.)**8/waist**8)

def By_profile(y,t):
    return aa0*space_envelope(y) * time_envelope(t) * carrier(t+math.pi/2)

def Bz_profile(y,t):
    return aa0*space_envelope(y) * time_envelope(t) * carrier(t)

def porte(y,t):
    return 1.

Main(
    geometry = "2Dcartesian",

    interpolation_order = 2 ,

    cell_length = [l0/resx,l0/resx],
    grid_length  = Lsim,

    patch_arrangement = "linearized_XY",

    number_of_patches = [ 1, 128],

    #timestep_over_CFL=0.5,

    timestep=dt,

    simulation_time = Tsim,

    EM_boundary_conditions = [
        ['silver-muller'],
        ['silver-muller'],
    ],

    #random_seed = smilei_mpi_rank,
    print_every = int( Tsim/dt/100. ),
    reference_angular_frequency_SI = wr
)

LaserOffset(
    box_side               = "xmin",
    space_time_profile     = [ By_profile, Bz_profile ],
    offset                 = edge1,
    #extra_envelope          = tconstant(),
    keep_n_strongest_modes = 256,
    angle = 0
)

Species(
    name = 'ions',
    position_initialization = 'random',
    momentum_initialization = 'cold',
    particles_per_cell = int(nn),
    mass = 1836.0*A,
    charge = Z,
    number_density = trapezoidal(nn/Z,xvacuum=edge1,xplateau=edge2-edge1, xslope1=0.1*(edge2-edge1), xslope2=0.1*(edge2-edge1), yvacuum=(Lsim[1]-plasma_width)/2., yplateau=plasma_width, yslope1=0., yslope2=0.),
    boundary_conditions = [
        ["remove", "remove"],
        ["remove", "remove"],
    ],
    #time_frozen = edge1-t0
)
Species(
    name = 'electrons',
    position_initialization = 'ions',
    momentum_initialization = 'cold',
    particles_per_cell = int(nn),
    mass = 1.0,
    charge = -1.0,
    number_density = trapezoidal(nn,xvacuum=edge1,xplateau=edge2-edge1, xslope1=0.1*(edge2-edge1), xslope2=0.1*(edge2-edge1), yvacuum=(Lsim[1]-plasma_width)/2., yplateau=plasma_width, yslope1=0., yslope2=0.),
    boundary_conditions = [
        ["remove", "remove"],
        ["remove", "remove"],
    ],
    #time_frozen = edge1-t0
)

MovingWindow(
    time_start = time_delay,
    velocity_x = 0.5,
    number_of_additional_shifts = 0.,
    additional_shifts_time = 0.,
)

globalEvery = 406.

DiagScalar(every=globalEvery)

DiagFields(
    every = globalEvery,
    fields = ['Ex','Ey','Ez','Rho_ions','Rho_electrons'],
    subgrid = s_[::ss, ::ss]
)

DiagParticleBinning(
    deposited_quantity = "weight",
    every = globalEvery,
    time_average = 1,
    species = ["ions"],
    axes = [ ["ekin",    0,    120000,   2000, "edge_inclusive"],
    ["moving_x",0,Lsim[0],int(Lsim[0]/l0)],["y",0,Lsim[1],int(Lsim[1]/l0)]],
)

DiagParticleBinning(
    deposited_quantity = "weight",
    every = globalEvery,
    time_average = 1,
    species = ["electrons"],
    axes = [ ["ekin",    0,    120000,   2000, "edge_inclusive"]],
)

DiagParticleBinning(
    deposited_quantity = "weight",
    every = globalEvery,
    time_average = 1,
    species = ["electrons"],
    axes = [["px",  -60000,    60000,    1000, "edge_inclusive"]],
)

DiagParticleBinning(
    deposited_quantity = "weight",
    every = globalEvery,
    time_average = 1,
    species = ["ions"],
    axes = [ ["ekin",    0,    120000,   2000, "edge_inclusive"]],
)

DiagParticleBinning(
    deposited_quantity = "weight",
    every = globalEvery,
    time_average = 1,
    species = ["ions"],
    axes = [ ["px",  -60000,    60000,   1000, "edge_inclusive"]],
)

Steps to reproduce the problem

If relevant, provide a step-by-step guide

Parameters

mccoys commented 4 years ago

The doc says:

Each “shift” consists in removing a column of patches from the x_min border and adding a new one after the x_max border, thus changing the physical domain that the simulation represents but keeping the same box size

If your patch takes the whole domain, then all will disappear at once. You need more patches.

egelfer commented 4 years ago

Dear mccoys, thank you very much, adding more patches in x direction solved the problem!