It seems that the methods isotherm_reinit() are different only at computation of scale_aqua and scale_sorbed. So it could be moved to SorptionDual and the only method which would be virtual would be compute_sorbing_scale(). It is prepared in comment code.
START_GLOBAL_TIMER(tag) - this calls the start_timer, which creates local timer on the correct place in the hierarchy, further this timer is added to the list of global timers, this contains groups of timers with same tag, and collect/sum data from these timers in the report.
Implement co-located vector class that takes reference of the original Vector and is of the same size. FullIterator should have method for accessing elements in the original vector.
in create_transport_matric_mpi, there there is condition edge_flow > ZERO this makes matrix sparser, but can lead to elements without outflow and other problems when there are big differences in fluxes, more over it doesn't work if overall flow is very small