Warwick-Plasma/epoch

Pukhov Solver - Boundary Conditions Query

Opened this issue · 2 comments

Hello, I am wanting to use the Pukhov solver with dt_multiplier =1.0, but when I have run my input deck, it appears to be giving me a c*dt/dx value of around 0.707 when I run using the following boundary conditions. I also tried dt_multiplier = 0.999 and it still gave around 0.707 when run

bc_x_min = simple_laser
bc_x_max = simple_outflow
bc_y_min = simple_outflow
bc_y_max = simple_outflow

I then tried the deck with cpml boundary conditions

bc_x_min = cpml_laser
bc_x_max = cpml_outflow
bc_y_min = cpml_outflow
bc_y_max = cpml_outflow

Along with dt_mulitplier = 1.0, this gave c*dt/dx = 1.0 as I was hoping, but the laser field does not appear to propagate normally and additionally, it appears my target is just very unstable and 'blowing up'.

This may be down to my lack of understanding of cpml boundaries, but is it possible to either ensure c*dt/dx =1.0 for simple outflow boundaries or can I tailor the other parameters of cpml to ensure this runs smoothly? I have my input deck attached here in case it is something else that I've missed in it that's causing the issue.

-------------------- EPOCH 2D -------------------

begin:constant

runtime = 120 * femto #total runtime of simulation

n_rotation = 45 * (pi/180) #plasma angle of rotation about its defined centre

x_res = 62.5 / micron #resolution - no. cells per micron in x
y_res = 62.5 / micron #resolution in y

grid_x_min = 0 * micron #minimum grid position in x
grid_y_min = 0 * micron #minimum grid position in y

xlength = 30 * micron #length of grid in x
ylength = 30 * micron #length of grid in y

las_lambda = 0.8 * micron #laser wavelength

las_polanglex = 0 #polarisation angles for p-polarised light
#las_polangley = pi/2

las_cep = pi/2 #carrier envelope phase of laser

las_fwhm = 4 * micron #full width at half maximum
las_fwhmt = 7 * femto #pulse duration in fs
las_a0 = 10 #normalised vector potential value
las_focus_time = 55 * femto #time to reach focus - should be at least ~ 3*fwhmt +10fs

n_crit = critical(2 * pi * c / las_lambda) #critical density of plasma

n_critx = 15 * micron #position of critical surface about which the plasma rotates
n_crity = 15 * micron

n_max = 200 * n_crit #target density
n_min = 0.05 * n_crit #minimum ramp density (to save computing time)
n_electron_per_cell = 100 #number of electrons per cell

scale_length = las_lambda / 15 #scale length - ensure resolution is smaller than this
n_thickness = 0.5 * micron #thickness of the target
n_length = 20 * micron #length of the target
n_width_shift = -0.5 * micron #shift parallel to surface
n_temp = 0 #target temp in ev

x_c = scale_length * loge ( n_max / n_crit ) #distance between n = n_crit and n = n_max in exponential ramp

x_rot = (x - n_critx) * cos(n_rotation) + (y - n_crity) * sin(n_rotation)
y_rot = -(x - n_critx) * sin(n_rotation) + (y - n_crity) * cos(n_rotation)

las_focus_x = n_critx #want laser to focus onto critical point of target surface
las_focus_y = n_crity

las_width = 0.894322 * las_fwhm #beam waist wo
las_widtht = 0.894322 * las_fwhmt #pulse duration for 1/e radii

z_lasx = x #axial direction of x_min laser
r_lasx = y #radial direction of x_min laser

las_lambdax = las_lambda
las_kx = (2 * pi ) / las_lambdax #wavenumber of the driver, corrected for the angle of propagation
las_rayleighx = pi * las_width ^ 2 / las_lambdax #rayleigh range - corrected for the wavelength
las_wx = las_width * sqrt( 1 + ( z_lasx / las_rayleighx ) ^2 ) #spot size - may need corrected for angle

las_rcx = z_lasx * (1 + (las_rayleighx / (z_lasx + 1e-20)) ^ 2) #radius of curvature
las_rc_shiftx = r_lasx ^ 2 / (2 * las_rcx + 1e-20) #phasefront shift due to the curvature of the wavefronts
las_gouyx = if(ndims eq 1,0,atan(z_lasx / las_rayleighx)) #gouy phase shift acquired by a gaussian beam as it propagates
las_intensity_focusx = las_a0 ^ 2 * 1.368167127e6 / las_lambdax ^ 2 #the intensity of the beam as it reaches focus

end:constant

begin:control

t_end = runtime #runtime
particle_tstart = 0 *femto #time to being particle pushing

nx = xlength * x_res #global no. cells in x
ny = ylength * y_res #global no. cells in y

x_min = grid_x_min #min grid positions
y_min = grid_y_min

x_max= x_min + xlength #max grid positions
y_max = y_min + ylength

stdout_frequency = 1

maxwell_solver = pukhov

dt_multiplier = 1.0
end:control

begin:boundaries
bc_x_min = cpml_laser
bc_x_max = cpml_outflow
bc_y_min = cpml_outflow
bc_y_max = cpml_outflow

end:boundaries

begin:laser
boundary = x_min
intensity_w_cm2 = (las_width / las_wx ) ^ (ndims -1 ) * las_intensity_focusx

profile = gauss(r_lasx, y_max/2, las_wx) * gauss(time- las_rc_shiftx/c ,las_focus_time + z_lasx/c, las_widtht)

phase = las_cep - las_kx * z_lasx - las_kx * las_rc_shiftx + las_gouyx - 2 * pi * c / las_lambdax * las_focus_time

lambda = las_lambdax
pol_angle = las_polanglex

end:laser

begin:species

name = electron
charge = -1.0
mass = 1.0
nparticles_per_cell = n_electron_per_cell

density = if ( (x_rot lt (x_c + n_thickness)) and (x_rot gt x_c) and (abs(y_rot - n_width_shift) lt (n_length/2)), n_max, 0.0)
density = if ( (x_rot lt x_c) and (abs(y_rot - n_width_shift) lt (n_length/2)), n_max * exp( (x_rot - x_c) / scale_length ), density(electron) )

density_min = n_min
density_max = n_max
temp_ev=n_temp
end:species

begin:species
name = O8
charge = 8
mass = 16 * 1836.2
nparticles_per_cell = n_electron_per_cell / 5

density = density(electron) / 30 * 2
density_min = n_min / 30 * 2
density_max = n_max / 30 * 2

temp_ev=n_temp
end:species

begin:species
name = Si14
charge = 14
mass = 28 * 1836.2
nparticles_per_cell = n_electron_per_cell / 5

density = density(electron) / 30
density_min = n_min / 30
density_max = n_max / 30

temp_ev=n_temp
end:species

begin:output

dt_snapshot = 10e-15
grid = always
ex = always
ey = always
ez = always
bx = always
by = always
bz = always
number_density = always + species

end:output

Hey @HollyHuddle,

I don't think this is a bug. The EPOCH timestep is set by the subroutine set_dt in the file housekeeping/setup.F90. Here we see that the standard Yee solver has a more restrictive time-step, which takes the form resulting in your 0.707 value. The important lines (in 2D) are:

    IF (maxwell_solver == c_maxwell_solver_yee) THEN
      ! Default maxwell solver with field_order = 2, 4 or 6
      ! cfl is a function of field_order
      dt = cfl * dx * dy / SQRT(dx**2 + dy**2) / c

    ELSE
      ! R. Lehe, PhD Thesis (2014)
      ! A. Pukhov, Journal of Plasma Physics 61, 425-433 (1999)
      dt = MIN(dx, dy) / c
    END IF

    IF (any_open) THEN
      dt_solver = dx * dy / SQRT(dx**2 + dy**2) / c
      dt = MIN(dt, dt_solver)
    END IF

We only use the dt=dx/c timestep in the Pukhov algorithm, but this is overwritten if we have an open/laser boundary in the simulation, as we use a traditional Yee FDTD solver for the boundaries.

Unfortunately this isn't something I know much about. I've checked a few references, and I'm not sure why the EPOCH time-step takes this form. I can leave the issue open as I look into this further.

Cheers,
Stuart

Hi @Status-Mirror,

Apologies for the late reply here. It would be great to have the Pukhov solver working for those boundary conditions but I understand it is no simple task. Thanks very much for taking a look.