Use a piecewise cubic polynomial for charateristic charge in GIEX adsorption model by sleweke · Pull Request #97 · cadet/CADET-Core

sleweke · 2021-07-21T00:10:15Z

The characteristic charge nu is extended to a piecewise cubic polynomial in the GIEX adsorption model. The nu of each component can have different polynomial coefficients and breaks, but their amount must be the same over all components.

If the breaks are not provided, a single global polynomial is applied. This maintains backwards compatibility.

The code compiles but is untested. Only the file format docs have been added.

schmoelder · 2021-07-21T14:28:30Z

Thanks for putting this together!

I pulled and it compiled fine. However, if I try to run it, it does not seem to read the parameters properly.

ERROR: GIEX_NU and GIEX_NU_LIN do not have the same length

I went in and printed the values and I got 120 for GIEX_NU, and 0 for GIEX_NU_LIN although I did provide the parameters correctly (see attached file). Maybe it's because of the 'skipConfig' flag?
https://github.com/modsim/CADET/blob/701166b0eceb63e359be229244b370e6465ecfe9/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp#L47

test_giex.zip

sleweke · 2021-07-21T23:06:42Z

There were some bugs (mostly copy & paste), as well as a conceptual problem:

Maybe it's because of the 'skipConfig' flag?

Indeed, that was the case. The validation code ran before the custom parameter reading code.

The provided test file has GIEX_NU_LIN, GIEX_NU_QUAD, GIEX_NU_CUBE, and GIEX_NU_BREAKS all 0.0. This cannot work!

schmoelder · 2021-07-22T07:29:03Z

+	// Check breaks
+	for (int i = 0; i < nComp; ++i)
+	{
+		cadet::active const* const b = _nuBreaks.get().data() + nPieces * i;


Since there are nPieces + 1 breaks for each component, I think this line should read

cadet::active const* const b = _nuBreaks.get().data() + (nPieces + 1) * i;

Otherwise the exception is thrown.

Well spotted!

schmoelder · 2021-07-22T07:31:27Z

-**Type:** double     **Length:** NCOMP
-===================  ========================= 
+===================  ===============================  
+**Type:** double     **Length:** multiples of NCOMP


Is this really NCOMP or rather NBOUND?

It's NCOMP. Each component has at most 1 bound state. Entries of non-binding components are ignored (but must be present just in the other adsorption models).

Ok, good! Otherwise setting up the config would be a nightmare! ^^

schmoelder · 2021-07-22T07:38:12Z

Maybe it's because of the 'skipConfig' flag?

Indeed, that was the case. The validation code ran before the custom parameter reading code.

Lucky guess! ^^

The provided test file has GIEX_NU_LIN, GIEX_NU_QUAD, GIEX_NU_CUBE, and GIEX_NU_BREAKS all 0.0. This cannot work!

Oh, I sent you the wrong file... Here is a new one that runs. However, I'm not sure yet that it works properly.

test_giex.zip

Here is another one where I tried to use the old interface but it does not run (returncode=-11).
test_giex_orig.zip

schmoelder · 2021-07-22T13:25:42Z

I noticed that the total resin capacity seems to be exceeded sometimes.

giex_check_lambda.zip

For example, if you check cadet.root.output.solution.unit_003.solution_solid[-1,0,:] (last solution time, first column cell, all components), it returns

array([3000., 0., 1026.06197414, 34.4423779, 34987893])

even though lambda is only 3000. I don't seem to have the same issue if I use only one global polynomial.

sleweke · 2021-07-22T14:18:58Z

Splines look really weird:

Note that the formula for the j-th spline piece is:

s_j(t) = const_j + (t - breaks_j) * (lin_j + (t - breaks_j) * (quad_j + (t - breaks_j) * cube_j))

Maybe the visualization code is wrong:

#!/usr/bin/env python3

import h5py
import matplotlib.pyplot as plt
import numpy as np

with h5py.File('giex_check_lambda.h5', 'r') as f:
	nu_const = f['/input/model/unit_003/adsorption/GIEX_NU'][()]
	nu_lin = f['/input/model/unit_003/adsorption/GIEX_NU_LIN'][()]
	nu_quad = f['/input/model/unit_003/adsorption/GIEX_NU_QUAD'][()]
	nu_cube = f['/input/model/unit_003/adsorption/GIEX_NU_CUBE'][()]
	nu_breaks = f['/input/model/unit_003/adsorption/GIEX_NU_BREAKS'][()]
	n_comp = f['/input/model/unit_003/NCOMP'][()]

n_pieces = int(len(nu_const) / n_comp)
print('Comp {} Pieces {}'.format(n_comp, n_pieces))
n_per_piece = 100

fig, axs = plt.subplots(2, 3)
axs = axs.reshape((6,))

for i in range(n_comp):
	breaks = nu_breaks[(n_pieces+1)*i:(n_pieces+1)*(i+1)]
	con = nu_const[n_pieces * i:n_pieces*(i+1)]
	lin = nu_const[n_pieces * i:n_pieces*(i+1)]
	quad = nu_const[n_pieces * i:n_pieces*(i+1)]
	cube = nu_const[n_pieces * i:n_pieces*(i+1)]
	
	x = np.empty((n_per_piece * n_pieces,))
	y = np.empty((n_per_piece * n_pieces,))
	for j in range(n_pieces):
		x[n_per_piece * j: n_per_piece * (j+1)] = np.linspace(breaks[j], breaks[j+1], n_per_piece)
		t = x[n_per_piece * j: n_per_piece * (j+1)] - breaks[j]
		y[n_per_piece * j: n_per_piece * (j+1)] = con[j] + t * (lin[j] + t * (quad[j] + t * cube[j]))

	axs[i].plot(x,y)
	axs[i].set_title('Component {}'.format(i))

plt.show()

schmoelder · 2021-07-22T14:45:49Z

I think there is an issue with the variable names here:

	con = nu_const[n_pieces * i:n_pieces*(i+1)]
	lin = nu_const[n_pieces * i:n_pieces*(i+1)]
	quad = nu_const[n_pieces * i:n_pieces*(i+1)]
	cube = nu_const[n_pieces * i:n_pieces*(i+1)]

If you update to:

con = nu_const[n_pieces * i:n_pieces*(i+1)]
lin = nu_lin[n_pieces * i:n_pieces*(i+1)]
quad = nu_quad[n_pieces * i:n_pieces*(i+1)]
cube = nu_cube[n_pieces * i:n_pieces*(i+1)]

It seems alright.

lieres · 2021-07-26T10:53:40Z

I assume most if not all applications will be based on cubic splines. However, the user needs to take care of the continuity and boundary conditions at the nodes, which can be tedious and error prone. In addition, it complicates parameter estimation from measurement data. Have you considered using a different parameterization, e.g. a series of x and y values of the interpolated function with natural boundary conditions?

schmoelder · 2021-07-29T10:16:41Z

Have you considered using a different parameterization, e.g. a series of x and y values of the interpolated function with natural boundary conditions?

As far as I know, this was not considered. However, to get the polynomial coefficients, I simply used scipy's cubic splines module which are already twice continuously differentiable.

Nevertheless, I still can't get it to run properly. It gets stuck shortly after a section boundary.

[Trace: linearSolveWrapper::313] ==> Solve at t = 825.241 alpha = 2.99227e+08 tol = 0.33
[Trace: residualDaeWrapper::217] ==> Residual at t = 825.241 sec = 3
[Trace: linearSolveWrapper::313] ==> Solve at t = 825.241 alpha = 1.49613e+08 tol = 0.33
[Trace: residualDaeWrapper::217] ==> Residual at t = 825.241 sec = 3
[Trace: linearSolveWrapper::313] ==> Solve at t = 825.241 alpha = 7.48066e+07 tol = 0.33
[Trace: residualDaeWrapper::217] ==> Residual at t = 825.241 sec = 3
[Trace: linearSolveWrapper::313] ==> Solve at t = 825.241 alpha = 2.99227e+08 tol = 0.33
[Trace: residualDaeWrapper::217] ==> Residual at t = 825.241 sec = 3
[Trace: linearSolveWrapper::313] ==> Solve at t = 825.241 alpha = 1.49613e+08 tol = 0.33
[Trace: residualDaeWrapper::217] ==> Residual at t = 825.241 sec = 3
[Trace: linearSolveWrapper::313] ==> Solve at t = 825.241 alpha = 7.48066e+07 tol = 0.33

giex_stuck.zip

@sleweke Is it already considered that negative nu should be ignored?

Adds a vector valued and a vector valued component dependent parameter class. The ordering of the latter is component-major (i.e., component index changes the slowest). The reaction index is used for the vector element index in the ParameterId.

The characteristic charge nu is extended to a piecewise cubic polynomial in the GIEX adsorption model. The nu of each component can have different polynomial coefficients and breaks, but their amount must be the same over all components. If the breaks are not provided, a single global polynomial is applied. This maintains backwards compatibility.

Renames the validateConfig() function in the parameter handler to validate() and makes it public. Splits configure() into configure() and configureAndValidate(). This allows more fine-grained control, especially in situations where adsorption models require custom configuration code.

Renames the validateConfig() function in the parameter handler to validate() and makes it public. Splits configure() into configure() and configureAndValidate(). This allows more fine-grained control, especially in situations where reaction models require custom configuration code.

Fixes copy and paste bugs in reading optional parameters of the GIEX adsorption model. Also adds a check for strict monotonicity of the polynomial breaks.

Fixes some more bugs in the GIEX adsorption model.

Number of spline pieces was fundamentally wrong.

sleweke requested a review from schmoelder July 21, 2021 00:10

schmoelder reviewed Jul 22, 2021

View reviewed changes

schmoelder force-pushed the feature/giex_spline branch 2 times, most recently from 1977401 to b1601b1 Compare August 19, 2021 07:48

schmoelder force-pushed the feature/giex_spline branch from b1601b1 to 72dea22 Compare January 12, 2022 08:49

sleweke added 8 commits July 21, 2023 10:54

Add vector valued parameters

eb9fde1

Adds a vector valued and a vector valued component dependent parameter class. The ordering of the latter is component-major (i.e., component index changes the slowest). The reaction index is used for the vector element index in the ParameterId.

Fix bugs in GIEX parameter configuration

25c2ea7

Fixes copy and paste bugs in reading optional parameters of the GIEX adsorption model. Also adds a check for strict monotonicity of the polynomial breaks.

Fix more bugs in GIEX adsorption model

c4b2a98

Fixes some more bugs in the GIEX adsorption model.

Fix bug in computing number of GIEX spline pieces

b834e8c

Number of spline pieces was fundamentally wrong.

Remove debug logging

d804b30

schmoelder force-pushed the feature/giex_spline branch from 72dea22 to d804b30 Compare July 21, 2023 08:54

jbreue16 added the binding model label Jun 17, 2024

Uh oh!

Conversation

sleweke commented Jul 21, 2021

Uh oh!

schmoelder commented Jul 21, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

sleweke commented Jul 21, 2021

Uh oh!

schmoelder Jul 22, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

sleweke Jul 22, 2021

Choose a reason for hiding this comment

Uh oh!

schmoelder Jul 22, 2021

Choose a reason for hiding this comment

Uh oh!

sleweke Jul 22, 2021

Choose a reason for hiding this comment

Uh oh!

schmoelder Jul 22, 2021

Choose a reason for hiding this comment

Uh oh!

schmoelder commented Jul 22, 2021

Uh oh!

schmoelder commented Jul 22, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

sleweke commented Jul 22, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

schmoelder commented Jul 22, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

lieres commented Jul 26, 2021

Uh oh!

schmoelder commented Jul 29, 2021

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

schmoelder commented Jul 21, 2021 •

edited

Loading

schmoelder Jul 22, 2021 •

edited

Loading

schmoelder commented Jul 22, 2021 •

edited

Loading

sleweke commented Jul 22, 2021 •

edited

Loading

schmoelder commented Jul 22, 2021 •

edited

Loading