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.