Skip to content

Commit

Permalink
new bc-s (WIP)
Browse files Browse the repository at this point in the history
  • Loading branch information
haykh committed Nov 21, 2024
1 parent cb7b179 commit 77d66a8
Show file tree
Hide file tree
Showing 11 changed files with 886 additions and 316 deletions.
18 changes: 13 additions & 5 deletions input.example.toml
Original file line number Diff line number Diff line change
Expand Up @@ -90,11 +90,11 @@
# Boundary conditions for fields:
# @required
# @type: 1/2/3-size array of string tuples, each of size 1 or 2
# @valid: "PERIODIC", "ABSORB", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON"
# @example: [["CUSTOM", "ABSORB"]] (for 2D spherical [[rmin, rmax]])
# @valid: "PERIODIC", "MATCH", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON"
# @example: [["CUSTOM", "MATCH"]] (for 2D spherical [[rmin, rmax]])
# @note: When periodic in any of the directions, you should only set one value: [..., ["PERIODIC"], ...]
# @note: In spherical, bondaries in theta/phi are set automatically (only specify bc @ [rmin, rmax]): [["ATMOSPHERE", "ABSORB"]]
# @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["ABSORB"]]
# @note: In spherical, bondaries in theta/phi are set automatically (only specify bc @ [rmin, rmax]): [["ATMOSPHERE", "match"]]
# @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["match"]]
fields = ""
# Boundary conditions for fields:
# @required
Expand All @@ -106,7 +106,7 @@
# @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["ABSORB"]]
particles = ""

[grid.boundaries.absorb]
[grid.boundaries.match]
# Size of the absorption layer in physical (code) units:
# @type: float
# @default: 1% of the domain size (in shortest dimension)
Expand All @@ -118,6 +118,14 @@
# @default: 1.0
coeff = ""

[grid.boundaries.absorb]
# Size of the absorption layer for particles in physical (code) units:
# @type: float
# @default: 1% of the domain size (in shortest dimension)
# @note: In spherical, this is the size of the layer in r from the outer wall
# @note: In cartesian, this is the same for all dimensions where applicable
ds = ""

[grid.boundaries.atmosphere]
# @required: if ATMOSPHERE is one of the boundaries
# Temperature of the atmosphere in units of m0 c^2
Expand Down
19 changes: 11 additions & 8 deletions setups/srpic/shock/pgen.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@
#include "archetypes/problem_generator.h"
#include "framework/domain/metadomain.h"

#include <utility>

namespace user {
using namespace ntt;

Expand Down Expand Up @@ -93,20 +95,21 @@ namespace user {

inline PGen() {}

auto FixField(const em& comp) const -> real_t {
auto FixFieldsConst(const bc_in&, const em& comp) const
-> std::pair<real_t, bool> {
if (comp == em::ex2) {
return init_flds.ex2({ ZERO });
return { init_flds.ex2({ ZERO }), true };
} else if (comp == em::ex3) {
return init_flds.ex3({ ZERO });
} else if (comp == em::bx1) {
return init_flds.bx1({ ZERO });
return { init_flds.ex3({ ZERO }), true };
} else {
raise::Error("Other components should not be requested when BC is in X",
HERE);
return ZERO;
return { ZERO, false };
}
}

auto MatchFields(real_t time) const -> InitFields<D> {
return init_flds;
}

inline void InitPrtls(Domain<S, M>& local_domain) {
const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
local_domain.random_pool,
Expand Down
170 changes: 170 additions & 0 deletions setups/wip/magpump/pgen.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,170 @@
#ifndef PROBLEM_GENERATOR_H
#define PROBLEM_GENERATOR_H

#include "enums.h"
#include "global.h"

#include "arch/traits.h"

#include "archetypes/particle_injector.h"
#include "archetypes/problem_generator.h"
#include "framework/domain/metadomain.h"

#include <plog/Log.h>

namespace user {
using namespace ntt;

template <Dimension D>
struct InitFields {
InitFields(real_t bsurf, real_t rstar) : Bsurf { bsurf }, Rstar { rstar } {}

Inline auto bx1(const coord_t<D>& x_Ph) const -> real_t {
return Bsurf * math::cos(x_Ph[1]) / CUBE(x_Ph[0] / Rstar);
}

Inline auto bx2(const coord_t<D>& x_Ph) const -> real_t {
return Bsurf * HALF * math::sin(x_Ph[1]) / CUBE(x_Ph[0] / Rstar);
}

private:
const real_t Bsurf, Rstar;
};

template <Dimension D>
struct DriveFields : public InitFields<D> {
DriveFields(real_t time, real_t bsurf, real_t rstar)
: InitFields<D> { bsurf, rstar }
, time { time } {}

using InitFields<D>::bx1;
using InitFields<D>::bx2;

Inline auto bx3(const coord_t<D>&) const -> real_t {
return ZERO;
}

Inline auto ex1(const coord_t<D>&) const -> real_t {
return ZERO;
}

Inline auto ex2(const coord_t<D>& x_Ph) const -> real_t {
return ZERO;
}

Inline auto ex3(const coord_t<D>&) const -> real_t {
return ZERO;
}

private:
const real_t time;
};

template <SimEngine::type S, class M>
struct Inflow : public arch::EnergyDistribution<S, M> {
Inflow(const M& metric, real_t vin)
: arch::EnergyDistribution<S, M> { metric }
, vin { vin } {}

Inline void operator()(const coord_t<M::Dim>&,
vec_t<Dim::_3D>& v_Ph,
unsigned short) const override {
v_Ph[0] = -vin;
}

private:
const real_t vin;
};

template <SimEngine::type S, class M>
struct Sphere : public arch::SpatialDistribution<S, M> {
Sphere(const M& metric, real_t r0, real_t dr)
: arch::SpatialDistribution<S, M> { metric }
, r0 { r0 }
, dr { dr } {}

Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t override {
return math::exp(-SQR((x_Ph[0] - r0) / dr)) *
(x_Ph[1] > 0.25 && x_Ph[1] < constant::PI - 0.25);
}

private:
const real_t r0, dr;
};

template <SimEngine::type S, class M>
struct PGen : public arch::ProblemGenerator<S, M> {
static constexpr auto engines { traits::compatible_with<SimEngine::SRPIC>::value };
static constexpr auto metrics {
traits::compatible_with<Metric::Spherical, Metric::QSpherical>::value
};
static constexpr auto dimensions { traits::compatible_with<Dim::_2D>::value };

using arch::ProblemGenerator<S, M>::D;
using arch::ProblemGenerator<S, M>::C;
using arch::ProblemGenerator<S, M>::params;

const real_t Bsurf, pump_period, pump_ampl, pump_radius, Rstar;
const real_t vin, drinj;
InitFields<D> init_flds;

inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
: arch::ProblemGenerator<S, M> { p }
, Bsurf { p.template get<real_t>("setup.Bsurf", ONE) }
, pump_period { p.template get<real_t>("setup.pump_period") }
, pump_ampl { p.template get<real_t>("setup.pump_ampl") }
, pump_radius { p.template get<real_t>("setup.pump_radius") }
, Rstar { m.mesh().extent(in::x1).first }
, vin { p.template get<real_t>("setup.vin") }
, drinj { p.template get<real_t>("setup.drinj") }
, init_flds { Bsurf, Rstar } {}

auto FieldDriver(real_t time) const -> DriveFields<D> {
return DriveFields<D> { time, Bsurf, Rstar };
}

void CustomPostStep(std::size_t, long double time, Domain<S, M>& domain) {
const real_t radius = pump_radius +
pump_ampl *
math::sin(time * constant::TWO_PI / pump_period);
const real_t dr = 1.0;
const auto& metric = domain.mesh.metric;
auto EM = domain.fields.em;
Kokkos::parallel_for(
"outerBC",
domain.mesh.rangeActiveCells(),
Lambda(index_t i1, index_t i2) {
const auto i1_ = COORD(i1), i2_ = COORD(i2);
const auto r = metric.template convert<1, Crd::Cd, Crd::Ph>(i1_);
if (r > radius - 5 * dr) {
const auto smooth = HALF * (ONE - math::tanh((r - radius) / dr));
EM(i1, i2, em::ex1) = smooth * EM(i1, i2, em::ex1);
EM(i1, i2, em::ex2) = smooth * EM(i1, i2, em::ex2);
EM(i1, i2, em::ex3) = smooth * EM(i1, i2, em::ex3);
EM(i1, i2, em::bx1) = smooth * EM(i1, i2, em::bx1);
EM(i1, i2, em::bx2) = smooth * EM(i1, i2, em::bx2);
EM(i1, i2, em::bx3) = smooth * EM(i1, i2, em::bx3);
}
});

if (time < pump_period * 0.25) {
const auto energy_dist = Inflow<S, M>(domain.mesh.metric, vin);
const auto spatial_dist = Sphere<S, M>(domain.mesh.metric, radius, drinj);
const auto injector = arch::NonUniformInjector<S, M, Inflow, Sphere>(
energy_dist,
spatial_dist,
{ 1, 2 });

arch::InjectNonUniform<S, M, arch::NonUniformInjector<S, M, Inflow, Sphere>>(
params,
domain,
injector,
ONE,
true);
}
}
};

} // namespace user

#endif
Loading

0 comments on commit 77d66a8

Please sign in to comment.