Oscillations like fortran #42
Conversation
…e. until 133.854 ka at which point it will give a RuntimeError: Time step below 1e-10.
…on dataclass. We want to use the scipy solver, i.e. solve_ivp, since it provides for implicit methods. Also we want to integrate over 250_000 years, so we should set the progress tracker interval accordingly.
…tions 35 from L'Heureux, but py-pde demands left and right boundary conditions. This means that 'curvature : 0' is a dummy boundary condition. We will make sure that it does not leak into our integration, by using only backward differencing for the spatial derivatives in the right-hand sides of the time derivative equations for CA and CC." This reverts commit 0d8191f.
…does indeed yield oscillations within 250 ka. Note that this commit includes a reversal of commit 0d8191f meaning that the 'dummy' bottom boundary conditions are reintroduced: self.bc_CA = [{'value': CA0}, {'curvature' : 0}]; self.bc_CC = [{'value': CC0}, {'curvature': 0}].
… a means of imposing boundary conditions. This seems a way to align with rhythmite when U changes sign, i.e. becomes negative (upward flow) which requires a bottom boundary condition for Araganonite and Calcite if one wants to keep applying the upwind solution. That bottom boundary condition is not available. The way rhythmite solves this is to apply a downwind solution, just for the deepest node, while all the other nodes retain upwind. Effectively this means that the spatial gradients for Aragonite and Calcite for the deepest node are equal to that gradient for the second deepest node. This is what I am trying to mimic here using 'adjacent_value'. Anyway, with the use of the 'typo' for the porosity diffusion coefficient this results in oscillations within 250 ka.
… boundary condition for Aragonite and Calcite. This reverts commit 73b91b4.
… a flat profile on aragonite and calcite at the bottom instead of a second spatial derivative=0 or a a first spatial derivative=adjacent_value.
…ss, added to the Fortran code, that inspired parameters.py. They are redundant for marlpde, so they were removed. Secondly, the dimension of the grid was doubled to check its possible effect on oscillations while retaining the depth resolution.
…ndition for aragonite and calcite makes oscilations vanish, similar to the slope=0 ad-hoc bottom boundary condition.
…nditions and porosity diffusion coefficient - that will yield oscillations independent of the grid size, retaining grid resolution.
…ed to an IDE is personal and should not be part of this repo.
…It includes removing the typo; any mainpulation of the porosity diffusion coefficient should be done through modifying b. Canonical values of phi0 and phi00 are now 0.8 and 0.8, respectively. These differ from Scenario A, but we are looking for oscillations. We will investigating larger grids and a higher grid resolution, but we will set for 5m and 2.5cm as default values. The default ad-hoc bottom boundary conditions for aragonite and calcite were chosen to align with rhythmite. This means that the spatial derivatives of the deepest and second deepest node are equal. If two adjacent spatial derivatives are equal, this means that the second spatial derivative is equal to 0. This is reflected by {'curvature': 0}.
|
Before getting to the discussion on F-V, I want to comment on this:
The "default value" of 5e-4 is clearly assigned so as to compensate for the typo. It has no other motivation. If we do not use the typo and use the correct expression for the porosity diffusion coefficient, this value has no meaning. I believe we can leave it there though as the choice of |
…TkAgg' which requires the 'tk' interactive framework, as 'headless' is currently running. This should fix that.
…ntegrating-diagenetic-equations-using-Python into Oscillations_like_Fortran
… It is set here equal to T* to make a comparison with ground truth data - needed for the regression tests - easier. But obviously, to see any oscillations, one will need to set it to somthing (much) larger.
…ve been changed, to facilitate the search for oscillations.
|
Confusing indeed!
What I am struggling with is that in Which will launch py-pde's native explicit solver with a Runge-Kutta scheme. However, Somewhat related: it may be good to add a line to the |
|
I think what you proposed is indeed the best solution. But it would be useful to include a warning somewhere if anyone (including our future selves) tries to change the solver. E.g. I already encountered that confusion that changing the
|
Sorry, not entirely correct, |
…ike 'backend' and 'ret_info' will always be propagated, but 'scheme' and 'adaptive'only for py-pde's native 'explicit' solver.
|
I ran after this change with the following options: and got probably not relevant but I also set |
…. Running marlpde as from the previous commit led to 'RuntimeError: Encountered Not-A-Number (NaN) in evolution' (this is with 'LSODA'). Stability was added by providing the equivalent of a Jacobian sparsity matrix for 'LSODA' which is done by setting 'lband' and 'uband'. This not only provided stability but also led to a substantial performance increase!
|
Yeah, |
|
But what should be set if I want to use |
|
You got me there... |
…ny solver other than 'scipy' is chosen: it is not a keyword argument in that case. This commit will fix that with a conditional --> two separate calls. Not the most elegant soluttion, it would be nicer to apply some filtering of the arguments, or by using a ternary operator, but I did not find a way to get that working for keyword arguments; for positional arguments a ternary operator should work.
|
Pls try again. |
|
I found a more elegant solution on SO, using the Btw, we already apply |
|
Let me try to fix that: it would be nice if Figure 3 would be reproduced with our default settings. |
Let me just add that commit. |
O wait, that is perhaps as it should be, since we now have 0.8 as initial and surface default porosity. |
|
You changed the phis to 0.8 because we wanted potential users to find the conditions that lead the oscillations, as this is the point of interest of this model. Fig. 3 is more trivial, so to say. I guess if anyone uses this model, it is primarily to look at the oscillations. But it boils down to what the documentation would be. Have you thought of what the documentation will be? Our SMP has this field empty. I guess a simple documentation generated from your excellent docstrings with some preamble will be sufficient, as we can refer to the paper on the model itself. |
|
Sounds like a plan! |
|
Then I think this branch can be merged and I'll assign you to another issue for the documentation |
…cationWarning: Jupyter is migrating its paths to use standard platformdirs given by the platformdirs library. To remove this warning and see the appropriate new directories, set the environment variable and then run .' To execute the second command, 'jupyter', as an executable, was needed, but it was not available, so I pip installed it.
…ices for 'solver', 'method', 'scheme' and 'adaptive'. Actually, it is about more than incompatible choices, it is about certain Solver() attributes that should be complete removed before they are passed as an argument to the eq.solve(...) call in Evolve_Scenario.py. E.g. for the 'explicit' solver, one cannot also provide a 'method', while this is possible for the 'scipy' solver. In a previous commit, a conditional was added before the eq.solve call, but that was considered not the most elegant solution; filtering out attributes in a __post_init__ method of the Solver dataclass seems more sustainable, since it can handle increasing complexity and that is what is achieved through this commit. Also some better formatting in the dataclasses (removing spaces). Next, the dict from asdict(Solver()) is passed 'as is' to the eq.solve call, where we previously had End_time = kwargs[tmax]/Tstar before the eq.solve(state, t_range=End_time, ...) call. That is why we no longer have in tmax in years in the Solver dataclass, but t_range: int = 1 instead, i.e. a default value of 1 times T*. Likewise retinfo--> ret_info (since we no longer have the 'conversion' of keyword names in the eq.solve call. Last but not least, the number of cells (nodes), denoted by N, is not an eq.solve argument and therefore had to be removed from the Solver dataclass and added to the Map_Scenario method, more specifically the 'post_init' therein.
…e 'method' atrribute from the Solver dataclass has been filtered out in its __post_init__, e.g. because the user has chosen solver: str = 'explicit', then we can no longer have a method=something, not even method=None keyword argument in the eq.solve call, unless we keep the conditional we have around the eq.solve call. However, that is not a sustainable solution with increasing incompatibility issues within the Solver dataclass. Instead we want an eq.solve call like this: eq.solve(state, **solver_parms, tracker_parms, liveplots, ...). That requires not concatenating all arguments of integrate_equations in the single 'kwargs' dict, but keeping 'solver_parms', 'tracking_parms' and 'pde_parms' separate. But we do want to record all metadata in the stored hdf5 file, that is why we still need a concatenation 'storage_parms=solver_parms | tracker_parms | pde_parms' and 'storage = FileStorage(stored_results, info=storage_parms)'.
… instead of kwargs, i.e. one (concatenated) dict, we need to adjust the regression tests to handle this.
…' implicit methods from solve_ivp, but not working, no clue why. The integration does not progress in time, but does not crash. The implementation should be identical to commit 2197188 from the Use_solve_ivp_without_py-pde_wrapper branch. In that branch this Jacobian sparsity matrix works succesfully. Lines 'jac_sparsity: csr_matrix = None' in the Solver dataclass and 'self.jac_sparsity = jacobian_sparsity()' in its __post_init__ method removed, to avoid users applying this Jacobian sparsity matrix implementation for now.
Fixes #39 and #40.
The 'typo' is not included here, any manipulation of the porosity diffusion coefficient should be done by modifying
b.However,
bis set here to a default value of 5e-4, not to the rhythmite value. I.e. after merging withmainand then running, marlpde will not yield oscillations within 250 ka.Initial and surface porosity are both set to 0.8, i.e. different from Scenario A.
There is an ongoing discussion on how to apply the FV scheme. It would be neat to settle that discussion asap. If the conclusion would be that we interpret the FV (1977) paper in such a way that this allows for the same spatial derivative to be calculated in two different ways for the same rhs, then an additional commit is required that would align marlpde with rhythmite in this respect. It would be nice if we could add that commit before this PR is merged. The FV paper is not an easy paper to interpret, though.