I was reviewing some of the less invasive changes that I had intertwined with the two-point continuation PR, and decided that I would try to get them together into a separate PR. I believe these changes do help with the overall readability of the StFlow methods regarding residual evaluation. One particular fix that I wanted to do was to move the continuity and right side residual functions away from the very bottom of the file closer to where they are called. After spending some time in the StFlow.cpp file, scrolling up and down to examine those functions was becoming tiresome.

@wandadars ... thanks for sorting through some of this code. Now that Cantera 3.0 is released, I don't anticipate any new merge conflicts.

While I am a little hesitant overall (the pre-existing code is somewhat brittle and can be easily broken), I think that there is potential for significant simplifications here. One change I am requesting is to move the for loop -

As an aside, I'd recommend starting a new commit for each new helper function introduced (and making sure scons test passes locally every time; if you want to save time, just run scons test-python-onedim as the Python test suite is the most comprehensive) to ensure that there are no inadvertent breakages.

PS: there are VScode settings that remove trailing whitespace automatically whenever you save.

Thanks for the suggestions @ischoegl. I just want to be sure I understand your suggestion before jumping back into this pull request again. Are you suggesting that we have a helper methods like:
evalContinuity()
evalEnergy()
evalLambda()
evalMomentum()

And that inside of each of these would be the for loop that would loop over every point in the 1D domain, and then update the appropriate entry in the rsd variable. Is that an accurate summary of your suggestion?

Is that an accurate summary of your suggestion?

👍 it indeed is!

The signature of evalContinuity would change from the current version, i.e. rather than

    virtual void evalContinuity(size_t j, double* x, double* r, int* diag, double rdt);

you'd need

    virtual void evalContinuity(double* x, double* rsd, int* diag, 
                                double rdt, size_t jmin, size_t jmax);

or equivalent. The other new methods, - evalEnergy(...), evalLambda(...), evalMomentum(...), evalSpecies(...), etc., - use the same interface. (on an unrelated note, the diag variable name is not an intuitive choice; Domain1d::eval uses mask, which is more appropriate)

I had stared at the code for a while and could not see any interactions that prevent this approach.

PS: I realize that this is more than a 'minor' update, although it's not too complicated - don't hesitate to reach out

I would just say, be careful about looping through the grid points multiple times. If you do that, I think you need to make sure all properties from the ThermoPhase and Kinetics object are calculated in a single pass and stored in the StFlow object. Otherwise, you risk doing a lot of extra work by resetting the phase state and destroying the ability to re-use any partial/cached data that those objects have. Handling this correctly may reduce the benefit of the suggested re-organization.

I would just say, be careful about looping through the grid points multiple times. If you do that, I think you need to make sure all properties from the ThermoPhase and Kinetics object are calculated in a single pass and stored in the StFlow object. Otherwise, you risk doing a lot of extra work by resetting the phase state and destroying the ability to re-use any partial/cached data that those objects have. Handling this correctly may reduce the benefit of the suggested re-organization.

@speth ... definitely a fair point, but the overall approach appears to be that property data are precalculated? (there is the separate call to updateProperties, which itself redirects to updateThermo, updateTransport and updateDiffFluxes). If this is executed properly, there shouldn't be an issue. But ...

It is absolutely true that the calls to getWdot in the species equation and setGas in the energy equation (which I'm not sure is even necessary) need to be addressed (I believe that these are the only ones); at the same time, there already is a m_wdot member variable, where the current implementation populates entries in a fashion that is not consistent with updateThermo (it almost looks like an oversight?).

@wandadars ... let me know if you want to give this a go. Otherwise, I may go ahead and take a stab myself.

@ischoegl I'm currently giving these extra updates that you suggested an attempt. I have them in, but I'm troubleshooting a failing Python one dim test.

@ischoegl I'm currently giving these extra updates that you suggested an attempt. I have them in, but I'm troubleshooting a failing Python one dim test.

Should have looked at the commits in a little more detail. Sorry!

I would just say, be careful about looping through the grid points multiple times. If you do that, I think you need to make sure all properties from the ThermoPhase and Kinetics object are calculated in a single pass and stored in the StFlow object. Otherwise, you risk doing a lot of extra work by resetting the phase state and destroying the ability to re-use any partial/cached data that those objects have. Handling this correctly may reduce the benefit of the suggested re-organization.

@speth ... definitely a fair point, but the overall approach appears to be that property data are precalculated? (there is the separate call to updateProperties, which itself redirects to updateThermo, updateTransport and updateDiffFluxes). If this is executed properly, there shouldn't be an issue. But ...

It is absolutely true that the calls to getWdot in the species equation and setGas in the energy equation (which I'm not sure is even necessary) need to be addressed (I believe that these are the only ones); at the same time, there already is a m_wdot member variable, where the current implementation populates entries in a fashion that is not consistent with updateThermo (it almost looks like an oversight?).

I think the source of my issue was when I re-ordered the evaluation of the governing equations to u,V,T,Y,L,E, and the species equation made a call to getWdot() that must have populated the m_wdot variable that is then used by the energy equation, thus moving the energy equation before the species equation results in bad values being used possibly. Should something like the getWdot() be called in the updateProperties() method or one of the methods that this method calls?

Edit: For the time being I just swapped the order that I called the residual methods evalSpecies() and evalEnergy() to put the species equation first, and it resolved my unit test issue.

We have:

//! Write the net production rates at point `j` into array `m_wdot`
    void getWdot(double* x, size_t j) {
        setGas(x,j);
        m_kin->getNetProductionRates(&m_wdot(0,j));
    }
    
    /**
     * Update the thermodynamic properties from point j0 to point j1
     * (inclusive), based on solution x.
     */
    void updateThermo(const double* x, size_t j0, size_t j1) {
        for (size_t j = j0; j <= j1; j++) {
            setGas(x,j);
            m_rho[j] = m_thermo->density();
            m_wtm[j] = m_thermo->meanMolecularWeight();
            m_cp[j] = m_thermo->cp_mass();
            m_thermo->getPartialMolarEnthalpies(&m_hk(0, j));
        }
    }

I know the production rate of a species isn't a thermodynamic quantity, but would it make sense to have the logic of the getWdot() method put into the updateProperties() method that is used to call the updateThermo() method?

For reference:

void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
{
    // properties are computed for grid points from j0 to j1
    size_t j0 = std::max<size_t>(jmin, 1) - 1;
    size_t j1 = std::min(jmax+1,m_points-1);

    updateThermo(x, j0, j1);
    if (jg == npos || m_force_full_update) {
        // update transport properties only if a Jacobian is not being
        // evaluated, or if specifically requested
        updateTransport(x, j0, j1);
    }
    if (jg == npos) {
        double* Yleft = x + index(c_offset_Y, jmin);
        m_kExcessLeft = distance(Yleft, max_element(Yleft, Yleft + m_nsp));
        double* Yright = x + index(c_offset_Y, jmax);
        m_kExcessRight = distance(Yright, max_element(Yright, Yright + m_nsp));
    }

    // update the species diffusive mass fluxes whether or not a
    // Jacobian is being evaluated
    updateDiffFluxes(x, j0, j1);
}

I know that the goal is to avoid having both the species and energy equations calling getWdot() as it just duplicates the work. I'm also wary of having implicit relation between the species and energy residual equations such that one must always be evaluated before the other.

Codecov Report

Merging #1595 (71db3af) into main (bc8ebe8) will increase coverage by 0.01%.
The diff coverage is 93.82%.

@@            Coverage Diff             @@
##             main    #1595      +/-   ##
==========================================
+ Coverage   72.70%   72.71%   +0.01%     
==========================================
  Files         370      370              
  Lines       56257    56286      +29     
  Branches    20351    20368      +17     
==========================================
+ Hits        40903    40931      +28     
- Misses      12357    12358       +1     
  Partials     2997     2997              

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

Should something like the getWdot() be called in the updateProperties() method or one of the methods that this method calls?

Yes, this is the right approach. getWdot() is one of the cases that @speth warned about. Populating m_wdot prior to calculating residuals should avoid any issues (if you dig deeper: there shouldn't be any instance that requires setGas, which is the problematic call). If everything is done correctly, the order of execution should not matter.

Regarding your question where to update m_wdot; my recommendation would be to place it within the loop in StFlow::updateThermo:

    void updateThermo(const double* x, size_t j0, size_t j1) {
        for (size_t j = j0; j <= j1; j++) {
            setGas(x, j);
            m_rho[j] = m_thermo->density();
            m_wtm[j] = m_thermo->meanMolecularWeight();
            m_cp[j] = m_thermo->cp_mass();
            m_thermo->getPartialMolarEnthalpies(&m_hk(0, j));
            m_kin->getNetProductionRates(&m_wdot(0, j)); // <----- here
        }
    }

The key is to not call setGas more than absolutely necessary (this is where all the caching happens). The updateThermo name of the method may not be ideal, but it's the best place.

I have made those changes. I also had to make a change to the IonFlow class, as it used the removed evalResidual() method. The IonFlow class now overloads the evalSpecies() and evalElectricField() methods, and in the overloaded methods, it calls the base StFlow class' methods before extending them. The unit tests pass.

With regards to the diag variable name, I had originally interpreted it to mean something like the coefficient on the diagonal entry. Should we change the name for all the diag's to mask in this PR?

With regards to the diag variable name, I had originally interpreted it to mean something like the coefficient on the diagonal entry. Should we change the name for all the diag's to mask in this PR?

If you look at the docstring of Domain1D::eval, it is actually a "Boolean mask indicating whether each solution component has a time derivative (1) or not (0)"; since the PR touches most of the interfaces where it is used, I think it makes sense to rename.

With regards to the diag variable name, I had originally interpreted it to mean something like the coefficient on the diagonal entry. Should we change the name for all the diag's to mask in this PR?

If you look at the docstring of Domain1D::eval, it is actually a "Boolean mask indicating whether each solution component has a time derivative (1) or not (0)"; since the PR touches most of the interfaces where it is used, I think it makes sense to rename.

I'm 👎 on mask as a replacement for diag. I don't think describing this as a "mask" improves the ability to understand what it means.

I'm 👎 on mask as a replacement for diag. I don't think describing this as a "mask" improves the ability to understand what it means.

@speth ... ok. The proposed mask would resolve the inconsistent interface; while there may be a rationale for diag it is undocumented/unintuitive? I think it's an opportunity to make things more intuitive - are there better terms?

Edit: what about transient? ... probably not as this is ultimately a steady-state solution, even if pseudo-time stepping is involved

@wandadars ... thanks. The docstrings themselves look good, but Cantera's convention is to put them in the *.h files. Could you move them over? (leaving the code itself in the *.cpp files)

Regarding @speth's dislike of mask, it unfortunately appears that there isn't a consensus (I had assumed that using what is already documented at the external API would be a good solution). Unless we find a better name, it probably is best to leave things at the status quo, i.e. diag. Sorry about this ...

I'm 👎 on mask as a replacement for diag. I don't think describing this as a "mask" improves the ability to understand what it means.

@speth ... ok. The proposed mask would resolve the inconsistent interface; while there may be a rationale for diag it is undocumented/unintuitive? I think it's an opportunity to make things more intuitive - are there better terms?

Sorry for the terse response earlier. I recognize that the solver is very poorly documented at present, but I think that makes it even more important to keep certain things in place where they do actually have a connection to the mathematics behind the solver. In the transient mode, where the problem being solved is a DAE, diag is in fact a diagonal matrix that's used in constructing the Jacobian. Borrowing from the SUNDIALS IDAS docs:
$$J = \frac{\partial G}{\partial y} = \frac{\partial F}{\partial y} + \alpha\frac{\partial F}{\partial \dot{y}}$$
diag is the matrix $\partial F/\partial\dot{y}$, where $\alpha$ is roughly the inverse step size, i.e. what is called rdt in Cantera's solver.

@speth ... thank you for providing the context. This is definitely a decisive reason to stick with diag. We absolutely should add your explanation somewhere, although this is likely going beyond the scope of this PR. As mentioned above, I extrapolated from the interface/docstring provided in Domain1D::eval, which - while itself not being wrong - uses a less appropriate variable name.

I know this is not the topic of this PR, but one thing that I noticed in the continuity residual equation was the following:

} else if (m_isFree) {
                // terms involving V are zero as V=0 by definition
                if (grid(j) > m_zfixed) {
                    rsd[index(c_offset_U,j)] =
                        - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1];
                } else if (grid(j) == m_zfixed) {
                    if (m_do_energy[j]) {
                        rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
                    } else {
                        rsd[index(c_offset_U,j)] = (rho_u(x,j)
                                                    - m_rho[0]*0.3); // why 0.3?
                    }
                } else if (grid(j) < m_zfixed) {
                    rsd[index(c_offset_U,j)] =
                        - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j];
                }

In particular, the comparison of grid(j) == m_zfixed, a floating point comparison testing for equality. Is there some underlying situation that makes this a permissible operation? I ran into this issue with the 2-point continuation algorithm and I ended up just using an epsilon to test for floating point equality. Should this be done here too?

In particular, the comparison of grid(j) == m_zfixed, a floating point comparison testing for equality. Is there some underlying situation that makes this a permissible operation? I ran into this issue with the 2-point continuation algorithm and I ended up just using an epsilon to test for floating point equality. Should this be done here too?

It is absolutely true that exact floating point comparisons are problematic. My interpretation for why what you point to works / is legitimate is that grid contains exactly one point that is identical to m_zfixed (this goes back to how things are set up in Sim1D::setFixedTemperature), so the comparison checks out as it is indeed the same value. Whether or not this is the same situation as in the 2-point continuation case, I cannot say.

In particular, the comparison of grid(j) == m_zfixed, a floating point comparison testing for equality. Is there some underlying situation that makes this a permissible operation? ...

It is absolutely true that exact floating point comparisons are problematic. My interpretation for why what you point to works / is legitimate is that grid contains exactly one point that is identical to m_zfixed (this goes back to how things are set up in Sim1D::setFixedTemperature), so the comparison checks out as it is indeed the same value. ...

To elaborate, the reason an exact comparison is possible is because the value of m_zfixed was copied from the grid in the first place, not the result of a floating point computation that's subject to rounding errors. The reason for storing the value of T at the fixed point, rather than the index, is because the index changes when we regrid, while the value stays constant.

Thanks for the clarification @speth !

Once this PR is completed, I plan to open a PR for the two-point control method that will be build upon the changes in this PR.

The diag vector is set to zero here, before calling the eval method on each domain.

cantera/src/oneD/OneDim.cpp

Lines 246 to 255 in 9ee13eb

	void OneDim::eval(size_t j, double* x, double* r, double rdt, int count)
	{
	clock_t t0 = clock();
	if (m_interrupt) {
	m_interrupt->eval(m_nevals);
	}
	fill(r, r + m_size, 0.0);
	if (j == npos) {
	fill(m_mask.begin(), m_mask.end(), 0);
	}

Given that this default is the opposite of what is needed for most of the flame governing equations, meaning we have to have the assignments of diag[i] = 1, I think there's something to be said for being explicit and setting it to zero for the algebraic variables.

Thanks @speth. It's good to know that it is defaulted to zero. My original thought was that if activating the pseudo-time was the action that had to be explicit, then it would remove a lot of lines of code. I checked and we set the diag variable to 1 in 3 places in the entire stagnation flow code. Is it that the solver is expecting to see diag of 1 and we have to set it to 0 as that is the less likely situation?

Thank you @ischoegl . I appreciate you walking me through the Github PR process for this one.

@speth ... in order to add a deprecation process for StFlow::evalResidual (which is not part of the official API, but would still be the obvious entry point for derivative work), I'd suggest to rename the StFlow class to Flow (or FlowDomain) as it really isn't a stagnation flow in many cases anyways. Then, recreate a legacy StFlow with a depreciation warning triggered in the constructor, as well as overloaded StFlow::eval that calls a StFlow::evalResidual as before (plus anything else that is defunct), so the old interface is still available even if the internals are refactored. That way, we could cover the eventuality of anyone using derived classes out there in the wild. I'd be happy to add this as a follow-up PR, once this is merged.

I've been hitting the "resolve conversation" button on review items that I have made the requested changes for. Is it that the proper way to use that Github feature?

I've been hitting the "resolve conversation" button on review items that I have made the requested changes for. Is it that the proper way to use that Github feature?

Yes, that's how we use it, at least, for anything that's straightforward.

I missed one of @speth suggestions for the docstrings, so I just pushed the update. I think other than that I believe I have covered all the suggestions.