Pukhov Solver - Boundary Conditions Query

[HollyHuddle]

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 ~ 3fwhmt +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

[Status-Mirror]

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

[HollyHuddle]

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.