Skip to content

Commit

Permalink
[oneD] Replace Sim1D::m_x by OneDim::m_data
Browse files Browse the repository at this point in the history
  • Loading branch information
@ischoegl @speth
ischoegl authored and speth committed Apr 9, 2023
1 parent 96ba557 commit b735e75
Show file tree
Hide file tree
Showing 4 changed files with 51 additions and 46 deletions.
2 changes: 2 additions & 0 deletions include/cantera/oneD/OneDim.h
View file
Expand Up @@ -339,6 +339,8 @@ class OneDim
//! factor time step is multiplied by if time stepping fails ( < 1 )
double m_tfactor = 0.5;

shared_ptr<vector<double>> m_data; //!< Solution vector

std::unique_ptr<MultiJac> m_jac; //!< Jacobian evaluator
std::unique_ptr<MultiNewton> m_newt; //!< Newton iterator
double m_rdt = 0.0; //!< reciprocal of time step
Expand Down
13 changes: 5 additions & 8 deletions include/cantera/oneD/Sim1D.h
View file
Expand Up @@ -196,20 +196,20 @@ class Sim1D : public OneDim
const doublereal* solution() {
warn_deprecated("Sim1D::solution",
"This method is unused and will be removed after Cantera 3.0.");
return m_x.data();
return m_data->data();
}

void setTimeStep(double stepsize, size_t n, const int* tsteps);

void solve(int loglevel = 0, bool refine_grid = true);

void eval(doublereal rdt=-1.0, int count = 1) {
OneDim::eval(npos, m_x.data(), m_xnew.data(), rdt, count);
OneDim::eval(npos, m_data->data(), m_xnew.data(), rdt, count);
}

// Evaluate the governing equations and return the vector of residuals
void getResidual(double rdt, double* resid) {
OneDim::eval(npos, m_x.data(), resid, rdt, 0);
OneDim::eval(npos, m_data->data(), resid, rdt, 0);
}

//! Refine the grid in all domains.
Expand Down Expand Up @@ -278,14 +278,14 @@ class Sim1D : public OneDim
void setSolution(const doublereal* soln) {
warn_deprecated("Sim1D::setSolution",
"This method is unused and will be removed after Cantera 3.0.");
std::copy(soln, soln + m_x.size(), m_x.data());
std::copy(soln, soln + m_data->size(), m_data->data());
}

// @deprecated To be removed after Cantera 3.0 (unused)
const doublereal* solution() const {
warn_deprecated("Sim1D::solution",
"This method is unused and will be removed after Cantera 3.0.");
return m_x.data();
return m_data->data();
}

doublereal jacobian(int i, int j);
Expand Down Expand Up @@ -317,9 +317,6 @@ class Sim1D : public OneDim
}

protected:
//! the solution vector
vector_fp m_x;

//! the solution vector after the last successful timestepping
vector_fp m_xlast_ts;

Expand Down
6 changes: 6 additions & 0 deletions src/oneD/OneDim.cpp
View file
Expand Up @@ -237,6 +237,12 @@ void OneDim::resize()
m_size = d->loc() + d->size();
}

if (m_data) {
m_data->resize(size());
} else {
m_data = make_shared<vector<double>>(size());
}

m_newt->resize(size());
m_mask.resize(size());

Expand Down
76 changes: 38 additions & 38 deletions src/oneD/Sim1D.cpp
View file
Expand Up @@ -30,7 +30,7 @@ Sim1D::Sim1D(vector<shared_ptr<Domain1D>>& domains) :
// domain-specific initialization of the solution vector.
resize();
for (size_t n = 0; n < nDomains(); n++) {
domain(n)._getInitialSoln(&m_x[start(n)]);
domain(n)._getInitialSoln(m_data->data() + start(n));
}

// set some defaults
Expand All @@ -50,7 +50,7 @@ Sim1D::Sim1D(vector<Domain1D*>& domains) :
// domain-specific initialization of the solution vector.
resize();
for (size_t n = 0; n < nDomains(); n++) {
domain(n)._getInitialSoln(&m_x[start(n)]);
domain(n)._getInitialSoln(m_data->data() + start(n));
}

// set some defaults
Expand All @@ -74,24 +74,24 @@ void Sim1D::setInitialGuess(const std::string& component, vector_fp& locs, vecto
void Sim1D::setValue(size_t dom, size_t comp, size_t localPoint, doublereal value)
{
size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
AssertThrowMsg(iloc < m_x.size(), "Sim1D::setValue",
"Index out of bounds: {} > {}", iloc, m_x.size());
m_x[iloc] = value;
AssertThrowMsg(iloc < m_data->size(), "Sim1D::setValue",
"Index out of bounds: {} > {}", iloc, m_data->size());
(*m_data)[iloc] = value;
}

doublereal Sim1D::value(size_t dom, size_t comp, size_t localPoint) const
{
size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
AssertThrowMsg(iloc < m_x.size(), "Sim1D::value",
"Index out of bounds: {} > {}", iloc, m_x.size());
return m_x[iloc];
AssertThrowMsg(iloc < m_data->size(), "Sim1D::value",
"Index out of bounds: {} > {}", iloc, m_data->size());
return (*m_data)[iloc];
}

doublereal Sim1D::workValue(size_t dom, size_t comp, size_t localPoint) const
{
size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
AssertThrowMsg(iloc < m_x.size(), "Sim1D::workValue",
"Index out of bounds: {} > {}", iloc, m_x.size());
AssertThrowMsg(iloc < m_data->size(), "Sim1D::workValue",
"Index out of bounds: {} > {}", iloc, m_data->size());
return m_xnew[iloc];
}

Expand Down Expand Up @@ -133,7 +133,7 @@ void Sim1D::save(const std::string& fname, const std::string& id,
if (extension == "h5" || extension == "hdf" || extension == "hdf5") {
SolutionArray::writeHeader(fname, id, desc, overwrite);
for (auto dom : m_dom) {
auto arr = dom->asArray(m_x.data() + dom->loc());
auto arr = dom->asArray(m_data->data() + dom->loc());
arr->writeEntry(fname, id, dom->id(), overwrite, compression);
}
return;
Expand All @@ -147,7 +147,7 @@ void Sim1D::save(const std::string& fname, const std::string& id,
SolutionArray::writeHeader(data, id, desc, overwrite);

for (auto dom : m_dom) {
auto arr = dom->asArray(m_x.data() + dom->loc());
auto arr = dom->asArray(m_data->data() + dom->loc());
arr->writeEntry(data, id, dom->id(), overwrite);
}

Expand All @@ -174,13 +174,14 @@ void Sim1D::saveResidual(const std::string& fname, const std::string& id,
void Sim1D::saveResidual(const std::string& fname, const std::string& id,
const std::string& desc, bool overwrite, int compression)
{
vector_fp res(m_x.size(), -999);
OneDim::eval(npos, &m_x[0], &res[0], 0.0);
// Temporarily put the residual into m_x, since this is the vector that the save()
// function reads.
std::swap(res, m_x);
vector_fp res(m_data->size(), -999);
OneDim::eval(npos, m_data->data(), &res[0], 0.0);
// Temporarily put the residual into m_data, since this is the vector that the
// save() function reads.
vector<double> backup(*m_data);
*m_data = res;
save(fname, id, desc, overwrite, compression);
std::swap(res, m_x);
*m_data = backup;
}

namespace { // restrict scope of helper function to local translation unit
Expand Down Expand Up @@ -302,7 +303,7 @@ AnyMap Sim1D::restore(const std::string& fname, const std::string& id)
resize();
m_xlast_ts.clear();
for (auto dom : m_dom) {
dom->restore(*arrs[dom->id()], m_x.data() + dom->loc());
dom->restore(*arrs[dom->id()], m_data->data() + dom->loc());
}
finalize();
} else if (extension == "yaml" || extension == "yml") {
Expand All @@ -319,7 +320,7 @@ AnyMap Sim1D::restore(const std::string& fname, const std::string& id)
resize();
m_xlast_ts.clear();
for (auto dom : m_dom) {
dom->restore(*arrs[dom->id()], m_x.data() + dom->loc());
dom->restore(*arrs[dom->id()], m_data->data() + dom->loc());
}
finalize();
} else {
Expand Down Expand Up @@ -356,7 +357,7 @@ void Sim1D::show(ostream& s)
{
for (size_t n = 0; n < nDomains(); n++) {
if (domain(n).type() != "empty") {
domain(n).show(s, &m_x[start(n)]);
domain(n).show(s, m_data->data() + start(n));
}
}
}
Expand All @@ -367,7 +368,7 @@ void Sim1D::show()
if (domain(n).type() != "empty") {
writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id()
+" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
domain(n).show(&m_x[start(n)]);
domain(n).show(m_data->data() + start(n));
}
}
}
Expand All @@ -378,7 +379,7 @@ void Sim1D::restoreTimeSteppingSolution()
throw CanteraError("Sim1D::restoreTimeSteppingSolution",
"No successful time steps taken on this grid.");
}
m_x = m_xlast_ts;
*m_data = m_xlast_ts;
}

void Sim1D::restoreSteadySolution()
Expand All @@ -387,7 +388,7 @@ void Sim1D::restoreSteadySolution()
throw CanteraError("Sim1D::restoreSteadySolution",
"No successful steady state solution");
}
m_x = m_xlast_ss;
*m_data = m_xlast_ss;
for (size_t n = 0; n < nDomains(); n++) {
vector_fp& z = m_grid_last_ss[n];
domain(n).setupGrid(z.size(), z.data());
Expand All @@ -397,14 +398,14 @@ void Sim1D::restoreSteadySolution()
void Sim1D::getInitialSoln()
{
for (size_t n = 0; n < nDomains(); n++) {
domain(n)._getInitialSoln(&m_x[start(n)]);
domain(n)._getInitialSoln(m_data->data() + start(n));
}
}

void Sim1D::finalize()
{
for (size_t n = 0; n < nDomains(); n++) {
domain(n)._finalize(&m_x[start(n)]);
domain(n)._finalize(m_data->data() + start(n));
}
}

Expand All @@ -419,9 +420,9 @@ void Sim1D::setTimeStep(double stepsize, size_t n, const int* tsteps)

int Sim1D::newtonSolve(int loglevel)
{
int m = OneDim::solve(m_x.data(), m_xnew.data(), loglevel);
int m = OneDim::solve(m_data->data(), m_xnew.data(), loglevel);
if (m >= 0) {
m_x = m_xnew;
*m_data = m_xnew;
return 0;
} else if (m > -10) {
return -1;
Expand Down Expand Up @@ -493,9 +494,8 @@ void Sim1D::solve(int loglevel, bool refine_grid)
if (loglevel > 0) {
writelog("Take {} timesteps ", nsteps);
}
dt = timeStep(nsteps, dt, m_x.data(), m_xnew.data(),
loglevel-1);
m_xlast_ts = m_x;
dt = timeStep(nsteps, dt, m_data->data(), m_xnew.data(), loglevel-1);
m_xlast_ts = *m_data;
if (loglevel > 6) {
save("debug_sim1d.yaml", "debug", "After timestepping");
}
Expand All @@ -506,7 +506,7 @@ void Sim1D::solve(int loglevel, bool refine_grid)

if (loglevel == 1) {
writelog(" {:10.4g} {:10.4g}\n", dt,
log10(ssnorm(m_x.data(), m_xnew.data())));
log10(ssnorm(m_data->data(), m_xnew.data())));
}
istep++;
if (istep >= m_steps.size()) {
Expand Down Expand Up @@ -551,7 +551,7 @@ int Sim1D::refine(int loglevel)
vector_fp znew, xnew;
std::vector<size_t> dsize;

m_xlast_ss = m_x;
m_xlast_ss = *m_data;
m_grid_last_ss.clear();

for (size_t n = 0; n < nDomains(); n++) {
Expand All @@ -562,7 +562,8 @@ int Sim1D::refine(int loglevel)
m_grid_last_ss.push_back(d.grid());

// determine where new points are needed
ianalyze = r.analyze(d.grid().size(), d.grid().data(), &m_x[start(n)]);
ianalyze = r.analyze(d.grid().size(), d.grid().data(),
m_data->data() + start(n));
if (ianalyze < 0) {
return ianalyze;
}
Expand Down Expand Up @@ -627,7 +628,7 @@ int Sim1D::refine(int loglevel)
}

// Replace the current solution vector with the new one
m_x = xnew;
*m_data = xnew;
resize();
finalize();
return np;
Expand Down Expand Up @@ -716,7 +717,7 @@ int Sim1D::setFixedTemperature(double t)
}

// Replace the current solution vector with the new one
m_x = xnew;
*m_data = xnew;

resize();
finalize();
Expand Down Expand Up @@ -813,7 +814,7 @@ doublereal Sim1D::jacobian(int i, int j)

void Sim1D::evalSSJacobian()
{
OneDim::evalSSJacobian(m_x.data(), m_xnew.data());
OneDim::evalSSJacobian(m_data->data(), m_xnew.data());
}

void Sim1D::solveAdjoint(const double* b, double* lambda)
Expand Down Expand Up @@ -843,7 +844,6 @@ void Sim1D::solveAdjoint(const double* b, double* lambda)
void Sim1D::resize()
{
OneDim::resize();
m_x.resize(size(), 0.0);
m_xnew.resize(size(), 0.0);
}

Expand Down

0 comments on commit b735e75

Please sign in to comment.