Commit f05d0c1d authored by Sauerwein Nick Jacob's avatar Sauerwein Nick Jacob
Browse files

last version of Benedikt

parents
import numpy as np
c = 299792458
#accelerator parts
def rotation(phi):
return np.array([[np.cos(phi),0,np.sin(phi),0,0,0],
[0,np.cos(phi),0,np.sin(phi),0,0],
[-np.sin(phi),0,np.cos(phi),0,0,0],
[0,-np.sin(phi),0,np.cos(phi),0,0],
[0,0,0,0,1,0],
[0,0,0,0,0,1]])
def driftSpace(para,setup):
L = para[0]
beta_0 = setup[1]
gamma_0 = setup[2]
#print L
dri = np.array([[1,L,0,0,0,0],
[0,1,0,0,0,0],
[0,0,1,L,0,0],
[0,0,0,1,0,0],
[0,0,0,0,1,L/(beta_0*gamma_0)**2],
[0,0,0,0,0,1]])
#print ("dri",dri)
return dri.T
def diPole(para,setup):
P_0 = setup[0]
q = setup[4]
m = setup[3]
beta_0 = setup[1]
gamma_0 = setup[2]
w = q/P_0*para[1]
L = para[0]
r = rotation(para[2])
r_i = rotation(-para[2])
K1 = -q/P_0*para[1]*np.tan(para[3]*L*q*para[1]/(2*np.pi*P_0))
fringe = np.array([[1,0,0,0,0,0],
[-K1,1,0,0,0,0],
[0,0,1,0,0,0],
[0,0,K1,1,0,0],
[0,0,0,0,1,0],
[0,0,0,0,0,1]])
di = np.array([[np.cos(w*L),np.sin(w*L)/w,0,0,0,(1-np.cos(w*L))/(w*beta_0)],
[-w*np.sin(w*L),np.cos(w*L),0,0,0,np.sin(w*L)/beta_0],
[0,0,1,L,0,0],
[0,0,0,1,0,0],
[-np.sin(w*L)/beta_0,-(1-np.cos(w*L))/(w*beta_0),0,0,1,L/(beta_0*gamma_0)**2-(w*L-np.sin(w*L))/(w*beta_0**2)],
[0,0,0,0,0,1]])
#print ("di",di)
return (np.dot(np.dot(np.dot(np.dot(r , fringe) , di) , fringe) , r_i)).T
def quadroPole(para,setup):
beta_0 = para[1]
P_0 = setup[0]
q = setup[4]
beta_0 = setup[1]
gamma_0 = setup[2]
w = np.sqrt(q/P_0*para[1])
#print ('qp',w)
L = para[0]
r = rotation(para[2])
r_i = rotation(-para[2])
#print L
quad = np.array([[np.cos(w*L),np.sin(w*L)/w,0,0,0,0],
[-w*np.sin(w*L),np.cos(w*L),0,0,0,0],
[0,0,np.cosh(w*L),np.sinh(w*L)/w,0,0],
[0,0,w*np.sinh(w*L),np.cosh(w*L),0,0],
[0,0,0,0,1,L/(beta_0*gamma_0)**2],
[0,0,0,0,0,1]])
#print ("quad",r)
return np.dot(np.dot(r , quad) , r_i).T
README
In order to get the results presented in the Semesterthesis by NICK SAUERWEIN follow the next steps.
1) Import the right gcc module and python:
module load gcc/gcc-4.7.2
alias qq='qstat -f -u "*" '
alias qbd='qstat -f -q prime_bd.q -u "*"'
export PATH=$HOME/miniconda2/bin/:$PATH
2) Set the simulation parameters a desired:
- open warp_script.py with editor and set the parameters
- lower the parameters N_steps and Nz by the same factor to get less resoluting
CAUTION: the results are then not physical anymore (numerical dispersion)
3) Run the simulation:
- set number of cores in run.sge:
line 4: #$ -pe orte 128 (for 128 cores)
line 15: ARGS=' warp_script.py -p 4 4 8' (4 * 4 * 8 = 128 spatial decomposition)
- commit job to merlin cluster:
qsub run.sge
- look at progress using:
tail -f LPWA-1.o*
4) Create plots and data analyis:
- open notebook.py and set which plots to generate:
line 13: overview = True
line 14: video = True
line 15: extracted = True
line 16: spray = True
- look at progess using:
tail -f notebook_out.txt
If problems occur, please contact me under nicksauerwein@hotmail.de
This diff is collapsed.
"""
This script is an optimization script for Warp
Author: Remi Lehe (rlehe@lbl.gov)
It performs 2-parameter optimization by launching different simulation
successively and picking the next points automatically
Usage:
-----
Modify the parameters below and run
`python optimize.py`
or for long optimization on a remote cluster
`nohup python optimize.py &`
To work properly, this script requires the files
control_sim.py and template to be in the same directory
"""
import os
import numpy as np
from scipy import constants
from opmd_viewer import OpenPMDTimeSeries
from sim_control import run_sim, prepare_sim
# ----------
# Parameters
# ----------
# System on which to run the simulation
code = "warp"
machine = "merlin_parallel"
# The number of procs in the x and z directions
nx = 2
nz = 3
# ----------------------
# The objective function
# ----------------------
def objective_func( X ) :
"""
Objective function
Creates a simulation, runs it and evaluates the output
Parameters
----------
X : 1darray
An array of shape (2,)
The first component corresponds to the log of the density in m-3
The second component the a0
Returns
-------
A single scalar : the maximum of the Ez field on axis at iteration 400
"""
# Parameters
sim_name = "runwithn_{:.3f}_w0_{:.3f}_reldens_{:.3f}".format(np.log(X[2]),np.log(X[1]),X[2])
n = X[3] # Density in m-3
reldens = X[2]
a0 = X[0]
w0 = X[1]
if os.path.exists( sim_name ) == False:
# Create a new simulation directory, with the proper input script
prepare_sim( sim_name, a0, w0,reldens , n, code, machine, nx, nz )
# Launch the simulation
result_directory = run_sim( sim_name, machine )
# Extract the maximum of the Ez field
ts = OpenPMDTimeSeries( os.path.join( sim_name, "diags/hdf5" ) )
# Return the result
return ts
"""
This is a restart script for Warp.
It can restart a simulation that was run with the script warp_script.py
Usage:
-----
- If warp_script.py was run in serial (python warp_script.py), then type
python restart_after_warp_script.py
- It warp_script.py was run in parallel (mpirun -np N python warp_script.py):
mpirun -np N python restart_after_warp_script.py
(where `N` should be replaced by the number of cores used, and should
be the same for the initial simulation and for the restart simulation.)
"""
from warp import *
import numpy as np
# Reload the simulation from the dump files
dump_name = 'warp_script011000'
# - Single proc case
if npes == 1:
restart( dump_name + '.dump')
# - Parallel case
else:
restart( dump_name )
# Proceed with 500 steps
N_steps = 11000
n_stepped=0
while n_stepped < N_steps:
step(10)
n_stepped = n_stepped + 10
#!/bin/bash
#$ -cwd
#$ -N LPWA-1
#$ -pe orte 64
#$ -q prime_bdl.q
#$ -v LD_LIBRARY_PATH
#$ -l s_rt=93:59:00,h_rt=94:00:00
###################################################
#
MPIEXEC=~sauerwein_n/miniconda2/bin/mpirun
CMD=~sauerwein_n/miniconda2/bin/python
#
###################################################
# The simulation to run with mpirun python
#
ARGS=' warp_script.py -p 4 4 4'
#
export OMP_NUM_THREADS=1
export OMPI_MCA_btl='self,sm,openib'
#
rm -rf diags *.pe *.po
#
##############
# BEGIN DEBUG
# Print the SGE environment on master host:
echo "================================================================"
echo "=== SGE job JOB_NAME=$JOB_NAME JOB_ID=$JOB_ID"
echo "================================================================"
echo DATE=`date`
echo HOSTNAME=`hostname`
echo PWD=`pwd`
echo "NSLOTS=$NSLOTS"
echo "PE_HOSTFILE=$PE_HOSTFILE"
cat $PE_HOSTFILE
echo "================================================================"
echo "Running environment:"
env
echo "================================================================"
echo "Loaded environment modules:"
module list 2>&1
echo
# Check that the libraries are available (on the master host):
echo "ldd $CMD"
ldd $CMD
echo "LD_LIBRARY_PATH=$LD_LIBRARY_PATH"
# Check the number of threads used by OpenMP:
echo "OMP_NUM_THREADS=$OMP_NUM_THREADS"
# END DEBUG
##############
# The MPI command to run under the control of SGE (PE orte):
MPICMD="$MPIEXEC -np $NSLOTS $CMD $ARGS"
echo "Command to run:"
echo "$MPICMD"
echo
$MPICMD
###################################################
#!/bin/bash
#$ -cwd
#$ -N LPWA-1
#$ -pe orte 128
#$ -v LD_LIBRARY_PATH
#$ -l s_rt=23:59:00,h_rt=24:00:00
###################################################
#
MPIEXEC=$HOME/miniconda2/bin/mpirun
CMD=$HOME/miniconda2/bin/python
#
###################################################
# The simulation to run with mpirun python
#
ARGS=' restart_after_warp_script.py -p 4 4 8'
#
export OMP_NUM_THREADS=1
export OMPI_MCA_btl='self,sm,openib'
#
rm -rf diags *.pe *.po
#
##############
# BEGIN DEBUG
# Print the SGE environment on master host:
echo "================================================================"
echo "=== SGE job JOB_NAME=$JOB_NAME JOB_ID=$JOB_ID"
echo "================================================================"
echo DATE=`date`
echo HOSTNAME=`hostname`
echo PWD=`pwd`
echo "NSLOTS=$NSLOTS"
echo "PE_HOSTFILE=$PE_HOSTFILE"
cat $PE_HOSTFILE
echo "================================================================"
echo "Running environment:"
env
echo "================================================================"
echo "Loaded environment modules:"
module list 2>&1
echo
# Check that the libraries are available (on the master host):
echo "ldd $CMD"
ldd $CMD
echo "LD_LIBRARY_PATH=$LD_LIBRARY_PATH"
# Check the number of threads used by OpenMP:
echo "OMP_NUM_THREADS=$OMP_NUM_THREADS"
# END DEBUG
##############
# The MPI command to run under the control of SGE (PE orte):
MPICMD="$MPIEXEC -np $NSLOTS $CMD $ARGS"
echo "Command to run:"
echo "$MPICMD"
echo
$MPICMD
###################################################
# ------------------
# Simulation control
# ------------------
import os
import time
import re
def prepare_sim(sim_name, a0, w0,reldens, n, code, machine, nx=1, nz=1) :
"""
Prepare a simulation directory.
Create a script with the proper values of the parameters.
Parameters
----------
sim_name: str
The name of the directory where the simulation will be run
n: double
The electron density in m-3
a0: double
The a0 in the simulation
code: str
Typically 'warp'
machine: str
The type of machine on which to run the simulation
- "local_serial" for a serial job on interactive shell
- "local_parallel" for a parallel job on interactive shell
nx, nz: int
The number of procs along each direction
"""
# Create the new directory
print '\nNew simulation!'
print 'Creating simulation {:s}'.format(sim_name)
if os.path.exists( sim_name ) == False:
os.mkdir( sim_name )
# Setup the simulation script
# ---------------------------
# Read the template simulation file
template_file = 'template/warp_script.py'
with open(template_file) as f :
template_txt = f.read()
# Replace the correct values and write the new file
# Important: The values are given in the same order as they appear
# in the template script (as %f)
script_txt = template_txt %(a0, w0, reldens, n)
with open( os.path.join( sim_name, "warp_script.py" ), 'w' ) as f:
f.write(script_txt)
def run_sim( sim_name, machine ) :
"""
Run the simulation and wait for the result
Parameter:
----------
sim_name: str
The name of the simulation
machine: str
The type of machine on which to run the simulation
- "local_serial" for a serial job on interactive shell (not batch)
- "local_parallel" for a parallel job on interactive shell
- "edison" for a batch job on edison
Returns:
--------
A string indicating the directory where the results are stored
"""
# Change directory to the simulation directory
os.chdir( sim_name )
# Run the script
# - Local machine
if machine == "local_serial":
err = os.system( 'ipython warp_script.py' )
elif machine == "local_parallel":
err = os.system( 'mpirun -np 10 python2.7 warp_script.py -p 2 1 5' )
elif machine == "merlin_parallel":
import pysftp
cinfo = {'host':'merlinl02.psi.ch', 'username':'sauerwein_n', 'password':'Shadow374Fields!'}
con = pysftp.Connection(**cinfo)
con.chdir('simulations')
if con.exists(sim_name):
raise ValueError('Simulation '+sim_name+' already exsits! '+str(con.listdir()))
con.mkdir( sim_name )
con.chdir( sim_name )
con.put('warp_script.py')
con.execute('load_warp \n mpirun -np 10 python warp_script.py -p 2 1 5')
else:
raise ValueError("Unsupported system: {:s}".format(machine))
# Handle code crashing
if err !=0:
raise ValueError(
"Last simulation crashed, with code {:d}".format(err))
# Change back
os.chdir("../")
# Indicate where the results are stored
return( sim_name )
"""
This is a typical input script that runs a simulation of
laser-wakefield acceleration using Warp 2D / Circ.
Usage
-----
- Modify the parameters below to suit your needs
- Type "python -i warp_script.py" in a terminal
- When the simulation finishes, the python session will *not* quit.
Therefore the simulation can be continued by running step()
Otherwise, one can just type exit()
"""
# Import warp-specific packages
from warp_init_tools import *
import numpy as np
# -----------------------------------------------------------------------------
# Parameters (Modify the values below to suit your needs)
# -----------------------------------------------------------------------------
# General parameters
# ------------------
# Dimension of simulation ("3d", "circ", "2d", "1d")
dim = "circ"
# Number of azimuthal modes beyond m=0, for "circ" (not used for "2d" and "3d")
circ_m = 2
# Total number of timesteps in the simulation
N_steps = 24000
# Whether to run the simulation interactively (0:off, 1:on)
interactive = 0
#Density ramp parameters
start = 0.e-6
ramp1 = 1.e-4
ramp2 = 1.e-4
ramp3 = 1.e-4
plateau1 = 0.8e-4
plateau2 = 0.5e-3
# Simulation box
# --------------
# Number of grid cells in the longitudinal direction
Nz = 1650
# Number of grid cells in transverse direction (represents Nr in "circ")
Nx = 400
# Number of grid cells in the 3rd dimension (not used for "2d" and "circ")
Ny = 50
# Dimension of the box in longitudinal direction (meters)
zmin = -4.e-5
zmax = 0.6e-5
# Dimension of the box in transverse direction (box ranges from -xmax to xmax)
xmax = 60.e-6
# Dimension of the box in 3rd direction (not used for "2d" and "circ")
ymax = 15.e-6
# Field boundary conditions (longitudinal and transverse respectively)
f_boundz = openbc
f_boundxy = openbc
if dim == "circ":
f_boundxy = dirichlet
# Particles boundary conditions (longitudinal and transverse respectively)
p_boundz = absorb
p_boundxy = absorb
# Moving window (0:off, 1:on)
use_moving_window = 1
# Speed of the moving window (ignored if use_moving_window = 0)
v_moving_window = clight
# Diagnostics
# -----------
# Period of diagnostics (in number of timesteps)
diag_period = 300
# Whether to write the fields
write_fields = 1
# Whether to write the particles
write_particles = 1
# Whether to write the diagnostics in parallel
parallel_output = False
# Numerical parameters
# --------------------
# Field solver (0:Yee, 1:Karkkainen on EF,B, 3:Lehe)
stencil = 0
# Particle shape (1:linear, 2:quadratic, 3:cubic)
depos_order = 2
# Gathering mode (1:from cell centers, 4:from Yee mesh)
efetch = 1
# Particle pusher (0:Boris, 1:Vay)
particle_pusher = 1
# Current smoothing parameters
# ----------------------------
# Turn current smoothing on or off (0:off; 1:on)
use_smooth = 1
# Number of passes of smoother and compensator in each direction (x, y, z)
npass_smooth = array([[ 0 , 0 ], [ 0 , 0 ], [ 1 , 1 ]])
# Smoothing coefficients in each direction (x, y, z)
alpha_smooth = array([[ 0.5, 3.], [ 0.5, 3.], [0.5, 3./2]])
# Stride in each direction (x, y, z)
stride_smooth = array([[ 1 , 1 ], [ 1 , 1 ], [ 1 , 1 ]])
# Laser parameters
# ----------------
# Initialize laser (0:off, 1:on)
use_laser = 1
# Position of the antenna (meters)
laser_source_z = -0.5e-5
# Polarization angle with respect to the x axis (rad)
laser_polangle = pi/2
# Laser file:
# When using a laser profile that was experimentally
# measured, provide a string with the path to an HDF5 laser file,
# otherwise provide None and a Gaussian pulse will be initialized
laser_file = None
laser_file_energy = 2. # When using a laser file, energy in Joule of the pulse
# Gaussian pulse:
# Laser amplitude at focus
Energy=30.
laser_a0 = 0.787662/np.sqrt(40./Energy)
# Waist at focus (meters)
laser_w0 = 0.000008
# Length of the pulse (length from the peak to 1/e of the amplitude ; meters)
laser_ctau = 15e-15*clight
# Initial position of the centroid (meters)
laser_z0 = -2 * laser_ctau
# Focal position
laser_zfoc = start+ramp1+plateau1+ramp2+1e-5
# Plasma macroparticles
# ---------------------
# Initialize some preexisting plasmas electrons (0:off, 1:on)
# (Can be used in order to neutralize pre-ionized ions, if any,
# or in order to simulate a plasma without having to initialize ions)
use_preexisting_electrons = 0
# Initialize plasma ions (0:off, 1:on)
use_ions = 1
# Number of macroparticles per cell in each direction
# In Circ, nppcelly is the number of particles along the
# azimuthal direction. Use a multiple of 4*circ_m
plasma_nx = 2
plasma_ny = 8
plasma_nz = 2
reldens = 0.5
levelofion = 8
# Plasma content and profile
# --------------------------
# Reference plasma density (in number of particles per m^3)
n_plasma = 1.4e25/levelofion/reldens
# Relative density of the preexisting electrons (relative to n_plasma)
rel_dens_preexisting_electrons = 0
# The different elements used. (Only used if use_ions is different than 0.)
# relative_density is the density relative to n_plasma.
# q_start is the ionization state of the ions at the beginning of the simulation
# q_max is the maximum ionization state
# If q_start is not equal to q_max, ionization between states will be computed.
ion_states = { 'Argon': {'relative_density':1., 'q_start':0, 'q_max':11 } }
# Positions between which the plasma is initialized
# (Transversally, the plasma is initialized between -plasma_xmax and
# plasma_xmax, along x, and -plasma_ymax and plasma_ymax along y)
plasma_zmin = 1.e-6
plasma_zmax = 1500.e-6
plasma_xmax = xmax
plasma_ymax = ymax
# Define your own profile and profile parameters below
def dens_profile(z):
if z < start or z > ramp1+plateau1+ramp2+plateau2+ramp3:
return 0.
if z < ramp1:
return z/ramp1
if z < ramp1+plateau1:
return 1.
if z < ramp1+plateau1+ramp2:
return (reldens-1)/ramp2*(z-(ramp1+plateau1))+1
if z < ramp1+plateau1+ramp2+plateau2:
return reldens
if z < ramp1+plateau1+ramp2+plateau2+ramp3:
return reldens/ramp3*(ramp1+plateau1+ramp2+plateau2+ramp3-z)
return 0.