Skip to content

Expression unit_integral to normalize sources - #3481

Open
bendudson wants to merge 3 commits into
nextfrom
normalise-expr
Open

Expression unit_integral to normalize sources#3481
bendudson wants to merge 3 commits into
nextfrom
normalise-expr

Conversation

@bendudson

Copy link
Copy Markdown
Contributor

In expressions unit_integral(expr) will normalise the given expr so that the volume integral over the domain is 1. This is intended to make specifying sources easier, removing the need to manually scale sources to the required power or fuelling rate. This is a common confusion and errors when running source-driven cases.

Some care will still be needed when using this: The integral will use the normalized metrics so the input expression needs to normalise the volume; If a fraction of the torus is being simulated then that needs to be corrected for.

An alternative would be to implement in Hermes-3 more custom inputs e.g.

full_torus_power = 1e6 # Heating power for a full torus i.e. a 1/4 torus simulation would receive 1/4 of this power
power_shape = is_y_periodic * exp(-x^2) # The shape of the power deposition profile.

The unit_integral() function is intended as a more general solution.

Stores the CELL_LOC and allows generators to retrieve
it using a `location() const` method.
In expressions `unit_integral(expr)` normalises the given `expr` so
that the volume integral over the domain is 1. This is intended
to make specifying sources easier, removing the need to manually
scale sources to the required power or fuelling rate.
Brief description in the table of available functions.

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

clang-tidy made some suggestions


void FieldUnitIntegral::populateCache(const Context& ctx) {
Mesh* localmesh = ctx.getMesh();
ASSERT0(localmesh != nullptr);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "ASSERT0" is directly included [misc-include-cleaner]

src/field/fieldgenerators.cxx:3:

- #include <bout/boutcomm.hxx>
+ #include "bout/assert.hxx"
+ #include <bout/boutcomm.hxx>

void FieldUnitIntegral::populateCache(const Context& ctx) {
Mesh* localmesh = ctx.getMesh();
ASSERT0(localmesh != nullptr);
Coordinates* coords = localmesh->getCoordinates(ctx.location());

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "Coordinates" is directly included [misc-include-cleaner]

src/field/fieldgenerators.cxx:3:

- #include <bout/boutcomm.hxx>
+ #include "bout/coordinates.hxx"
+ #include <bout/boutcomm.hxx>

void FieldUnitIntegral::populateCache(const Context& ctx) {
Mesh* localmesh = ctx.getMesh();
ASSERT0(localmesh != nullptr);
Coordinates* coords = localmesh->getCoordinates(ctx.location());

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: pointee of variable 'coords' of type 'Coordinates *' can be declared 'const' [misc-const-correctness]

Suggested change
Coordinates* coords = localmesh->getCoordinates(ctx.location());
Coordinates const* coords = localmesh->getCoordinates(ctx.location());

Coordinates* coords = localmesh->getCoordinates(ctx.location());
if (coords == nullptr) {
throw BoutException("unit_integral function needs coordinates at {}",
toString(ctx.location()));

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "toString" is directly included [misc-include-cleaner]

                        toString(ctx.location()));
                        ^


cached_values = Field3D(localmesh).setLocation(ctx.location()).allocate();

BOUT_FOR(i, cached_values.getRegion("RGN_ALL")) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "BOUT_FOR" is directly included [misc-include-cleaner]

src/field/fieldgenerators.cxx:3:

- #include <bout/boutcomm.hxx>
+ #include "bout/region.hxx"
+ #include <bout/boutcomm.hxx>


BoutReal FieldUnitIntegral::generate(const Context& ctx) {
if (!cacheMatches(ctx)) {
std::lock_guard<std::mutex> guard(cache_mutex);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: variable 'guard' of type 'std::lock_guardstd::mutex' can be declared 'const' [misc-const-correctness]

Suggested change
std::lock_guard<std::mutex> guard(cache_mutex);
std::lock_guard<std::mutex> const guard(cache_mutex);

}
}

ASSERT1(ctx.location() == cached_location);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "ASSERT1" is directly included [misc-include-cleaner]

  ASSERT1(ctx.location() == cached_location);
  ^

public:
explicit FieldUnitIntegral(FieldGeneratorPtr g = nullptr) : gen(std::move(g)) {}

FieldGeneratorPtr clone(const std::list<FieldGeneratorPtr> args) override;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: parameter 'args' is const-qualified in the function declaration; const-qualification of parameters only has an effect in function definitions [readability-avoid-const-params-in-decls]

Suggested change
FieldGeneratorPtr clone(const std::list<FieldGeneratorPtr> args) override;
FieldGeneratorPtr clone(std::list<FieldGeneratorPtr> args) override;

explicit FieldUnitIntegral(FieldGeneratorPtr g = nullptr) : gen(std::move(g)) {}

FieldGeneratorPtr clone(const std::list<FieldGeneratorPtr> args) override;
BoutReal generate(const bout::generator::Context& pos) override;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: function 'FieldUnitIntegral::generate' has a definition with different parameter names [readability-inconsistent-declaration-parameter-name]

  BoutReal generate(const bout::generator::Context& pos) override;
           ^
Additional context

src/field/fieldgenerators.cxx:101: the definition seen here

BoutReal FieldUnitIntegral::generate(const Context& ctx) {
                            ^

src/field/fieldgenerators.hxx:273: differing parameters are named here: ('pos'), in definition: ('ctx')

  BoutReal generate(const bout::generator::Context& pos) override;
           ^

bool cacheMatches(const bout::generator::Context& ctx) const;

FieldGeneratorPtr gen;
Field3D cached_values{};

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: initializer for member 'cached_values' is redundant [readability-redundant-member-init]

Suggested change
Field3D cached_values{};
Field3D cached_values;

@bendudson bendudson added enhancement feature A PR that adds new capabilities / features and removed enhancement labels Aug 27, 2026

@ZedThree ZedThree left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good, a few comments

BoutReal local = 0.0;
BOUT_FOR(i, field.getRegion("RGN_NOBNDRY")) {
local +=
field[i] * coords->J()[i] * coords->dx()[i] * coords->dy()[i] * coords->dz()[i];

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
field[i] * coords->J()[i] * coords->dx()[i] * coords->dy()[i] * coords->dz()[i];
field[i] * coords->cell_volume[i];

Field3D ones{localmesh};
ones.allocate();
BOUT_FOR(i, ones.getRegion("RGN_NOBNDRY")) {
local += coords->J()[i] * coords->dx()[i] * coords->dy()[i] * coords->dz()[i];

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
local += coords->J()[i] * coords->dx()[i] * coords->dy()[i] * coords->dz()[i];
local += coords->cell_volume[i];

using bout::generator::Context;

namespace {
BoutReal volumeIntegral(const Field3D& field) {

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm surprised that we don't already have a sum(Field3D) given that we have mean -- if we add that, then these two functions are just sum(field * coords->cell_volume()) and sum(coords->cell_volume()).

We could then use this in the generator implementation?


BoutReal FieldUnitIntegral::generate(const Context& ctx) {
if (!cacheMatches(ctx)) {
std::lock_guard<std::mutex> guard(cache_mutex);

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this AI by any chance? We don't have usual threads, we use OpenMP, so this should be BOUT_OMP(single) instead?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

feature A PR that adds new capabilities / features

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants