sot / xija

Thermal modeling framework for Chandra X-ray Observatory
https://sot.github.io/xija
BSD 3-Clause "New" or "Revised" License
9 stars 5 forks source link

New ODE solver #75

Closed jzuhone closed 4 years ago

jzuhone commented 4 years ago

This PR implements a new ODE solver for xija whilst leaving the old one as-is and as the default.

The new solver performs an Runge-Kutta integration on every timestep, using half-steps in between steps. This is in contrast to the current solver which performs the integration on every two steps and treats the step in between as the half-step. This requires us to interpolate quantities which go into forming the derivative to the half-step, which is simply done linearly here. By default, the integration is RK2, but RK4 is given as an option. Both of these options can be controlled by input into the XijaModel class or by setting them in the model specification file.

This solver will remove the issue of starting at an odd or even index since every step is now evaluated in the same fashion.

So far I have tested this by refitting the 1DEAMZT model with the new solver. Other ideas for tests are welcome. A PR for the new 1DEAMZT model is incoming. Also open to suggestions for better naming schemes.

Unrelated, I've also quieted down a FutureWarning which is associated with the use of numpy.hstack.

taldcroft commented 4 years ago

@jzuhone - I'm planning to take the rest of the week for a sprint on annie, but will have a look when I get back.

One quick question is on speed performance, e.g. time to make a prediction for a year of data.

jeanconn commented 4 years ago

PR reviewed and approved by TWG. Could probably go in 2020_205 matlab release to make the solver available for future thermal model updates via the expedited thermal model promotion procedure. Will plan to hold for review by @taldcroft

jskrist commented 4 years ago

Could probably go in 2020_205 matlab release

small clarification, the numbering for releases in 2020 will start with 2020_005 and increment by 5 for each major release (which leaves room for up to 4 "quick releases" between major releases).

jeanconn commented 4 years ago

Yep, my 2020_205 was a typo. Meant 2020_005.

jeanconn commented 4 years ago

Don't know yet if I need to care about the new warnings on build of calc_model_new this throws on CentOS5

ska3-presto: python setup.py build_ext --inplace
running build_ext
building 'xija.core_new' extension
gcc -pthread -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/proj/sot/ska3/flight/include/python3.6m -c xija/core_new.c -o build/temp.linux-x86_64-3.6/xija/core_new.o
xija/core_new.c: In function ‘calc_model_new’:
xija/core_new.c:167:19: warning: ‘k4’ may be used uninitialized in this function [-Wmaybe-uninitialized]
             k4[i] = dt * deriv[i];
                   ^
xija/core_new.c:156:19: warning: ‘k3’ may be used uninitialized in this function [-Wmaybe-uninitialized]
             k3[i] = dt * deriv[i];
                   ^
gcc -pthread -shared -L/proj/sot/ska3/flight/lib -Wl,-rpath=/proj/sot/ska3/flight/lib,--no-as-needed build/temp.linux-x86_64-3.6/xija/core_new.o -L/proj/sot/ska3/flight/lib -lpython3.6m -o build/lib.linux-x86_64-3.6/xija/core_new.cpython-36m-x86_64-linux-gnu.so
/usr/bin/ld: skipping incompatible /usr/lib/../lib/libpthread.so when searching for -lpthread
/usr/bin/ld: skipping incompatible /usr/lib/../lib/libpthread.a when searching for -lpthread
/usr/bin/ld: skipping incompatible /usr/lib/../lib/libc.so when searching for -lc
/usr/bin/ld: skipping incompatible /usr/lib/../lib/libc.a when searching for -lc
copying build/lib.linux-x86_64-3.6/xija/core.cpython-36m-x86_64-linux-gnu.so -> xija
copying build/lib.linux-x86_64-3.6/xija/core_new.cpython-36m-x86_64-linux-gnu.so -> xija
jzuhone commented 4 years ago

@jeanconn I can probably get rid of those

jzuhone commented 4 years ago

@jeanconn those warnings should be gone now.

jeanconn commented 4 years ago

Thanks @jzuhone. @taldcroft I can confirm that I diffed the ska_testr regress outputs for our runs of the thermal models and, as expected, don't see any diffs in the temperatures.dat files that use the previous solver. I think "core_new" perhaps isn't a great name, but I don't have a better one in mind.

I think we should merge and make a release.

jzuhone commented 4 years ago

@jeanconn I'm not excited about the name either, to be frank. Happy to brainstorm about alternatives.

jzuhone commented 4 years ago

@jeanconn @taldcroft just a process note--my feeling is that this should go in before PR #69, because either way some rebasing will be needed and I think that it will be better for this to go first.

taldcroft commented 4 years ago

For name core_rk4 seems an obvious choice.

jeanconn commented 4 years ago

I thought it added rk4 as an option but the bigger difference was the interpolated time?

jzuhone commented 4 years ago

@taldcroft the problem with that is that RK4 is optional and it's not even set by default--the real difference here is the way it handles the steps by treating each step the same and evaluating derivatives at half-steps. The default is RK2.

taldcroft commented 4 years ago

So core_2020 for the 2020-and-beyond enhancement. Or core_jayz as a silly hip ref to the J.Z. developer of the new core.

jskrist commented 4 years ago

my vote is for core_rk2, which is descriptive, and hopefully distinctive enough. I just noticed that there is a folder in this repo called alt_core, is that intended to be a location for alternative cores? the files in there seem to be named numerically (core2, core3, etc.) so perhaps core7 should be the new name?

jeanconn commented 4 years ago

I think core_rk2 is not distinct enough when core_new can do either rk4 or rk2, right? Or @jskrist are you suggesting we rename the previous core? I have no thoughts on core7 vs core_2020.

taldcroft commented 4 years ago

The current core is supposed to be RK2 (albeit with a special "feature"), so that's why I didn't go for core_rk2. The versions in alt_core are basically just WIPs or explorations that I tossed there for reference.

jskrist commented 4 years ago

It seems like the functionality to swap cores opens up the possibility to separate the various algorithms into distinct 'cores', so perhaps the naming convention could be something like core_[algorithm], but that would mean more work to separate the rk2 and rk4 implementations in the core_new core.

@jeanconn - I was originally thinking about a name for the new core, but if renaming the old core is on the table, then something like core_rk2_dubstep or core_rk2_ish 😋 work.

jzuhone commented 4 years ago

The evolve_method parameter is either 1 or 2. So maybe we rename the cores to core_1 and core_2?

jzuhone commented 4 years ago

Hearing no objections, I'm going to make the change I just suggested.

jzuhone commented 4 years ago

@jskrist I’d be hesitant to split the RK2/RK4 parts from the second core since there would be a lot of code duplication.

jzuhone commented 4 years ago

Ok, I've made this change and tested on my end and all looks good, but I will assume @jeanconn will want to test also.

jeanconn commented 4 years ago

This looks fine in building and ska_testr regress tests on my end. Last question.. I've forgotten, does core.pyx belong in dev or in alt_core with the others?

jeanconn commented 4 years ago

I'm holding the ska update for the matlab tools for this.