max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-16 08:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'ml/glic' into 'poet'

Browse Source 
Code refactoring to allow better support for older compilers & using C=++ features for type-safety and disadvantages of using preprocessor macros

See merge request naaice/iphreeqc!15
This commit is contained in:
Max Lübke 2024-11-07 14:21:49 +01:00
commit 61214a01ad
4 changed files with 26 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
poet/src/Wrapper/SolutionWrapper.hpp
View File

@ -1,6 +1,5 @@
#pragma once #pragma once
#include "NameDouble.h"
#include "Solution.h" #include "Solution.h"
#include "WrapperBase.hpp" #include "WrapperBase.hpp"
#include <array> #include <array>
@ -26,7 +25,7 @@ private:
cxxSolution *solution; cxxSolution *solution;
const std::vector<std::string> solution_order; const std::vector<std::string> solution_order;
static constexpr std::array<std::string, 5> ESSENTIALS = {"H", "O", "Charge", static constexpr std::array<const char *, 5> ESSENTIALS = {"H", "O", "Charge",
"H(0)", "O(0)"}; "H(0)", "O(0)"};
static constexpr std::size_t NUM_ESSENTIALS = ESSENTIALS.size(); static constexpr std::size_t NUM_ESSENTIALS = ESSENTIALS.size();

11
src/phreeqcpp/basicsubs.cpp
View File

@ -1,4 +1,5 @@
#include "Phreeqc.h" #include "Phreeqc.h"
#include "global_structures.h"
#include "phqalloc.h" #include "phqalloc.h"
#include "Utils.h" #include "Utils.h"
@ -173,15 +174,15 @@ sa_declercq(double sa_type, double Sa, double d, double m, double m0, double gfw
double mass0 = m0 * gfw; double mass0 = m0 * gfw;
double V0 = mass0 / d; // volume double V0 = mass0 / d; // volume
double St0 = mass0 * Sa; // total surface double St0 = mass0 * Sa; // total surface
double a0 = pow(3.0 * V0/(4.0 * pi), 1.0/3.0); // ((3*V0)/(4 * 3.14159265359))^(1/3) double a0 = pow(3.0 * V0/(4.0 * piConstant), 1.0/3.0); // ((3*V0)/(4 * 3.14159265359))^(1/3)
double Sp0 = (4.0 * pi) * a0 * a0; // surface particle double Sp0 = (4.0 * piConstant) * a0 * a0; // surface particle
double np = St0 / Sp0; // number of particles double np = St0 / Sp0; // number of particles
double RATS = Sa / St0; double RATS = Sa / St0;
double mass = m * gfw; double mass = m * gfw;
double V = mass / d; double V = mass / d;
double a = pow(3.0 * V/(4.0 * pi), 1.0/3.0); //((3*V)/(4 * 3.14159265359))^(1/3) double a = pow(3.0 * V/(4.0 * piConstant), 1.0/3.0); //((3*V)/(4 * 3.14159265359))^(1/3)
double St = 4.0 * pi * a * a * np; double St = 4.0 * piConstant * a * a * np;
return St * RATS; // total current surface return St * RATS; // total current surface
} }
error_string = sformatf( "Unknown surface area type in SA_DECLERCQ %d.", (int) sa_type); error_string = sformatf( "Unknown surface area type in SA_DECLERCQ %d.", (int) sa_type);
@ -421,7 +422,7 @@ calc_SC(void)
sqrt_q = sqrt(q); sqrt_q = sqrt(q);
// B1 = relaxtion, B2 = electrophoresis in ll = (ll0 - B2 * sqrt(mu) / f2(1 + ka)) * (1 - B1 * sqrt(mu) / f1(1 + ka)) // B1 = relaxtion, B2 = electrophoresis in ll = (ll0 - B2 * sqrt(mu) / f2(1 + ka)) * (1 - B1 * sqrt(mu) / f1(1 + ka))
a = 1.60218e-19 * 1.60218e-19 / (6 * pi); a = 1.60218e-19 * 1.60218e-19 / (6 * piConstant);
B1 = a / (2 * 8.8542e-12 * eps_r * 1.38066e-23 * tk_x) * q / (1 + sqrt_q) * DH_B * 1e10 * z_plus * z_min; // DH_B is per Angstrom (*1e10) B1 = a / (2 * 8.8542e-12 * eps_r * 1.38066e-23 * tk_x) * q / (1 + sqrt_q) * DH_B * 1e10 * z_plus * z_min; // DH_B is per Angstrom (*1e10)
B2 = a * AVOGADRO / viscos_0 * DH_B * 1e17; // DH_B per Angstrom (*1e10), viscos in mPa.s (*1e3), B2 in cm2 (*1e4) B2 = a * AVOGADRO / viscos_0 * DH_B * 1e17; // DH_B per Angstrom (*1e10), viscos in mPa.s (*1e3), B2 in cm2 (*1e4)
//B1 = a / (2 * 8.8542e-12 * eps_c * 1.38066e-23 * tk_x) * q / (1 + sqrt_q) * DH_B * 1e10 * z_plus * z_min; // DH_B is per Angstrom (*1e10) //B1 = a / (2 * 8.8542e-12 * eps_c * 1.38066e-23 * tk_x) * q / (1 + sqrt_q) * DH_B * 1e10 * z_plus * z_min; // DH_B is per Angstrom (*1e10)

29
src/phreeqcpp/global_structures.h
View File

@ -17,19 +17,22 @@
#define MISSING -9999.999 #define MISSING -9999.999
#include "NA.h" /* NA = not available */ #include "NA.h" /* NA = not available */
#define F_C_MOL 96493.5 /* C/mol or joule/volt-eq */ constexpr double F_C_MOL = 96493.5; /* C/mol or joule/volt-eq */
#define F_KJ_V_EQ 96.4935 /* kJ/volt-eq */ constexpr double F_KJ_V_EQ = 96.4935; /* kJ/volt-eq */
#define F_KCAL_V_EQ 23.0623 /* kcal/volt-eq */ constexpr double F_KCAL_V_EQ = 23.0623; /* kcal/volt-eq */
#define R_LITER_ATM 0.0820597 /* L-atm/deg-mol */ constexpr double R_LITER_ATM = 0.0820597; /* L-atm/deg-mol */
#define R_KCAL_DEG_MOL 0.00198726 /* kcal/deg-mol */ constexpr double R_KCAL_DEG_MOL = 0.00198726; /* kcal/deg-mol */
#define R_KJ_DEG_MOL 0.00831470 /* kJ/deg-mol */ constexpr double R_KJ_DEG_MOL = 0.00831470; /* kJ/deg-mol */
#define EPSILON 78.5 /* dialectric constant, dimensionless. Is calculated as eps_r(P, T) in calc_dielectrics. Update the code?? */ constexpr double EPSILON =
#define EPSILON_ZERO 8.854e-12 /* permittivity of free space, C/V-m = C**2/m-J */ 78.5; /* dialectric constant, dimensionless. Is calculated as eps_r(P, T) in
#define JOULES_PER_CALORIE 4.1840 calc_dielectrics. Update the code?? */
#define PASCAL_PER_ATM 1.01325E5 /* conversion from atm to Pa */ constexpr double EPSILON_ZERO =
#define AVOGADRO 6.02252e23 /* atoms / mole */ 8.854e-12; /* permittivity of free space, C/V-m = C**2/m-J */
#define pi 3.14159265358979 constexpr double JOULES_PER_CALORIE = 4.1840;
#define AH2O_FACTOR 0.017 constexpr double PASCAL_PER_ATM = 1.01325E5; /* conversion from atm to Pa */
constexpr double AVOGADRO = 6.02252e23; /* atoms / mole */
constexpr double AH2O_FACTOR = 0.017;
constexpr double piConstant = 3.14159265358979;
#define TRUE 1 #define TRUE 1
#define FALSE 0 #define FALSE 0

3
src/phreeqcpp/utilities.cpp
View File

@ -1,5 +1,6 @@
#include "Utils.h" #include "Utils.h"
#include "Phreeqc.h" #include "Phreeqc.h"
#include "global_structures.h"
#include "phqalloc.h" #include "phqalloc.h"
#include "NameDouble.h" #include "NameDouble.h"
#include "Exchange.h" #include "Exchange.h"
@ -111,7 +112,7 @@ calc_dielectrics(LDBLE tc, LDBLE pa)
= 1.671008e-3 (esu^2 / (erg/K)) / (eps_r * T) */ = 1.671008e-3 (esu^2 / (erg/K)) / (eps_r * T) */
LDBLE e2_DkT = 1.671008e-3 / (eps_r * T); LDBLE e2_DkT = 1.671008e-3 / (eps_r * T);
DH_B = sqrt(8 * pi * AVOGADRO * e2_DkT * rho_0 / 1e3); // Debye length parameter, 1/cm(mol/kg)^-0.5 DH_B = sqrt(8 * piConstant * AVOGADRO * e2_DkT * rho_0 / 1e3); // Debye length parameter, 1/cm(mol/kg)^-0.5
DH_A = DH_B * e2_DkT / (2. * LOG_10); //(mol/kg)^-0.5 DH_A = DH_B * e2_DkT / (2. * LOG_10); //(mol/kg)^-0.5