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

Adding a gradient algorithm #694

Open
wants to merge 14 commits into
base: develop
Choose a base branch
from
Open

Adding a gradient algorithm #694

wants to merge 14 commits into from

Conversation

taobrienlbl
Copy link
Collaborator

No description provided.

@taobrienlbl taobrienlbl requested review from burlen and amandasd May 17, 2022 23:49
@taobrienlbl taobrienlbl changed the title A new algorithm for calculating the gradient of a scalar field Adding a gradient algorithm May 17, 2022
Copy link
Collaborator

@burlen burlen left a comment

Choose a reason for hiding this comment

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

Hi Travis, The calculation looks correct for a uniform Cartesian mesh. I found a few metadata related issues that should be fixed.

alg/teca_gradient.cxx Show resolved Hide resolved
alg/teca_gradient.cxx Outdated Show resolved Hide resolved
alg/teca_gradient.cxx Outdated Show resolved Hide resolved
alg/teca_gradient.cxx Show resolved Hide resolved
@taobrienlbl
Copy link
Collaborator Author

I believe I've addressed all your comments @burlen

@amandasd - do you have any?

num_t *dest = grad_y;
num_t *src = grad_y + n_lon;
for (unsigned long i = 0; i < n_lon; ++i)
dest[i] = src[i+n_lon];
Copy link
Collaborator

@amandasd amandasd May 18, 2022

Choose a reason for hiding this comment

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

Wouldn't it be dest[i] = src[i] instead of dest[i] = src[i+n_lon] in line 130?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I see what you're saying. It does look like that, but I'm confused because this code reproduced a python version of the same boundary condition (implemented totally differently, so it's unlikely the same mistake was made twice).

@burlen, this is based on teca_vorticity.cxx (see line 117), so if there's a problem here, there is a problem there too. I always get turned around when thinking about pointer math though. What do you think?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I can confirm that this was a pointer math bug. (And it turns out there was a bug in the code I used to compare with the python-generated version, such that I wasn't actually comparing TECA's gradient at all.)

I've also submitted a fix to teca_vorticity.cxx.

Copy link
Collaborator

Choose a reason for hiding this comment

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

@taobrienlbl sorry for late reply, agree. Good catch @amandasd!

num_t dx = delta_x[j];

gx[0] = (f_x2[0] - f_x0[0]) / dx;
gy[0] = (f_y2[0] - f_y0[0]) / dx;
Copy link
Collaborator

Choose a reason for hiding this comment

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

Wouldn't it be dy instead of dx in line 92?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

yes, thank you on this and the one below!! The test wouldn't have caught this because dx=dy

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

this has been fixed, thanks!

num_t dx = delta_x[j];

gx[0] = (f_x2[0] - f_x0[0]) / dx;
gy[0] = (f_y2[0] - f_y0[0]) / dx;
Copy link
Collaborator

@amandasd amandasd May 18, 2022

Choose a reason for hiding this comment

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

Wouldn't it be dy instead of dx in line 107?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

yes!

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

this has been fixed, thanks!

unsigned long n_lat, bool periodic_lon=true)
{
size_t n_bytes = n_lat*sizeof(num_t);
num_t *delta_x = static_cast<num_t*>(malloc(n_bytes));
Copy link
Collaborator

Choose a reason for hiding this comment

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

Isn't better to use std::vector<num_t> instead of explicitly malloc?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

That would be more in the C++ style, yes. I copied the approach from teca_vorticity.cxx, so I'm wondering if there's a performance reason for using a dynamically allocated array rather than a vector: perhaps the overhead from bounds-checking in a vector? @burlen - do you have an opinion here on which approach should be used?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I ended up going with std::vector since the [] operator doesn't have any overheads.

Copy link
Collaborator

@burlen burlen May 19, 2022

Choose a reason for hiding this comment

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

@taobrienlbl malloc is preferable. std::vector<T> makes an extra pass over the array initializing all elements to T(). For POD types this zeros the values out. This is wasted cycles because you are immediately initializing after. There's also substantially more overhead when resizing the std::vector vs. using realloc, since the std::vector will always allocate a new buffer and copy the contents where as realloc can extend the memory region in place. The latter doesn't impact you here.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

OK, I'll revert it.

Copy link
Collaborator

Choose a reason for hiding this comment

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

@taobrienlbl It's fine to use vector here, I didn't mean for you to have to revert this, just wanted to give an answer as to why I use malloc.

@taobrienlbl taobrienlbl force-pushed the teca_gradient branch 2 times, most recently from 5fc44f9 to c82afe3 Compare May 19, 2022 19:20
@taobrienlbl
Copy link
Collaborator Author

taobrienlbl commented May 19, 2022 via email

@taobrienlbl
Copy link
Collaborator Author

@burlen - when and if you think this is ready, would you please approve for merging?

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

Successfully merging this pull request may close these issues.

3 participants