Commit 5650eb51 authored by adelmann's avatar adelmann 🎗
Browse files

Merge branch 'master' of gitlab.psi.ch:OPAL/src

parents 3876559f 05ec58ef
......@@ -12,14 +12,16 @@ size_t AmrPartBunch::nAttributes = 10;
AmrPartBunch::AmrPartBunch(const Geometry& geom,
const DistributionMapping & dmap,
const BoxArray & ba)
: ParticleContainer(geom, dmap, ba)
: ParticleContainer(geom, dmap, ba),
nLocalParticles_m(0)
{ }
AmrPartBunch::AmrPartBunch(const Array<Geometry>& geom,
const Array<DistributionMapping>& dmap,
const Array<BoxArray>& ba,
const Array<int>& rr)
: ParticleContainer(geom, dmap, ba, rr)
: ParticleContainer(geom, dmap, ba, rr),
nLocalParticles_m(0)
{ }
......@@ -248,4 +250,4 @@ void AmrPartBunch::setB(Vector_t Bf, int i) {
for (int d = 0; d < 3; ++d)
m_particles[l][g][dq].m_data[d + 7] = Bf(d);
}
\ No newline at end of file
}
......@@ -24,10 +24,18 @@ void Distribution::uniform(double lower, double upper, size_t nloc, int seed) {
nloc_m = nloc;
std::mt19937_64 mt(seed);
std::mt19937_64 mt(0/*seed*/ /*0*/);
std::uniform_real_distribution<> dist(lower, upper);
// // assume that seed == rank of node
// mt.discard(6 * (nloc + 1) * seed);
// assume that seed == rank of node
// inefficient but only way to make sure that parallel distribution is equal to sequential
for (size_t i = 0; i < 6 * nloc_m * seed; ++i)
dist(mt);
x_m.resize(nloc);
y_m.resize(nloc);
z_m.resize(nloc);
......@@ -55,10 +63,18 @@ void Distribution::gaussian(double mean, double stddev, size_t nloc, int seed) {
nloc_m = nloc;
std::mt19937_64 mt(seed);
std::mt19937_64 mt(0/*seed*/ /*0*/);
std::normal_distribution<double> dist(mean, stddev);
// // assume that seed == rank of node
// mt.discard(6 * (nloc + 1) * seed);
// assume that seed == rank of node
// inefficient but only way to make sure that parallel distribution is equal to sequential
for (size_t i = 0; i < 6 * nloc_m * seed; ++i)
dist(mt);
x_m.resize(nloc);
y_m.resize(nloc);
z_m.resize(nloc);
......
#include "Solver.h"
#include "Ippl.h"
void
Solver::solve_for_accel(container_t& rhs,
container_t& phi,
......@@ -9,7 +11,8 @@ Solver::solve_for_accel(container_t& rhs,
int finest_level,
Real offset)
{
static IpplTimings::TimerRef edge2centerTimer = IpplTimings::getTimer("grad-edge2center");
Real tol = 1.e-10;
Real abs_tol = 1.e-14;
......@@ -48,6 +51,7 @@ Solver::solve_for_accel(container_t& rhs,
abs_tol);
// Average edge-centered gradients to cell centers and fill the values in ghost cells.
IpplTimings::startTimer(edge2centerTimer);
for (int lev = base_level; lev <= finest_level; lev++)
{
#ifdef UNIQUE_PTR
......@@ -79,6 +83,7 @@ Solver::solve_for_accel(container_t& rhs,
grad_phi[lev].mult(-1.0, 0, 3);
}
#endif
IpplTimings::stopTimer(edge2centerTimer);
}
......@@ -92,6 +97,12 @@ Solver::solve_with_f90(container_t& rhs,
Real tol,
Real abs_tol)
{
static IpplTimings::TimerRef initSolverTimer = IpplTimings::getTimer("init-solver");
static IpplTimings::TimerRef doSolveTimer = IpplTimings::getTimer("do-solve");
static IpplTimings::TimerRef gradientTimer = IpplTimings::getTimer("gradient");
IpplTimings::startTimer(initSolverTimer);
int nlevs = finest_level - base_level + 1;
int mg_bc[2*BL_SPACEDIM];
......@@ -99,7 +110,7 @@ Solver::solve_with_f90(container_t& rhs,
// This tells the solver that we are using Dirichlet bc's
if (Geometry::isAllPeriodic()) {
// if ( ParallelDescriptor::IOProcessor() )
std::cerr << "Periodic BC" << std::endl;
// std::cerr << "Periodic BC" << std::endl;
for (int dir = 0; dir < BL_SPACEDIM; ++dir) {
// periodic BC
......@@ -108,7 +119,7 @@ Solver::solve_with_f90(container_t& rhs,
}
} else {
// if ( ParallelDescriptor::IOProcessor() )
std::cerr << "Dirichlet BC" << std::endl;
// std::cerr << "Dirichlet BC" << std::endl;
for (int dir = 0; dir < BL_SPACEDIM; ++dir) {
// Dirichlet BC
......@@ -177,8 +188,14 @@ Solver::solve_with_f90(container_t& rhs,
int always_use_bnorm = 0;
int need_grad_phi = 1;
fmg.set_verbose(1);
IpplTimings::stopTimer(initSolverTimer);
IpplTimings::startTimer(doSolveTimer);
fmg.solve(phi_p, rhs_p, tol, abs_tol, always_use_bnorm, need_grad_phi);
IpplTimings::stopTimer(doSolveTimer);
IpplTimings::startTimer(gradientTimer);
#ifdef UNIQUE_PTR
const Array<Array<MultiFab*> >& g_phi_edge = BoxLib::GetArrOfArrOfPtrs(grad_phi_edge);
for (int ilev = 0; ilev < nlevs; ++ilev) {
......@@ -191,7 +208,7 @@ Solver::solve_with_f90(container_t& rhs,
fmg.get_fluxes(grad_phi_edge[amr_level], ilev);
}
#endif
IpplTimings::stopTimer(gradientTimer);
}
#ifdef USEHYPRE
......
##
# @file field.py
# @author Matthias Frey
# @date 4. Jan. 2017
# @version 1.0
# @brief Plot vector and scalar fields
# @details When configuring OPAL with the flag
# -DDBG_SCALARFIELD=1 it dumps the
# scalar and vector fields of every
# time step to "data" directory of the
# simulation. These are visualized by
# this script.\n
# Call: python field.py [file] \n
# where [file] either *field*
# or *scalar*
#
from tools import match, doGridPlot
import sys
import numpy as np
# regular expression pattern, where parentheses () represent groups that can be accessed later
# group(0) is whole match
# format: i j k ( x-component , y-component , z-component )
# where i, j and k are the grid points
vector_field_pattern = r'(\d*) (\d*) (\d*) \( (-*.*\d*) \, (-*.*\d*) \, (-*.*\d*) \)'
# format: i j k value
# where i, j and k are the grid points
scalar_field_pattern = r'(\d*) (\d*) (\d*) (-*.*\d*)'
try:
f = sys.argv[1]
if 'field' in f:
data = match(f, vector_field_pattern)
# electric field in horizontal direction
mid = int(0.5 * max(data[:, 2]))
i = np.extract(data[:, 2] == mid, data[:, 0])
j = np.extract(data[:, 2] == mid, data[:, 1])
val = np.extract(data[:, 2] == mid, data[:, 3])
plt = doGridPlot(i, j, val, "horizontal grid", "vertical grid", r'$E_x$ [V/m]')
plt.savefig(f + '_slice_z_Ex.png')#, bbox_inches='tight')
print ("max. Ex: ", max(val))
print ("min. Ex: ", min(val))
# electric field in vertical direction
val = np.extract(data[:, 2] == mid, data[:, 4])
plt = doGridPlot(i, j, val, "horizontal grid", "vertical grid", r'$E_y$ [V/m]')
plt.savefig(f + '_slice_z_Ey.png')#bbox_inches='tight')
print ("max. Ey: ", max(val))
print ("min. Ey: ", min(val))
# center in z
mid = int(0.5 * max(data[:, 0]))
i = np.extract(data[:, 0] == mid, data[:, 1])
j = np.extract(data[:, 0] == mid, data[:, 2])
val = np.extract(data[:, 0] == mid, data[:, 5])
plt = doGridPlot(i, j, val, "vertical grid", "longitudinal grid", r'$E_z$ [V/m]')
plt.savefig(f + '_slice_x_Ez.png')# bbox_inches='tight')
print ("max. Ez: ", max(val))
print ("min. Ez: ", min(val))
elif 'scalar' in f:
data = match(f, scalar_field_pattern)
clab = ''
if 'phi' in f:
clab = r'$\phi$ [V]'
elif 'rho' in f:
clab = r'$\rho$ [C/m**3]'
# center in x
mid = int(0.5 * max(data[:, 0]))
i = np.extract(data[:, 0] == mid, data[:, 1])
j = np.extract(data[:, 0] == mid, data[:, 2])
val = np.extract(data[:, 0] == mid, data[:, 3])
plt = doGridPlot(i, j, val, "vertical grid", "longitudinal grid", clab)
plt.savefig(f + '_slice_x.png')#, bbox_inches='tight')
# center in y
mid = int(0.5 * max(data[:, 1]))
i = np.extract(data[:, 1] == mid, data[:, 0])
j = np.extract(data[:, 1] == mid, data[:, 2])
val = np.extract(data[:, 1] == mid, data[:, 3])
plt = doGridPlot(i, j, val, "horizontal grid", "longitudinal grid", clab)
plt.savefig(f + '_slice_y.png')#bbox_inches='tight')
# center in z
mid = int(0.5 * max(data[:, 2]))
i = np.extract(data[:, 2] == mid, data[:, 0])
j = np.extract(data[:, 2] == mid, data[:, 1])
val = np.extract(data[:, 2] == mid, data[:, 3])
plt = doGridPlot(i, j, val, "horizontal grid", "vertical grid", clab)
plt.savefig(f + '_slice_z.png')# bbox_inches='tight')
print ("max. value (slice z): ", max(val))
print ("min. value (slice z): ", min(val))
else:
raise RuntimeError('This type of file is not supported.')
except Exception as e:
print (e.strerror)
\ No newline at end of file
##
# @file tools.py
# @author Matthias Frey
# @date 22. Dec. 2016
# @version 1.0
# @brief This module contains several helper functions
# for data analysis
#
import yt
import os
import re
import numpy as np
import matplotlib.pyplot as plt
##
# @param ds is the data
# @param direct is the direction 'x', 'y' or 'z' (normal)
# @param field to plot
# @param unit the data should be converted to (otherwise it
# takes the default given by the data)
# @param col is the color for the time stamp and scale annotation
def doSlicePlot(ds, direct, field, unit, col = 'white'):
slc = yt.SlicePlot(ds, normal=direct, fields=field)
if unit is not None:
slc.set_unit(field, unit)
slc.annotate_grids()
slc.annotate_timestamp(corner='upper_left', redshift=False, draw_inset_box=True)
slc.annotate_scale(corner='upper_right', size_bar_args={'color':col})
slc.save()
##
# @param ds is the data
# @param direct is the direction 'x', 'y' or 'z' (normal)
# @param field to plot
# @param unit the data should be converted to (otherwise it
# takes the default given by the data)
# @param col is the color for the time stamp and scale annotation
def doProjectionPlot(ds, direct, field, unit, col = 'white'):
slc = yt.ProjectionPlot(ds, direct, fields=field)
if unit is not None:
slc.set_unit(field, unit)
slc.annotate_grids()
slc.annotate_timestamp(corner='upper_left', redshift=False, draw_inset_box=True)
slc.annotate_scale(corner='upper_right', size_bar_args={'color':col})
slc.save()
##
# Take an integer and transform it
# to a string of four characters, e.g.\n
# \f$ 1 \rightarrow 0001 \f$ \n
# \f$ 12 \rightarrow 0012 \f$ \n
# \f$ 586 \rightarrow 0586 \f$ \n
# @param step is an integer
def concatenate(step):
res = str(step)
while len(res) < 4:
res = '0' + res
return res
##
# Count subdirectories
# @param parent is the path to the parent directory
# @param substr specifies the substring that should be contained
# @returns the number of subdirectories containing
# a given substring
def countSubdirs(parent, substr):
nDirs = 0
# 22. Dec. 2016
# http://stackoverflow.com/questions/10377998/how-can-i-iterate-over-files-in-a-given-directory
for filename in os.listdir(parent):
if substr in filename:
nDirs = nDirs + 1
return nDirs
##
# Read field data written by OPAL (-DDBG_SCALARFIELD=1) using
# regular expression
# @param f is the file that is read line by line
# @param pattern is a string specifying rule for matching
# @returns a matrix where each row is matched line in the file
def match(f, pattern):
regex = re.compile(pattern)
data = np.empty((0, regex.groups))
# go through file line by line and match
with open(f) as ff:
for line in ff:
res = regex.match(line)
row = []
for i in range(1, regex.groups + 1):
row = np.append(row, [float(res.group(i))])
data = np.append(data, np.array([row]), axis=0)
return data
##
def doGridPlot(i, j, val, xlab, ylab, clab):
imin = int(i[0])
imax = int(i[-1])
jmin = int(j[0])
jmax = int(j[-1])
i = np.reshape(i, (imax, jmax))
j = np.reshape(j, (imax, jmax))
val = np.reshape(val, (imax, jmax))
plt.figure()
plt.pcolor(i, j, val, cmap='YlGnBu')
plt.xlim(imin, imax)
plt.xlabel(xlab)
plt.ylim(imin, jmax)
plt.ylabel(ylab)
plt.colorbar(label=clab)
return plt
##
# @file tracking_vis.py
# @author Matthias Frey
# @date 14. October 2016, LBNL
# @version 1.1 (22. Dec. 2016)
# @brief Plot the electric self-field, density and self-field
# potential using the yt framework.
# @details It reads in the all generated BoxLib plotfiles "plt_amr_*"
# (* are four digits) of a directory and saves the content as slice
# plots. It reads in the output of testReal.cpp.
# Make sure you sourced the yt directory $YT_DIR/yt-x86_64/bin/activate.
#
import os
import sys
import yt
from tools import countSubdirs, concatenate, doSlicePlot
try:
opal = os.environ['OPAL_BUILD']
# leading name of plotfiles
substr = 'plt_amr_'
nFiles = countSubdirs(opal + '/ippl/test/AMR/', substr)
print ('Found ' + str(nFiles) + ' plotfiles')
for i in range(0, nFiles):
# do a string of four characters, e.g. 12 --> 0012
res = concatenate(i)
ds = yt.load(opal + '/ippl/test/AMR/' + substr + res, unit_system='mks')
doSlicePlot(ds, 'z', 'rho', 'C/m**3', 'gray')
doSlicePlot(ds, 'z', 'potential', 'V')
doSlicePlot(ds, 'z', 'Ex', 'V/m')
doSlicePlot(ds, 'z', 'Ey', 'V/m')
doSlicePlot(ds, 'x', 'Ez', 'V/m')
except KeyError:
print ("Please export the environment variable 'OPAL_BUILD'.")
except IOError as e:
print (e.strerror)
except TypeError:
print ("No subdirectories containing the given substring " +
"found in " + opal + "'/ippl/test/AMR/'.")
\ No newline at end of file
##
# @file visualize.py
# @author Matthias Frey
# @date 14. October 2016, LBNL
# @version 1.1 (21. Dec. 2016)
#
# @pre Environment variable OPAL_BUILD has to be set.
# @details Plot the electric self-field, density and self-field
# potential using the yt framework.
# 1. mpirun -np #cores testSolver
# 2. python visualize.py (make sure you sourced the yt directory $YT_DIR/yt-x86_64/bin/activate)
# @brief Slice plots of plotfiles generated by writePlotFile.H
import os
import yt
from tools import doSlicePlot, doProjectionPlot
try:
opal = os.environ['OPAL_BUILD']
ds = yt.load(opal + "ippl/test/AMR/plt0000", dataset_type='opal')
ds.print_stats()
print ("Field list:", ds.field_list)
print ("Derived field list:", ds.derived_field_list)
doSlicePlot(ds, 'z', 'rho', 'C/m**3', 'gray')
doSlicePlot(ds, 'y', 'rho', 'C/m**3', 'gray')
doSlicePlot(ds, 'x', 'rho', 'C/m**3', 'gray')
doProjectionPlot(ds, 'x', 'rho', 'C/m**2', 'gray')
doProjectionPlot(ds, 'y', 'rho', 'C/m**2', 'gray')
doProjectionPlot(ds, 'z', 'rho', 'C/m**2', 'gray')
doSlicePlot(ds, 'z', 'Ex', 'V/m')
doSlicePlot(ds, 'z', 'Ey', 'V/m')
doSlicePlot(ds, 'x', 'Ez', 'V/m')
doSlicePlot(ds, 'z', 'potential', 'V')
doSlicePlot(ds, 'y', 'potential', 'V')
doSlicePlot(ds, 'x', 'potential', 'V')
ad = ds.all_data()
print ( ad.quantities.extrema("rho").in_units('C/m**3') )
print ( ad.quantities.extrema("Ex").in_units('V/m') )
print ( ad.quantities.extrema("potential").in_units('V') )
except KeyError:
print ("Please export the environment variable 'OPAL_BUILD'.")
except IOError as e:
print (e.strerror)
##
# @file analyse.py
# @author Matthias Frey
# @date 26. Oct. 2016
# @details The files generated by error.cpp can be visualized using this script.
# It plots the number of levels vs. the Euclidean error norm where the
# error is computed between the single- and multi-level solution.
# @pre Environment variable OPAL_BUILD has to be set.
# @brief Plot the Euclidean error norm vs. the level of refinement
import numpy as np
import os
import matplotlib.pyplot as plt
## Absolute path to OPAL build directory
opal = os.environ['OPAL_BUILD']
##
# Create an error plot
# @param filename is the absolute path
# @param xlab is the label on the x-axis
# @param ylab is the label on the y-axis
def doPlot(filename, xlab, ylab):
# error with respect to single-level solve.
nLevels, l2err = np.loadtxt(filename, unpack=True)
plt.figure()
plt.plot(nLevels, l2err, "-o")
plt.xlabel(xlab)
plt.ylabel(ylab)
plt.show()
doPlot(opal + "/ippl/test/AMR/l2_error_NOPARTICLES.dat",
'#levels', r'$L_{2}$-error')
doPlot(opal + "/ippl/test/AMR/l2_error_UNIFORM.dat",
'#levels', r'$L_{2}$-error')
doPlot(opal + "/ippl/test/AMR/l2_error_GAUSSIAN.dat",
'#levels', r'$L_{2}$-error')
......@@ -293,9 +293,9 @@ inline double totalCharge(const container_t& rhs,
for (int i = 0; i <= finest_level; ++i) {
Real vol = (*(geom[i].CellSize()) * *(geom[i].CellSize()) * *(geom[i].CellSize()) );
#ifdef UNIQUE_PTR
Real sum = rhs[i]->sum(0) * vol * scale;
Real sum = (scale) ? rhs[i]->sum(0) * vol : rhs[i]->sum(0);
#else
Real sum = rhs[i].sum(0) * vol * scale;
Real sum = (scale) ? rhs[i].sum(0) * vol : rhs[i].sum(0);
#endif
totCharge += sum;
}
......
......@@ -53,25 +53,33 @@ void doSolve(AmrOpal& myAmrOpal, PartBunchBase* bunch,
int nLevels,
Inform& msg)
{
static IpplTimings::TimerRef solvTimer = IpplTimings::getTimer("solv");
static IpplTimings::TimerRef allocTimer = IpplTimings::getTimer("alloc-memory-grid");
static IpplTimings::TimerRef assignTimer = IpplTimings::getTimer("assign-charge");
// =======================================================================
// 4. prepare for multi-level solve
// =======================================================================
rhs.resize(nLevels);
phi.resize(nLevels);
grad_phi.resize(nLevels);
IpplTimings::startTimer(allocTimer);
for (int lev = 0; lev < nLevels; ++lev) {
initGridData(rhs, phi, grad_phi, myAmrOpal.boxArray()[lev], lev);
}
IpplTimings::stopTimer(allocTimer);
// Define the density on level 0 from all particles at all levels
int base_level = 0;
int finest_level = myAmrOpal.finestLevel();
IpplTimings::startTimer(assignTimer);
dynamic_cast<AmrPartBunch*>(bunch)->AssignDensity(0, false, rhs, base_level, 1, finest_level);
IpplTimings::stopTimer(assignTimer);
double totCharge = totalCharge(rhs, finest_level, geom);
msg << "Total Charge: " << totCharge << " C" << endl
......@@ -98,7 +106,6 @@ void doSolve(AmrOpal& myAmrOpal, PartBunchBase* bunch,
// solve
Solver sol;
IpplTimings::startTimer(solvTimer);
sol.solve_for_accel(rhs,
phi,
grad_phi,