Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change PI density controller integrator #10

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Change PI density controller integrator #10

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 19 additions & 31 deletions sd1d.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -215,11 +215,10 @@ class SD1D : public PhysicsModel {
OPTION(opt, density_integral_positive, false);
OPTION(opt, density_source_positive, true);

density_error_lasttime = -1.0; // Signal no value

// Save and load error integral from file, since
// this determines the source function
restart.add(density_error_integral, "density_error_integral");
// Evolve the error integral
// Note: Only the first point in the domain
// (upstream density) is actually evolved
SOLVE_FOR(density_error_integral);

if (!restarting) {
density_error_integral = 0.0;
Expand Down Expand Up @@ -786,6 +785,8 @@ class SD1D : public PhysicsModel {

TRACE("Density upstream");

ddt(density_error_integral) = 0.0; // Most of the domain not evolving

BoutReal source;
for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) {
int jz = 0;
Expand All @@ -796,33 +797,22 @@ class SD1D : public PhysicsModel {
ASSERT2(finite(error));
ASSERT2(finite(density_error_integral));

// PI controller, using crude integral of the error
if (density_error_lasttime < 0.0) {
// First time
density_error_lasttime = time;
density_error_last = error;
}

// Integrate using Trapezium rule
if (time > density_error_lasttime) { // Since time can decrease
density_error_integral += (time - density_error_lasttime) * 0.5 *
(error + density_error_last);
}

if ((density_error_integral < 0.0) && density_integral_positive) {
// Limit density_error_integral to be >= 0
density_error_integral = 0.0;
// The error integral in the first cell
BoutReal density_error_integral_upstream = density_error_integral(r.ind, mesh->ystart);

if (!density_integral_positive || (density_error_integral_upstream > 0.0) || (error > 0)) {
// If density_integral_positive = true
// and density_error_integral < 0
// and error < 0
// then don't decrease the error integral further
//
// Otherwise (the default), evolve the first cell
ddt(density_error_integral)(r.ind, mesh->ystart) = error;
}

// Calculate source from combination of error and integral
source = density_controller_p * error +
density_controller_i * density_error_integral;

// output.write("\n Source: %e, %e : %e, %e -> %e\n", time, (time -
// density_error_lasttime), error, density_error_integral, source);

density_error_last = error;
density_error_lasttime = time;
density_controller_i * density_error_integral_upstream;

if (!volume_source) {
// Convert source into a flow velocity
Expand Down Expand Up @@ -1884,9 +1874,7 @@ class SD1D : public PhysicsModel {
bool density_integral_positive; // Limit the i term to be positive
bool density_source_positive; // Limit the source to be positive

BoutReal density_error_lasttime,
density_error_last; // Value and time of last error
BoutReal density_error_integral; // Integral of error
Field2D density_error_integral; // Integral of error

///////////////////////////////////////////////////////////////
// Numerical dissipation
Expand Down