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.

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

Refiner documentation #1786

Open
wants to merge 29 commits into
base: main
Choose a base branch
from
Open

Conversation

wandadars
Copy link
Contributor

This is a WIP

Mostly documentation changes, with a few variable name adjustments for clarity.

Checklist

  • The pull request includes a clear description of this code change
  • Commit messages have short titles and reference relevant issues
  • Build passes (scons build & scons test) and unit tests address code coverage
  • Style & formatting of contributed code follows contributing guidelines
  • The pull request is ready for review

@wandadars
Copy link
Contributor Author

wandadars commented Aug 22, 2024

This getNewGrid method doesn't seem to be used anywhere in the refiner or elsewhere.

Also, this failing test is a bit perplexing to me as it doesn't seem that I actually changed anything.

@wandadars
Copy link
Contributor Author

One strange thing about the way that certain points get marked as being in the "keep" list is that the curvature check will only keep the point next to the point of interest if it decides to split the two intervals to the right of the point of interest. This is different from how the slope criterion check marks points. I haven't quite grasped any underlying logic for which points get added to the "keep" list during each refinement check.

@wandadars
Copy link
Contributor Author

@speth I went ahead of created a new page that talks about the grid refinement.

Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

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

@wandadars ... thank you for taking this on! While this is still work-in-progress, I wanted to provide some feedback that may be a little nit-picky, but should help with improving consistency of our documentation as well as our code base.

Edit: I am not sure why unit tests broke, but it does appear that you inadvertently changed criteria for regridding, which can hopefully be rolled back.

vector<bool> m_active;

//! Refinement criteria
Copy link
Member

Choose a reason for hiding this comment

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

This docstring only applies to the first of the subsequent member variables.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ah, my vs-code editor was leading me astray there, it applied it to all the subsequent items.

Copy link
Member

Choose a reason for hiding this comment

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

If this is the case, there may be a setting I’m unaware of. The doxygen output will tell.

include/cantera/oneD/refine.h Outdated Show resolved Hide resolved
size_t m_nv;
size_t m_npmax = 1000;

Domain1D* m_domain; //! Pointer to the domain to be refined
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
Domain1D* m_domain; //! Pointer to the domain to be refined
Domain1D* m_domain; //!< Pointer to the domain to be refined.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

We want two spaces between the semicolon and the start of the comment?

Copy link
Member

Choose a reason for hiding this comment

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

I don’t think we have a guideline on this. On the Python side, it’s a recommendation that I have seen, and I think it’s slightly more readable. It’s really your choice; the < was definitely missing.

//! Names of components that require the addition of new grid points
set<string> m_c;
set<string> m_component_name;
Copy link
Member

@ischoegl ischoegl Sep 1, 2024

Choose a reason for hiding this comment

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

Suggested change
set<string> m_component_name;
set<string> m_componentNames;

See C++ Style Guide. I believe that for variable names that do not appear as part of the public interface we are not as concerned about descriptive names, which have to be balanced with readability of the implementation itself. m_c can certainly approved upon, but m_cNames would be a shorter compromise.

set<size_t> m_loc;
set<size_t> m_insertion_points;

//! Map of grid point indices that should be kept, 1=keep, -1=remove, 0=unset
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
//! Map of grid point indices that should be kept, 1=keep, -1=remove, 0=unset
//! Map of grid point indices that should be kept. Values are ...

Detailed description should not be part of the first statement, which forms Doxygen's @brief description. Also, which integer is what?

As an aside, using an enum as the value instead of these sentinel values would further clarify things.

include/cantera/oneD/refine.h Outdated Show resolved Hide resolved
src/oneD/refine.cpp Outdated Show resolved Hide resolved
double prune() {
return m_prune;
}

protected:
//! Indices of grid points that need new grid points added after them
set<size_t> m_loc;
set<size_t> m_insertion_points;
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
set<size_t> m_insertion_points;
set<size_t> m_insertionPoints;

I am torn about changing member variable names here: m_loc is used in many other contexts within the Cantera code base where the meaning is well known. At the same time, it is a more meaningful name.

While we aren't completely consistent on this, the default naming style would be to use camelCase after the initial m_, which in this case would yield m_insertionPoints - see C++ Style Guidelines

Copy link
Contributor Author

Choose a reason for hiding this comment

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

How about m_insertPts ?

@ischoegl
Copy link
Member

ischoegl commented Sep 1, 2024

This getNewGrid method doesn't seem to be used anywhere in the refiner or elsewhere.

If it's unused, please go ahead and deprecate it. For this, we usually issue a warn_deprecated warning, and mark the method with Doxygen's @deprecated tag. Search for these and you will find instances which can serve as a template.

Also, this failing test is a bit perplexing to me as it doesn't seem that I actually changed anything.

I disagree - some of the failing tests show

E   *******************************************************************************
E   CanteraError thrown by Refiner::analyze:
E   max number of grid points reached (1000).
E   -------------------------------------------------------------------------------

which indicates that regridding was inadvertently changed. Same goes for TestImpingingJet.test_restore_legacy_hdf et al., where differences in array length are explained by regridding that should not take place as far as I remember.

@wandadars
Copy link
Contributor Author

Thank you for the feedback @ischoegl . I'm still looking into how exactly I changed the algorithm.

@wandadars
Copy link
Contributor Author

The source of the change that was causing the tests to fail was the change of:

double ratio = fabs(slope[j+1] - slope[j]) / (max_change + m_thresh/dz[j]);
to
double ratio = fabs(slope[j+1] - slope[j]) / max_change;

Where I had just moved the m_thresh up into the max_change variable, as the purpose of that variable seems to be a protection against divide-by-zero FPEs.

double max_change = m_curve*(slopeMax - slopeMin) + m_thresh;

The thing that I'm not sure about is why this mattered. There's the division of m_thresh by dz[j], which doesn't make much sense to me given that we're just protecting the denominator from going to zero.

m_thresh is roughly 1.49011745e-8. For cases where the grid spacing is close, then the division by the dz term would make the effective FPE threshold larger.

The failing test is shown below. The width of the domain is 0.05m, and if there's say 100 points, then the dz is on the order of 5e-4m, which would make the effective FPE threshold something like 2.98e-5. Doing it the way I had adjusted it to, the threshold would still be 1.4e-8 for guarding against the FPEs.

def test_unity_lewis(self):
        self.create_sim(ct.one_atm, 300, 'H2:1.1, O2:1, AR:5.3')
        self.sim.transport_model = 'unity-Lewis-number'
        self.sim.set_refine_criteria(ratio=3.0, slope=0.08, curve=0.12)
        self.sim.solve(loglevel=0, auto=True)
        dh_unity_lewis = np.ptp(self.sim.enthalpy_mass)

        self.sim.transport_model = 'mixture-averaged'
        self.sim.solve(loglevel=0)
        dh_mix = np.ptp(self.sim.enthalpy_mass)

        # deviation of enthalpy should be much lower for unity Le model (tends
        # towards zero as grid is refined)
        self.assertLess(dh_unity_lewis, 0.1 * dh_mix)

@wandadars
Copy link
Contributor Author

wandadars commented Sep 18, 2024

I put some print statements into that part of the code to compare the two, and it does seem like my change makes the ratio larger than the other way. This makes sense because for a smaller value of the max_change, the m_thresh by itself would be smaller. I don't really see the logic of making the FPE guard larger though.

Edit: I'm inclined to think that the old way of handling the ratio criterion was permitting un-resolved points to pass erroneously.

Old demoninator: 0.00191138
New demoninator: 4.04799e-06
Old Ratio: 0.000436741
New Ratio: 0.20622

Old demoninator: 0.00191138
New demoninator: 4.04799e-06
Old Ratio: 0.000578154
New Ratio: 0.272993

Old demoninator: 0.00191138
New demoninator: 4.04799e-06
Old Ratio: 0.000713674
New Ratio: 0.336983

Old demoninator: 0.00191138
New demoninator: 4.04799e-06
Old Ratio: 0.000835034
New Ratio: 0.394286

Old demoninator: 0.00191138
New demoninator: 4.04799e-06
Old Ratio: 0.000930045
New Ratio: 0.439149

Old demoninator: 0.00191138
New demoninator: 4.04799e-06
Old Ratio: 0.000978374
New Ratio: 0.461968

Old demoninator: 0.00191138
New demoninator: 4.04799e-06
Old Ratio: 0.00211619
New Ratio: 0.999222

Old demoninator: 0.000957707
New demoninator: 4.04799e-06
Old Ratio: 0.00484523
New Ratio: 1.14632

Old demoninator: 0.000957707
New demoninator: 4.04799e-06
Old Ratio: 0.00504515
New Ratio: 1.19362

@wandadars
Copy link
Contributor Author

wandadars commented Sep 18, 2024

I looked at that failing test and it does look like for the case when the denominator is simply using the m_thresh, a ton of points get put in a region of almost constant slope. So maybe the guard against regions of constant slope isn't working properly causing the algorithm to be activated?

I printed out all the cases where the ratio was greater than 1 using the newer ratio definition and it looks like there's a sensitivity going on where maybe the more points that are added causes noise that can make the slope vary significantly, which drives more points to be added.

run.log

@wandadars
Copy link
Contributor Author

wandadars commented Sep 19, 2024

I think one issue was that the refiner was considering (for the failing test) the species Argon, which was a constant, but there were small fluctuations in the value, and because the slope criteria only tested for regions of constant slope, it entered into the slope refinement section and inserted points, which seems to have generated more noise in the slope, which generated more points. If I have understood the nature of the issue, basically the zero slope case seems to have had positive/negative values or just noise that the first if-statement catch seems to have deemed as being significant enough of a variation from a constant to let it enter the clause. Maybe a case where because the value under consideration was zero, this is different from a case of a non-zero constant slope.

I put an extra clause on the slope constant-value-check clause where it just uses the previous constant-value-clause for the previous if-statement that refines based on values. Doing this removed the excessive refinement ( it was putting 3500 points in the free flame case that was failing). It returned back to putting around 100 points over the flame after adding this additional check.

Looking at Argon, we have the following outputs for the free-flame test case.

Component: AR
valMax-valMin: 2.80705e-11
min_range*valMagnitude: 0.00860884
slopeMax-slopeMin: 2.79602e-08
min_range*slopeMagnitude: 2.71097e-10
Slope[j]: 2.25352e-08
Slope[j+1]: 9.84457e-10
SlopeMin and SlopeMax: -8.5052e-10 2.71097e-08
Old demoninator: 4.76871e-05
New demoninator: 1.82564e-08
Old Ratio: 0.00045192
New Ratio: 1.18045

Looking at the valMax-ValMin relative to the min_range*valMagnitude, it is clear this is a constant value. But looking at the slope quantity tells a different story because the magnitude of the slope is close to zero, so the noise is considered to be significant. We can also see the old ratio and new ratio have very different opinions on the ratio. The old ratio had that m_thresh term, which added just enough to the denominator so as to ignore slopes that were very close it seems. I don't think that was the intention of the original coders, as they also used the m_thresh in the previous check for the value refinement. I may be incorrect though. It seems like the expected behavior is to not let the "bad" points even get into the refinement insertion algorithm instead of letting them in and then within the algorithm ignoring ones that fail to meet that m_thresh/dz[j] quantity.

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

Successfully merging this pull request may close these issues.

3 participants