In [1]:
import pencilnew as pcn
import os

Warning: no h5py library found.
!! ERR in diag/fixed_points.py: Dependency of h5py not fullfilled.
Warning: no h5py library found.

This tutorial will introduce you to simulation objects. They are used to most easily set up, compile, modify, access and copy simulations. Many more functionallity may be added soon.

Basics of creating simulations by copying

Simulation objects in PCN are structured as :

  • The simulation object and its methods can be found in pencilnew/sim/simulation.py
  • If not therein methods for handling simulations, examples below, are found in pencilnew/sim as well.
  • Methods which understand simulation objects as inputs, i.e. reading, calculating, etc. have the optional parameter 'sim' where you can hand over your simulation object, or 'sims' where you can hand over a list of simulation objects.

Lets start with getting a simulation object and make a copy of it.

In [2]:
SIM = pcn.get_sim('../sample_simulations/2d_streaming_instability')

# Check the identity of this simulation by using:
print('Simulation with name '+SIM.name+' is located at')
print(SIM.path)

Simulation with name 2d_streaming_instability is located at
/home/user/pencil-code/python/tutorials/sample_simulations/2d_streaming_instability

Beware, quiet=True is default. So if something seems to go wrong, set quiet=False check for problems. If the following output is to long, click on the output cell and use 'O' to hide its content

In [3]:
SIM = pcn.get_sim('../sample_simulations/2d_streaming_instability', quiet=False)

# Creating Simulation object for /home/user/pencil-code/python/tutorials/sample_simulations/2d_streaming_instability
#
# Provide parameters from Pencil code namelist files (typically param.nml
# and param2.nml).
#
# This file was generated by nl2python, so you will probably not want to
# edit it.
#

class Params:
    cvsid = '$Id: start.in 9840 2008-09-05 07:29:37Z ajohan $'
    ip = 14
    xyz0 = numpy.repeat(-5.0000000000000003E-2,3)
    xyz1 = numpy.repeat(5.0000000000000003E-2,3)
    lxyz = numpy.repeat(0.10000000000000001,3)
    lperi = numpy.repeat(True,3)
    lshift_origin = numpy.repeat(False,3)
    lshift_origin_lower = numpy.repeat(False,3)
    coord_system = 'cartesian'
    lpole = numpy.repeat(False,3)
    lequidist = numpy.repeat(True,3)
    coeff_grid = numpy.repeat(1.0,3)
    zeta_grid0 = 0.0
    grid_func = numpy.repeat('linear',3)
    xyz_star = numpy.repeat(0.0,3)
    lwrite_ic = True
    lwrite_avg1d_binary = False
    lnowrite = False
    luniform_z_mesh_aspect_ratio = False
    unit_system = 'cgs'
    unit_length = 1.0
    lmodify = False
    modify_filename = 'modify.dat'
    unit_velocity = 1.0
    unit_density = 1.0
    unit_temperature = 1.2027220936797545E-8
    unit_magnetic = 3.5449077018110318
    c_light = 29979245800.0
    g_newton = 6.6742000000000003E-8
    hbar = 1.054571596E-27
    random_gen = 'min_std'
    seed0 = 1812
    nfilter = 0
    lserial_io = False
    der2_type = 'standard'
    lread_oldsnap = False
    lread_oldsnap_nomag = False
    lread_oldsnap_nopscalar = False
    lread_oldsnap_notestfield = False
    lread_oldsnap_notestscalar = False
    ireset_tstart = 2
    tstart = 0.0
    lghostfold_usebspline = False
    lread_aux = False
    lwrite_aux = True
    pretend_lntt = False
    lprocz_slowest = True
    lcopysnapshots_exp = False
    bcx = numpy.repeat('she',5)
    bcy = numpy.repeat('p',5)
    bcz = numpy.repeat('p',5)
    r_int = 0.0
    r_ext = 3.9084999999999999E+37
    r_ref = 1.0
    rsmooth = 0.0
    r_int_border = 3.9084999999999999E+37
    r_ext_border = 3.9084999999999999E+37
    mu0 = 0.99999999999999989
    force_lower_bound = ''
    force_upper_bound = ''
    lseparate_persist = False
    ldistribute_persist = True
    lpersist = True
    fbcx1 = numpy.repeat(0.0,5)
    fbcx2 = numpy.repeat(0.0,5)
    fbcx1_2 = numpy.repeat(0.0,5)
    fbcx2_2 = numpy.repeat(0.0,5)
    fbcy1 = numpy.repeat(0.0,5)
    fbcy2 = numpy.repeat(0.0,5)
    fbcy1_1 = numpy.repeat(0.0,5)
    fbcy1_2 = numpy.repeat(0.0,5)
    fbcy2_1 = numpy.repeat(0.0,5)
    fbcy2_2 = numpy.repeat(0.0,5)
    fbcz1 = numpy.repeat(0.0,5)
    fbcz2 = numpy.repeat(0.0,5)
    fbcz1_1 = numpy.repeat(0.0,5)
    fbcz1_2 = numpy.repeat(0.0,5)
    fbcz2_1 = numpy.repeat(0.0,5)
    fbcz2_2 = numpy.repeat(0.0,5)
    fbcx_bot = numpy.repeat(0.0,5)
    fbcx_top = numpy.repeat(0.0,5)
    fbcy_bot = numpy.repeat(0.0,5)
    fbcy_top = numpy.repeat(0.0,5)
    fbcz_bot = numpy.repeat(0.0,5)
    fbcz_top = numpy.repeat(0.0,5)
    xyz_step = numpy.repeat(1.0,6)
    xi_step_frac = numpy.repeat(1.0,6)
    xi_step_width = numpy.repeat(1.5,6)
    dxi_fact = numpy.repeat(1.0,3)
    trans_width = numpy.repeat(1.0,3)
    lcylinder_in_a_box = False
    lsphere_in_a_box = False
    llocal_iso = False
    init_loops = 1
    lwrite_2d = False
    lcylindrical_gravity = False
    border_frac_x = numpy.repeat(0.0,2)
    border_frac_y = numpy.repeat(0.0,2)
    border_frac_z = numpy.repeat(0.0,2)
    lborder_hyper_diff = True
    luse_latitude = False
    lshift_datacube_x = False
    lfargo_advection = False
    yequator = 0.0
    lequatory = False
    lequatorz = False
    zequator = 0.0
    lav_smallx = False
    xav_max = 3.9084999999999999E+37
    niter_poisson = 30
    lforce_shear_bc = True
    lread_from_other_prec = False
    pipe_func = 'error_function'
    glncrosssec0 = 0.0
    crosssec_x1 = -1.0
    crosssec_x2 = 1.0
    crosssec_w = 0.10000000000000001
    lcorotational_frame = False
    rcorot = 1.0
    lproper_averages = False
    ldirect_access = False
    ltolerate_namelist_errors = False
    lyinyang = False
    lall_onesided = False
    xhe = 0.0
    mu = 1.0
    cp = 1.0
    cs0 = 1.0
    rho0 = 1.0
    gamma = 1.0
    error_cp = 9.9999999999999995E-7
    cs2top_ini = 3.9084999999999999E+37
    dcs2top_ini = 3.9084999999999999E+37
    sigmasbt = 1.0
    lanelastic_lin = False
    lcs_as_aux = False
    lcs_as_comaux = False
    fac_cs = 1.0
    isothmid = 0
    lstratz = False
    gztype = 'zero'
    gz_coeff = 0.0
    ampluu = numpy.repeat(0.0,5)
    ampl_ux = numpy.repeat(0.0,5)
    ampl_uy = numpy.repeat(0.0,5)
    ampl_uz = numpy.repeat(0.0,5)
    phase_ux = numpy.repeat(0.0,5)
    phase_uy = numpy.repeat(0.0,5)
    phase_uz = numpy.repeat(0.0,5)
    inituu = numpy.array(['zero']+list(numpy.repeat('nothing',4)))
    widthuu = 0.10000000000000001
    radiusuu = 1.0
    urand = 0.0
    urandi = 0.0
    lpressuregradient_gas = True
    uu_xz_angle = numpy.repeat(0.0,5)
    relhel_uu = 1.0
    coefuu = numpy.repeat(complex(0.0,0.0),3)
    r_omega = 0.0
    w_omega = 0.0
    uu_left = 0.0
    uu_right = 0.0
    uu_lower = 1.0
    uu_upper = 1.0
    kx_uu = 1.0
    ky_uu = 1.0
    kz_uu = 1.0
    kx_ux = numpy.repeat(0.0,5)
    ky_ux = numpy.repeat(0.0,5)
    kz_ux = numpy.repeat(0.0,5)
    kx_uy = numpy.repeat(0.0,5)
    ky_uy = numpy.repeat(0.0,5)
    kz_uy = numpy.repeat(0.0,5)
    kx_uz = numpy.repeat(0.0,5)
    ky_uz = numpy.repeat(0.0,5)
    kz_uz = numpy.repeat(0.0,5)
    uy_left = 0.0
    uy_right = 0.0
    uu_const = numpy.repeat(0.0,3)
    omega = 1.0
    u_out_kep = 0.0
    initpower = 1.0
    initpower2 = -1.6666666666666667
    cutoff = 0.0
    ncutoff = 1.0
    kpeak = 10.0
    kgaussian_uu = 0.0
    lskip_projection = False
    lno_second_ampl = True
    z1_uu = 0.0
    z2_uu = 0.0
    n_modes_uu = 0
    lcoriolis_force = True
    lcentrifugal_force = False
    ladvection_velocity = True
    lprecession = False
    omega_precession = 0.0
    alpha_precession = 0.0
    velocity_ceiling = -1.0
    loo_as_aux = False
    luut_as_aux = False
    loot_as_aux = False
    mu_omega = 0.0
    nb_rings = 0
    om_rings = numpy.repeat(0.0,5)
    gap = 0.0
    lscale_tobox = True
    ampl_omega = 0.0
    omega_ini = 0.0
    r_cyl = 1.0
    skin_depth = 0.10000000000000001
    incl_alpha = 0.0
    rot_rr = 0.0
    xsphere = 0.0
    ysphere = 0.0
    zsphere = 0.0
    neddy = 0
    amp_meri_circ = 0.0
    rnoise_int = 3.9084999999999999E+37
    rnoise_ext = 3.9084999999999999E+37
    lreflecteddy = False
    louinit = False
    hydro_xaver_range = [-5.0000000000000003E-2,5.0000000000000003E-2]
    max_uu = 0.0
    amp_factor = 0.0
    kx_uu_perturb = 0.0
    llinearized_hydro = False
    hydro_zaver_range = [-1.7976931348623158E+307,1.7976931348623158E+307]
    index_rsh = 1
    ll_sh = numpy.repeat(0,5)
    mm_sh = numpy.repeat(0,5)
    ampllnrho = numpy.repeat(0.0,5)
    initlnrho = numpy.array(['const_rho']+list(numpy.repeat('nothing',4)))
    widthlnrho = numpy.repeat(0.10000000000000001,5)
    rho_left = numpy.repeat(1.0,5)
    rho_right = numpy.repeat(1.0,5)
    lnrho_const = 0.0
    hrho = 1.0
    rho_const = 1.0
    cs2bot = 1.0
    cs2top = 1.0
    radius_lnrho = numpy.repeat(0.5,5)
    eps_planet = 0.5
    xblob = 0.0
    yblob = 0.0
    zblob = 0.0
    b_ell = 1.0
    q_ell = 5.0
    hh0 = 0.0
    rbound = 1.0
    lwrite_stratification = False
    mpoly = 1.5
    ggamma = 1.6666666666666665
    strati_type = 'lnrho_ss'
    beta_glnrho_global = numpy.array([-0.10000000000000001]+list(numpy.repeat(0.0,2)))
    kx_lnrho = numpy.repeat(1.0,5)
    ky_lnrho = numpy.repeat(1.0,5)
    kz_lnrho = numpy.repeat(1.0,5)
    amplrho = numpy.repeat(0.0,5)
    phase_lnrho = numpy.repeat(0.0,5)
    coeflnrho = complex(0.0,0.0)
    kxx_lnrho = numpy.repeat(0.0,5)
    kyy_lnrho = numpy.repeat(0.0,5)
    kzz_lnrho = numpy.repeat(0.0,5)
    co1_ss = 0.0
    co2_ss = 0.0
    sigma1 = 150.0
    idiff = numpy.repeat('',4)
    ldensity_nolog = True
    wdamp = 0.0
    lcontinuity_gas = True
    lisothermal_fixed_hrho = False
    density_floor = -1.0
    lanti_shockdiffusion = False
    lmassdiff_fixmom = False
    lmassdiff_fixkin = False
    lrho_as_aux = False
    ldiffusion_nolog = False
    lnrho_z_shift = 0.0
    powerlr = 3.0
    zoverh = 1.5
    hoverr = 5.0000000000000003E-2
    lffree = False
    ffree_profile = 'none'
    rzero_ffree = 0.0
    wffree = 0.0
    rho_top = 1.0
    rho_bottom = 1.0
    r0_rho = 3.9084999999999999E+37
    invgrav_ampl = 0.0
    datafile = 'dens_temp.dat'
    mass_cloud = 0.0
    t_cloud = 0.0
    cloud_mode = 'isothermal'
    t_cloud_out_rel = 1.0
    xi_coeff = 1.0
    density_xaver_range = [-5.0000000000000003E-2,5.0000000000000003E-2]
    dens_coeff = 1.0
    temp_coeff = 1.0
    temp_trans = 0.0
    temp_coeff_out = 1.0
    reduce_cs2 = 1.0
    lreduced_sound_speed = False
    lscale_to_cs2top = False
    density_zaver_range = [-1.7976931348623158E+307,1.7976931348623158E+307]
    lconserve_total_mass = False
    total_mass = -1.0
    ireference_state = 'nothing'
    qshear = 1.5
    qshear0 = 0.0
    sshear = -1.5
    sshear1 = -1.5
    deltay = 0.0
    u0_advec = numpy.repeat(0.0,3)
    lshearadvection_as_shift = False
    shear_method = 'fft'
    lrandomx0 = False
    x0_shear = 0.0
    norder_poly = 3
    ltvd_advection = False
    lposdef_advection = False
    lposdef = numpy.repeat(False,5)
    lmagnetic_stretching = True
    sini = 0.0
    initxxp = numpy.array(['random']+list(numpy.repeat('nothing',4)))
    initvvp = numpy.array(['dragforce_equilibrium']+list(numpy.repeat('nothing',4)))
    xp0 = 0.0
    yp0 = 0.0
    zp0 = 0.0
    vpx0 = 0.0
    vpy0 = 0.0
    vpz0 = 0.0
    delta_vp0 = 1.0
    ldragforce_gas_par = False
    ldragforce_dust_par = False
    bcpx = 'p'
    bcpy = 'p'
    bcpz = 'p'
    tausp = 0.10000000000000001
    beta_dpdr_dust = 0.0
    np_swarm = 0.0
    mp_swarm = 3.9062500000000009E-6
    mpmat = 0.0
    rhop_swarm = 0.10000000000000001
    eps_dtog = 1.0
    nu_epicycle = 0.0
    rp_int = -3.9084999999999999E+37
    rp_ext = 3.9084999999999999E+37
    gravx_profile = ''
    gravz_profile = ''
    gravr_profile = ''
    gravx = 0.0
    gravz = 0.0
    gravr = 1.0
    gravsmooth = 0.0
    kx_gg = 1.0
    kz_gg = 1.0
    ri0 = 0.25
    eps1 = 0.5
    lmigration_redo = False
    ldragforce_equi_global_eps = False
    coeff = numpy.repeat(complex(0.0,0.0),7)
    kx_vvp = 0.0
    ky_vvp = 0.0
    kz_vvp = 0.0
    amplvvp = 0.0
    kx_xxp = 0.0
    ky_xxp = 0.0
    kz_xxp = 0.0
    amplxxp = 0.0
    kx_vpx = 0.0
    kx_vpy = 0.0
    kx_vpz = 0.0
    ky_vpx = 0.0
    ky_vpy = 0.0
    ky_vpz = 0.0
    kz_vpx = 0.0
    kz_vpy = 0.0
    kz_vpz = 0.0
    phase_vpx = 0.0
    phase_vpy = 0.0
    phase_vpz = 0.0
    lcoldstart_amplitude_correction = False
    particle_mesh = 'tsc'
    lparticlemesh_cic = False
    lparticlemesh_tsc = True
    linterpolate_spline = True
    tstart_dragforce_par = 0.0
    tstart_grav_par = 0.0
    lparticle_gravity = True
    tstart_grav_x_par = 0.0
    tstart_grav_z_par = 0.0
    tstart_grav_r_par = 0.0
    taucool = 0.0
    lcollisional_cooling_taucool = False
    lcollisional_cooling_rms = False
    lcollisional_cooling_twobody = False
    tausp_species = 0.10000000000000001
    tau_coll_min = 0.0
    ltau_coll_min_courant = True
    coeff_restitution = 0.5
    tstart_collisional_cooling = 0.0
    tausg_min = 0.0
    l_hole = 0
    m_hole = 0
    n_hole = 0
    epsp_friction_increase = 100.0
    lcollisional_dragforce_cooling = False
    ldragforce_heat = False
    lcollisional_heat = False
    lcompensate_friction_increase = False
    lmigration_real_check = True
    ldraglaw_epstein = True
    ldraglaw_simple = False
    ldraglaw_epstein_stokes_linear = False
    mean_free_path_gas = 0.0
    ldraglaw_epstein_transonic = False
    lcheck_exact_frontier = False
    ldraglaw_eps_stk_transonic = False
    dustdensity_powerlaw = 0.0
    rad_sphere = 0.0
    pos_sphere = numpy.repeat(0.0,3)
    ldragforce_stiff = False
    a_ellipsoid = 0.0
    b_ellipsoid = 0.0
    c_ellipsoid = 0.0
    pos_ellipsoid = numpy.repeat(0.0,3)
    ldraglaw_steadystate = False
    tstart_liftforce_par = 0.0
    ldraglaw_purestokes = False
    rpbeta_species = 0.0
    rpbeta = 0.0
    gab_width = 3.0
    tstart_brownian_par = 0.0
    tstart_sink_par = 0.0
    lbrownian_forces = False
    lthermophoretic_forces = False
    lenforce_policy = False
    interp_pol_uu = 'cic'
    interp_pol_oo = 'ngp'
    interp_pol_tt = 'ngp'
    interp_pol_rho = 'ngp'
    interp_pol_pp = 'ngp'
    interp_pol_species = 'ngp'
    brownian_t0 = 0.0
    thermophoretic_t0 = 0.0
    lnostore_uu = True
    ldt_grav_par = True
    ldragforce_radialonly = False
    lsinkpoint = False
    xsinkpoint = 0.0
    ysinkpoint = 0.0
    zsinkpoint = 0.0
    rsinkpoint = 0.0
    lcoriolis_force_par = True
    lcentrifugal_force_par = False
    ldt_adv_par = True
    lx0 = 0.0
    ly0 = 0.0
    lz0 = 0.0
    lglobalrandom = False
    linsert_particles_continuously = False
    lrandom_particle_pencils = False
    lnocalc_np = False
    lnocalc_rhop = False
    np_const = 0.0
    rhop_const = 0.0
    particle_radius = 0.0
    lignore_rhop_swarm = False
    ldragforce_equi_noback = False
    rhopmat = 1.0
    deltauy_gas_friction = 0.0
    xp1 = 0.0
    yp1 = 0.0
    zp1 = 0.0
    vpx1 = 0.0
    vpy1 = 0.0
    vpz1 = 0.0
    xp2 = 0.0
    yp2 = 0.0
    zp2 = 0.0
    vpx2 = 0.0
    vpy2 = 0.0
    vpz2 = 0.0
    xp3 = 0.0
    yp3 = 0.0
    zp3 = 0.0
    vpx3 = 0.0
    vpy3 = 0.0
    vpz3 = 0.0
    lsinkparticle_1 = False
    rsinkparticle_1 = 0.0
    lcalc_uup = False
    temp_grad0 = numpy.repeat(0.0,3)
    thermophoretic_eq = 'nothing'
    cond_ratio = 0.0
    interp_pol_gradtt = 'ngp'
    lreassign_strat_rhom = False
    lparticlemesh_pqs_assignment = False
    lwithhold_init_particles = False
    frac_init_particles = 1.0
    lvector_gravity = False
    birthring_r = 1.0
    birthring_width = 0.10000000000000001
    lgaussian_birthring = False
    ldraglaw_stokesschiller = False
    lbirthring_depletion = False
    remove_particle_at_time = -1.0
    remove_particle_criteria = 'all'
    remove_particle_criteria_size = 0.0
    lhydro = True
    ldensity = True
    lentropy = False
    ltemperature = False
    lshock = True
    lmagnetic = False
    lforcing = False
    llorenz_gauge = False
    ldustvelocity = False
    ldustdensity = False
    ltestscalar = False
    ltestfield = False
    ltestflow = False
    linterstellar = False
    lcosmicray = False
    lcosmicrayflux = False
    lshear = True
    lpscalar = False
    lsupersat = False
    lradiation = False
    leos = True
    lchiral = False
    lneutralvelocity = False
    lneutraldensity = False
    lpolymer = False
    lpointmasses = False
    lsolid_cells = False
    lpower_spectrum = True
    lparticles = True
    lparticles_drag = False
    lcollective_io = False
    io_strategy = 'dist'

This is the most convinient way! Same one could achieve by using the simulation object constructor directly. Beware, quiet=False is default for the constructor.

In [4]:
SIM = pcn.sim.simulation('../sample_simulations/2d_streaming_instability', quiet=True)

Lets make a copy of this simulation and modify it a bit. This needs two things: a) where to put the new simulation (path_root) and what is its name, which is its simulation dir as well!

In [5]:
SIM_copy = SIM.copy('../tmp/0002', 'copy_of_2d_si')

Now be careful! We stated name and path for the new simulation seperatly to the copy constructor. The reasons are:

  • Simulations should always rest in a folder which is the simulation name! This is also where simulation names are gathered from.
In [6]:
print(SIM_copy.name)
print(SIM_copy.path)

copy_of_2d_si
/home/user/pencil-code/python/tutorials/tmp/0002/copy_of_2d_si
  • If you want to create many simulations from one main simulation you should be able to do so in a few rows. We will do that below in an example.
In [7]:
root_path = '../tmp/0002'
for name in ['a', 'b', 'c']:
    SIM.copy(root_path, name)
os.listdir(root_path)
Out[7]:
['b', 'copy_of_2d_si', 'c', 'a', '2d_streaming_instability_copy']
  • One should be able to make a copy of a simulation without knowing its name!
In [8]:
SIM.copy(root_path)
os.listdir(root_path)

? Warning: No name specified, will alter old simulation name to 2d_streaming_instability_copy
Out[8]:
['b', 'copy_of_2d_si', 'c', 'a', '2d_streaming_instability_copy']

Changing simulation parameters

So far, we havnt implemented a way to modify simulation modules from python! But you can modify 'src/cparam.local' and 'start.in' and 'run.in'. So, also ways to modify 'print.in', 'video.in', etc. are missing. Feel free to add them!

The reason they are missing is, that valid input parameters depend on modules and its sometimes unclear where and how to add parameters correcly. But, changing them is rather easy and anyway what we do most of the time!

So, you workflow should be:

  1. Find a suiting master simulation or set up one on your own.
  2. Frome this master produce your simulation copys.
  3. Change parameters of the simulations.

You know already how to do the first two things. Lets explore how we can do the last step from python.

Do changes in cparam.local

Frist, lets get rid of the simulation we created beforehand. Therefor, we use a the ability to get all simulations that are within a folder! This is very mighty, as you will see in the following. Since deleting should always considered seriously, we first do a dry removal run.

In [9]:
SIMS = pcn.get_sims(root_path)
for SIM in SIMS:
    print('! Removing '+SIM.name)
    #SIM.remove()

! Removing b
! Removing copy_of_2d_si
! Removing c
! Removing a
! Removing 2d_streaming_instability_copy

Do changes in start.in and run.in

In [10]:
SIM_copy.compile()

! ERROR: Compiling ended with error code 2!
! Please check log file in
! /home/user/pencil-code/python/tutorials/tmp/0002/copy_of_2d_si/pc/compiler_log_2016-12-22_16:36:35.032424
Out[10]:
2

Compile and run simulations

In [ ]: