EPM ODE state is lagged by one time step

[lrobion]

In Integrate.cpp for the original EPM, Original::RunMicrophysics tracks its state via a mix of a StateObserver object and the state vector x. Time evolution of the EPM is done two ways: an outer loop of log spaced timesteps with 300 values (hard coded) and an inner loop: from $t_i$ to $t_{i+1}$ of the outer loop the internal evolution is done by an adaptive integration using a Boost ODE solver. We do not control the timesteps there intentionally and let the ODE solver decide.

The observer's operator() gets called once per inner timestep to record the state during the inner steps. This is its implementation:

APCEMM/Code.v05-00/src/EPM/Models/Original/StateObserver.cpp

Lines 50 to 56 in adf714b

	void StateObserver::operator()(const Vector_1D &x, double t) {
	if ( ( m_count % m_write_every ) == 0 ) {
	m_states.push_back( x );
	m_times.push_back( t );
	m_count++;
	}
	}

m_write_every = 2 is hard coded to record every other inner timestep. m_count++; being inside that if check means that only the first inner timestep (the outer timestep $t_i$) is recorded. This is because m_count is 0 after the first inner timestep -> enter the if branch, record the initial state of the inner solve which is the outer timestep $t_i$), then increment m_count. Now m_count = 1 and the if check never passes for all subsequent timesteps, so the other inner steps are never recorded.

As it's written m_states (which really is just a reference to obs_Var owned by RunMicrophysics) is a vector of size 300 (the hard coded number of outer timesteps) whereas the intention with the state observer seems for it to be a vector of size N where N is the total number of inner steps taken to go from $t_0$ to $t_{max}$

This also means that the state in m_states[m_states.size() -1] is always an outer timestep behind: you'd think it's $t_{i+1}$ because that's the last element of the state array after running the inner ODE solve, except it's not it's just the inital state of the inner solve $t_i$. This means that the physics calculation that happen in a timestep run with a mix of data from the current and previous timestep, see this snippet right after the adaptive solve:

APCEMM/Code.v05-00/src/EPM/Models/Original/Integrate.cpp

Lines 319 to 324 in adf714b

	SO4l = observer.m_states[observer.m_states.size()-1][EPM_ind_SO4l] ;
	n_air = Na * x[EPM_ind_P]/(R * x[EPM_ind_T] * 1.0e6);
	n_air_prev = Na * P_b/(R * T_b * 1.0e6); // P_b and T_b seems to be used uninitialized during the first timestep with iTime = 0
	nPDF_new = ( SO4l*n_air - SO4l_b*n_air_prev);

nPDF_new = ( SO4l*n_air - SO4l_b*n_air_prev); is really nPDF_new = ( SO4l(t_i)*n_air(t_i+1) - SO4l_b(t_i-1)*n_air_prev(t_i-1); when what we want is $t_{i+1} - t_i$. This happens because x is truly the final state of the inner solve so that one is correctly at $t_{i+1}$

A side note is that using StateObserver.m_count or StateObserver::getLastElement() would also fail because StateObserver is passed by value to the adaptive solver so the increment to m_count are done on a copy which never resurfaces to the StateObserver that is used in runMicrophysics.

I imagine this has more repercussions on rest the EPM physics.