Review, take 2

Follow-up to #1 (closed)

The code looks very nice. My complaints would be mainly personal preference. Some things are bit more important, though:

Example

The incidence matrices seem to have the wrong dimensions. Given that there are 6 vertices and 4 edges, the incidence matrices corresponding to Fig. 2 should be the following:

A_B = \begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 \\
0 & 0 & 1 & 0
\end{bmatrix}\\
A_C = \begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0
\end{bmatrix}\\
A_I = \begin{bmatrix}
0 & 1 & 1 & 1
\end{bmatrix}

Incidence matrices

I could not find where the code handles above incidence matrices. Is the network of Fig. 1 hard-coded at the moment?

Demos two* and four*

The corresponding .cpp files are too similar for my liking. If we have a look at their diff, it seems that all of it could be done by slightly adjusting some config.json files.

diff two*/*.cpp four*/*.cpp

Click to expand

74c74,76
<   const double compr_ratio       = config["compressor"]["compression_ratio"].get<double>();
---
>   const std::string compr_type   = config["compressor"]["type"].get<std::string>();
>   const std::string compr_model  = config["compressor"]["model"].get<std::string>();
>   const double compr_spec        = config["compressor"]["specification"].get<double>();
83c85
<         PHModel::Compressor("FC", "AV", compr_ratio, kappa)
---
>         PHModel::Compressor(compr_type, compr_model, compr_spec, kappa)
101a104
>   // initial guess
104c107,114
<   pipeL_init_momentum.setConstant(mom0/compressors[0].momentum_scale);
---
>   pipeL_init_momentum.setConstant(mom0/network.compressors[0].momentum_scale);
> 
>   if (network.compressors[0].type == "FP") {
>     network.compressors[0].update_compression_ratio(compr_spec/p0);
>   }
>   else if (network.compressors[0].type == "FC") {
>     network.compressors[0].update_compression_ratio(compr_spec);
>   }
107c117
<   pipeR_init_density.setConstant(p0*compr_ratio/(R*outlet_temperature));
---
>   pipeR_init_density.setConstant(p0*network.compressors[0].compression_ratio/(R*outlet_temperature));
114a125
>   // input at boundary
117,118c128,129
<     1.0/compressors[0].momentum_scale,
<     compressors[0].compression_ratio,
---
>     0.0,
>     network.compressors[0].specification,
121a133,140
>   if (network.compressors[0].model == "AV") {
>     u_b(1) = 1.0/std::pow(network.compressors[0].specification, 1/network.compressors[0].isentropic_exponent);
>   }
>   else if (network.compressors[0].model == "AM") {
>     u_b(1) = 1.0;
>   }
> 
>   // configure solver
135c154
<   options.function_tolerance = 1e-8;
---
>   options.function_tolerance = 1e-12;
144d162
< 
173c191
<     guess(0)                 = inlet_pressure/ (R*inlet_temperature);
---
>     guess(0)                 = inlet_pressure/(R*inlet_temperature);
191c209
<     options.function_tolerance = 1e-8;
---
>     options.function_tolerance = 1e-12;

Handling of unknown configurations

... could be improved. For example, Network::set_gas_state would silently don't apply anything if compressors[0].type and compressors[0].model have unknown values. You should at least print a warning but ideally refuse to assemble/configure the network, e.g. by throwing an exception.