Skip to content

Commit 9c7bce8

Browse files
rochSmetsrochSmets
committed
introduce charge, mass and particle densities
1 parent ce77d6b commit 9c7bce8

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

46 files changed

+717
-268
lines changed

‎.gdbinit

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
break __assert_fail
2+
break abort
3+
catch throw

‎pyphare/pyphare/pharein/diagnostics.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -222,6 +222,7 @@ def population_in_model(population):
222222
class FluidDiagnostics_(Diagnostics):
223223
fluid_quantities = [
224224
"density",
225+
"charge_density",
225226
"mass_density",
226227
"flux",
227228
"bulkVelocity",

‎pyphare/pyphare/pharesee/hierarchy/__init__.py

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,10 +11,11 @@
1111

1212

1313
def hierarchy_from(
14-
simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, **kwargs
14+
simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, func=None, **kwargs
1515
):
1616
from .fromh5 import hierarchy_fromh5
1717
from .fromsim import hierarchy_from_sim
18+
from .fromfunc import hierarchy_from_func
1819

1920
"""
2021
this function reads an HDF5 PHARE file and returns a PatchHierarchy from
@@ -39,4 +40,7 @@ def hierarchy_from(
3940
if simulator is not None and qty is not None:
4041
return hierarchy_from_sim(simulator, qty, pop=pop)
4142

43+
if func is not None and hier is not None:
44+
return hierarchy_from_func(func, hier, **kwargs)
45+
4246
raise ValueError("can't make hierarchy")

‎pyphare/pyphare/pharesee/hierarchy/fromfunc.py

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from
2+
3+
import numpy as np
4+
5+
6+
def ions_mass_density_func1d(x, **kwargs):
7+
masses = kwargs["masses"] # list of float : the ion pop masses
8+
densities = kwargs["densities"] # list of callable : the ion pop density profiles
9+
10+
assert len(masses) == len(densities)
11+
funcs = np.zeros((x.size, len(masses)))
12+
13+
for i, (mass, density) in enumerate(zip(masses, densities)):
14+
funcs[:,i] = mass*density(x)
15+
16+
return funcs.sum(axis=1)
17+
18+
19+
def ions_charge_density_func1d(x, **kwargs):
20+
charges = kwargs["charges"] # list of float : the ion pop charges
21+
densities = kwargs["densities"] # list of callable : the ion pop density profiles
22+
23+
assert len(charges) == len(densities)
24+
25+
funcs = np.zeros((x.size, len(charges)))
26+
27+
for i, (charge, density) in enumerate(zip(charges, densities)):
28+
funcs[:,i] = charge*density(x)
29+
30+
return funcs.sum(axis=1)
31+
32+
33+
def hierarchy_from_func1d(func, hier, **kwargs):
34+
assert hier.ndim == 1
35+
36+
def compute_(patch_datas, **kwargs):
37+
ref_name = next(iter(patch_datas.keys()))
38+
x_ = patch_datas[ref_name].x
39+
40+
return (
41+
{"name": "value", "data": func(x_, **kwargs), "centering": patch_datas[ref_name].centerings},
42+
)
43+
44+
return compute_hier_from(compute_, hier, **kwargs)
45+
46+
47+
def ions_mass_density_func2d(x, y, **kwargs):
48+
masses = kwargs["masses"] # list of float : the ion pop masses
49+
densities = kwargs["densities"] # list of callable : the ion pop density profiles
50+
51+
yv, xv = np.meshgrid(y, x)
52+
53+
assert len(masses) == len(densities)
54+
funcs = np.zeros((x.size, y.size, len(masses)))
55+
56+
for i, (mass, density) in enumerate(zip(masses, densities)):
57+
funcs[:,:,i] = mass*density(xv, yv)
58+
59+
return funcs.sum(axis=2)
60+
61+
62+
def ions_charge_density_func2d(x, y, **kwargs):
63+
charges = kwargs["charges"] # list of float : the ion pop charges
64+
densities = kwargs["densities"] # list of callable : the ion pop density profiles
65+
66+
yv, xv = np.meshgrid(y, x)
67+
68+
assert len(charges) == len(densities)
69+
funcs = np.zeros((x.size, y.size, len(charges)))
70+
71+
for i, (charge, density) in enumerate(zip(charges, densities)):
72+
funcs[:,:,i] = charge*density(xv, yv)
73+
74+
return funcs.sum(axis=2)
75+
76+
77+
def hierarchy_from_func2d(func, hier, **kwargs):
78+
assert hier.ndim == 2
79+
80+
def compute_(patch_datas, **kwargs):
81+
ref_name = next(iter(patch_datas.keys()))
82+
x_ = patch_datas[ref_name].x
83+
y_ = patch_datas[ref_name].y
84+
85+
return (
86+
{"name": "value", "data": func(x_, y_, **kwargs), "centering": patch_datas[ref_name].centerings},
87+
)
88+
89+
return compute_hier_from(compute_, hier, **kwargs)
90+
91+
92+
def hierarchy_from_func(func, hier, **kwargs):
93+
if hier.ndim == 1:
94+
return hierarchy_from_func1d(func, hier, **kwargs)
95+
if hier.ndim == 2:
96+
return hierarchy_from_func2d(func, hier, **kwargs)
97+

‎pyphare/pyphare/pharesee/hierarchy/hierarchy.py

Lines changed: 39 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -360,6 +360,8 @@ def box_to_Rectangle(self, box):
360360
return Rectangle(box.lower, *box.shape)
361361

362362
def plot_2d_patches(self, ilvl, collections, **kwargs):
363+
from matplotlib.patches import Rectangle
364+
363365
if isinstance(collections, list) and all(
364366
[isinstance(el, Box) for el in collections]
365367
):
@@ -370,35 +372,47 @@ def plot_2d_patches(self, ilvl, collections, **kwargs):
370372
level_domain_box = self.level_domain_box(ilvl)
371373
mi, ma = level_domain_box.lower.min(), level_domain_box.upper.max()
372374

373-
fig, ax = kwargs.get("subplot", plt.subplots(figsize=(6, 6)))
375+
fig, ax = kwargs.get("subplot", plt.subplots(figsize=(16, 16)))
376+
377+
i0, j0 = level_domain_box.lower
378+
i1, j1 = level_domain_box.upper
379+
ij = np.zeros((i1 - i0 + 1, j1 - j0 + 1)) + np.nan
380+
ix = np.arange(i0, i1 + 1)
381+
iy = np.arange(j0, j1 + 1)
374382

375383
for collection in collections:
376-
facecolor = collection.get("facecolor", "none")
377-
edgecolor = collection.get("edgecolor", "purple")
378-
alpha = collection.get("alpha", 1)
379-
rects = [self.box_to_Rectangle(box) for box in collection["boxes"]]
380-
381-
ax.add_collection(
382-
PatchCollection(
383-
rects, facecolor=facecolor, alpha=alpha, edgecolor=edgecolor
384-
)
385-
)
384+
value = collection.get("value", np.nan)
385+
for box in collection["boxes"]:
386+
i0, j0 = box.lower
387+
i1, j1 = box.upper
388+
ij[i0 : i1 + 1, j0 : j1 + 1] = value
389+
if "coords" in collection:
390+
for coords in collection["coords"]:
391+
ij[coords] = collection["value"]
392+
393+
ax.pcolormesh(ix, iy, ij.T, edgecolors="k", cmap="jet")
394+
ax.set_xticks(ix)
395+
ax.set_yticks(iy)
396+
397+
for patch in self.level(ilvl).patches:
398+
box = patch.box
399+
r = Rectangle(box.lower - 0.5, *(box.upper + 0.5))
400+
401+
r.set_edgecolor("r")
402+
r.set_facecolor("none")
403+
r.set_linewidth(2)
404+
ax.add_patch(r)
386405

387406
if "title" in kwargs:
388407
from textwrap import wrap
389408

390409
xfigsize = int(fig.get_size_inches()[0] * 10) # 10 characters per inch
391410
ax.set_title("\n".join(wrap(kwargs["title"], xfigsize)))
392411

393-
major_ticks = np.arange(mi - 5, ma + 5 + 5, 5)
394-
ax.set_xticks(major_ticks)
395-
ax.set_yticks(major_ticks)
396-
397-
minor_ticks = np.arange(mi - 5, ma + 5 + 5, 1)
398-
ax.set_xticks(minor_ticks, minor=True)
399-
ax.set_yticks(minor_ticks, minor=True)
400-
401-
ax.grid(which="both")
412+
if "xlim" in kwargs:
413+
ax.set_xlim(kwargs["xlim"])
414+
if "ylim" in kwargs:
415+
ax.set_ylim(kwargs["ylim"])
402416

403417
return fig
404418

@@ -433,7 +447,7 @@ def plot1d(self, **kwargs):
433447
qty = pdata_names[0]
434448

435449
layout = patch.patch_datas[qty].layout
436-
nbrGhosts = layout.nbrGhostFor(qty)
450+
nbrGhosts = patch.patch_datas[qty].ghosts_nbr
437451
val = patch.patch_datas[qty][patch.box]
438452
x = patch.patch_datas[qty].x[nbrGhosts[0] : -nbrGhosts[0]]
439453
label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip)
@@ -544,9 +558,10 @@ def plot2d(self, **kwargs):
544558
if "ylim" in kwargs:
545559
ax.set_ylim(kwargs["ylim"])
546560

547-
divider = make_axes_locatable(ax)
548-
cax = divider.append_axes("right", size="5%", pad=0.08)
549-
plt.colorbar(im, ax=ax, cax=cax)
561+
if kwargs.get("cbar", True):
562+
divider = make_axes_locatable(ax)
563+
cax = divider.append_axes("right", size="5%", pad=0.08)
564+
fig.colorbar(im, ax=ax, cax=cax)
550565

551566
if kwargs.get("legend", None) is not None:
552567
ax.legend()

‎pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@
3636
"momentum_tensor_zz": "Mzz",
3737
"density": "rho",
3838
"mass_density": "rho",
39+
"charge_density": "rho",
3940
"tags": "tags",
4041
}
4142

‎pyphare/pyphare/pharesee/plotting.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -274,7 +274,7 @@ def finest_field_plot(run_path, qty, **kwargs):
274274
time = times[0]
275275
interpolator, finest_coords = r.GetVi(time, merged=True, interp=interp)[qty]
276276
elif qty == "rho":
277-
file = os.path.join(run_path, "ions_density.h5")
277+
file = os.path.join(run_path, "ions_charge_density.h5")
278278
if time is None:
279279
times = get_times_from_h5(file)
280280
time = times[0]

‎pyphare/pyphare/pharesee/run/run.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@
2727
"EM_B": "B",
2828
"EM_E": "E",
2929
"ions_bulkVelocity": "Vi",
30-
"ions_density": "Ni",
30+
"ions_charge_density": "Ni",
3131
"particle_count": "nppc",
3232
}
3333

@@ -114,7 +114,7 @@ def GetMassDensity(self, time, merged=False, interp="nearest", **kwargs):
114114
return ScalarField(self._get(hier, time, merged, interp))
115115

116116
def GetNi(self, time, merged=False, interp="nearest", **kwargs):
117-
hier = self._get_hierarchy(time, "ions_density.h5", **kwargs)
117+
hier = self._get_hierarchy(time, "ions_charge_density.h5", **kwargs)
118118
return ScalarField(self._get(hier, time, merged, interp))
119119

120120
def GetN(self, time, pop_name, merged=False, interp="nearest", **kwargs):
@@ -151,7 +151,7 @@ def GetPi(self, time, merged=False, interp="nearest", **kwargs):
151151
return self._get(Pi, time, merged, interp) # should later be a TensorField
152152

153153
def GetPe(self, time, merged=False, interp="nearest", all_primal=True):
154-
hier = self._get_hierarchy(time, "ions_density.h5")
154+
hier = self._get_hierarchy(time, "ions_charge_density.h5")
155155

156156
Te = hier.sim.electrons.closure.Te
157157

‎pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -177,11 +177,11 @@ def test_large_patchoverlaps(self, expected):
177177
collections=[
178178
{
179179
"boxes": overlap_boxes,
180-
"facecolor": "yellow",
180+
"value": 1,
181181
},
182182
{
183183
"boxes": [p.box for p in hierarchy.level(ilvl).patches],
184-
"facecolor": "grey",
184+
"value": 2,
185185
},
186186
],
187187
)
@@ -390,11 +390,11 @@ def test_level_ghost_boxes(self, refinement_boxes, expected):
390390
collections=[
391391
{
392392
"boxes": ghost_area_box_list,
393-
"facecolor": "yellow",
393+
"value": 1,
394394
},
395395
{
396396
"boxes": [p.box for p in hierarchy.level(ilvl).patches],
397-
"facecolor": "grey",
397+
"value": 2,
398398
},
399399
],
400400
title="".join(
@@ -502,7 +502,7 @@ def test_patch_periodicity_copy(self, refinement_boxes, expected):
502502
collections=[
503503
{
504504
"boxes": [p.box for p in periodic_list],
505-
"facecolor": "grey",
505+
"value": 2,
506506
},
507507
],
508508
)

‎pyphare/pyphare_tests/test_pharesee/test_hierarchy.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -145,7 +145,7 @@ def vthz(x, y):
145145
for quantity in ["E", "B"]:
146146
ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
147147

148-
for quantity in ["density", "bulkVelocity"]:
148+
for quantity in ["charge_density", "bulkVelocity"]:
149149
ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
150150

151151
return sim

‎src/amr/level_initializer/hybrid_level_initializer.hpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -98,7 +98,7 @@ namespace solver
9898
core::depositParticles(ions, layout, interpolate_, core::LevelGhostDeposit{});
9999
}
100100

101-
ions.computeDensity();
101+
ions.computeChargeDensity();
102102
ions.computeBulkVelocity();
103103
}
104104

‎src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -385,20 +385,21 @@ namespace amr
385385
{
386386
// first thing to do is to project patchGhostParitcles moments
387387
auto& patchGhosts = pop.patchGhostParticles();
388-
auto& density = pop.density();
388+
auto& particleDensity = pop.particleDensity();
389+
auto& chargeDensity = pop.chargeDensity();
389390
auto& flux = pop.flux();
390391

391-
interpolate_(makeRange(patchGhosts), density, flux, layout);
392+
interpolate_(makeRange(patchGhosts), particleDensity, chargeDensity, flux, layout);
392393

393394
if (level.getLevelNumber() > 0) // no levelGhost on root level
394395
{
395396
// then grab levelGhostParticlesOld and levelGhostParticlesNew
396397
// and project them with alpha and (1-alpha) coefs, respectively
397398
auto& levelGhostOld = pop.levelGhostParticlesOld();
398-
interpolate_(makeRange(levelGhostOld), density, flux, layout, 1. - alpha);
399+
interpolate_(makeRange(levelGhostOld), particleDensity, chargeDensity, flux, layout, 1. - alpha);
399400

400401
auto& levelGhostNew = pop.levelGhostParticlesNew();
401-
interpolate_(makeRange(levelGhostNew), density, flux, layout, alpha);
402+
interpolate_(makeRange(levelGhostNew), particleDensity, chargeDensity, flux, layout, alpha);
402403
}
403404
}
404405
}
@@ -538,7 +539,7 @@ namespace amr
538539

539540
auto& J = hybridModel.state.J;
540541
auto& Vi = hybridModel.state.ions.velocity();
541-
auto& Ni = hybridModel.state.ions.density();
542+
auto& Ni = hybridModel.state.ions.chargeDensity();
542543

543544
Jold_.copyData(J);
544545
ViOld_.copyData(Vi);

‎src/amr/physical_models/hybrid_model.hpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -154,7 +154,7 @@ void HybridModel<GridLayoutT, Electromag, Ions, Electrons, AMR_Types, Grid_t>::f
154154

155155
hybridInfo.modelMagnetic = core::VecFieldNames{state.electromag.B};
156156
hybridInfo.modelElectric = core::VecFieldNames{state.electromag.E};
157-
hybridInfo.modelIonDensity = state.ions.densityName();
157+
hybridInfo.modelIonDensity = state.ions.chargeDensityName();
158158
hybridInfo.modelIonBulkVelocity = core::VecFieldNames{state.ions.velocity()};
159159
hybridInfo.modelCurrent = core::VecFieldNames{state.J};
160160

‎src/amr/tagging/default_hybrid_tagger_strategy.hpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ void DefaultHybridTaggerStrategy<HybridModel>::tag(HybridModel& model,
3636
auto& By = model.state.electromag.B.getComponent(PHARE::core::Component::Y);
3737
auto& Bz = model.state.electromag.B.getComponent(PHARE::core::Component::Z);
3838

39-
auto& N = model.state.ions.density();
39+
auto& N = model.state.ions.chargeDensity();
4040

4141
// we loop on cell indexes for all qties regardless of their centering
4242
auto const& [start_x, _]

0 commit comments

Comments
 (0)