commit bc81ddd74ed837ea35424ef27d80985fc737e905 Author: Thomas Mabit Date: Tue Nov 19 12:50:18 2024 +0100 Initial commit diff --git a/src/Main.cpp b/src/Main.cpp new file mode 100644 index 0000000..8c1d09f --- /dev/null +++ b/src/Main.cpp @@ -0,0 +1,56 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +struct JNICaller {}; +#include "Simulateur/Simulateur.h" + + +int main() { + // -------- Parameters -------- // + Configuration configuration; + // System dimensions [m] + configuration.L[0] = 2.64e-3; + configuration.L[1] = 4.2e-3; + // Current density [A/m2] + configuration.J = 500; + // Simulation duration [s] + configuration.tf = 4.0;//24.25; + // --------- Topology --------- // + Topology topology; + // Polygon count + topology.polygonCount[0] = 1; + topology.polygonCount[1] = 1; + // Polygon indexs + topology.polygonIndex[0] = new int[] {0, 6}; + topology.polygonIndex[1] = new int[] {6, 12}; + // Polygon defs + Point points[] = { + {0.0, 0.0}, {0.27, 0.0}, {0.44, 0.63}, {0.43, 0.8}, {0.27, 1.0}, {0.0, 1.0}, + {1.0, 0.0}, {1.0, 1.0}, {0.60, 1.0}, {0.72, 0.8}, {0.71, 0.46}, {0.60, 0.0}, + }; topology.polygonDefinitions = points; + + // Simulation + if (!fs::exists(outputFolder)) {fs::create_directory(outputFolder);} outputDelete(outputFolder.c_str()); + Performances performances = simulate(configuration, topology, nullptr); + cout << endl << "Simulation complete (" << (performances.collected[H2]+performances.lost[H2]+performances.remaining[H2]) << "mol of H2 produced) (n=" << (performances.collected[H2] / (performances.collected[H2]+performances.lost[H2])) << "%, R=" << performances.resistor << "ohm.m2, p=" << (performances.collected[H2] / (performances.collected[H2]+performances.contamination[H2])) << ") !" << endl; + + // Cleanup memory + delete[] topology.polygonIndex[0]; + delete[] topology.polygonIndex[1]; + + // Exit + return 0; +} diff --git a/src/Simulateur/Simulateur.h b/src/Simulateur/Simulateur.h new file mode 100644 index 0000000..a5f9790 --- /dev/null +++ b/src/Simulateur/Simulateur.h @@ -0,0 +1,157 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace std; + +#include +#include +#include +#include +typedef Eigen::MatrixXd Mat; +typedef Eigen::SparseMatrix SpMat; +typedef Eigen::Triplet Trip; +using Eigen::VectorXd; +using Eigen::MatrixXd; +using Eigen::MatrixXi; + +#include "params/params.h" +#include "utils/typesdef.h" +#include "utils/math2.h" +#include "params/consts.h" +#include "params/loader.h" +#include "utils/polygon.h" + +#include "utils/progress.h" + +#include "solver/Tension.h" +#include "solver/Current.h" +#include "solver/Topology.h" +#include "solver/Nucleation.h" +#include "solver/Timestep.h" +#include "solver/Dynamics.h" +#include "solver/Counter.h" + +#include "utils/output.h" +#include "utils/loader.h" + +void capture(Configuration& configuration, Topology& topology, const char *path) { + // Load parameters and output mask + Parameters parameters = loadParameters(configuration); + MatrixXi MASK = generateMask(parameters, topology); + outputMatrixI(path, MASK, 0, 2); +} + +Performances simulate(Configuration& configuration, Topology& topology, JNICaller* jniCaller) { + // Load parameters and compute pre-valus + Parameters parameters = loadParameters(configuration); + + // Prepare progressbar + progressInit(jniCaller); + + // Electrical simulation + MatrixXi MASK = generateMask(parameters, topology); + MatrixTension TENSION = solveTension(parameters, MASK, false); + MatrixXd CURRENT = computeCurrent(parameters, MASK, TENSION); + MatrixXd SOURCES = computeSourceNormalized(parameters, MASK, TENSION); + vector Sites = buildNucleationSite(parameters, topology, SOURCES); + + // Bubble dynamics initialization + TimeSteps timesteps = {0, 0, 0, 0}; + ListBubbles Bubbles; Bubbles.flagAttached = false; Bubbles.flagDetached = false; + Bubbles.countAttached = 0; Bubbles.countDetached = 0; + Bubbles.tNextNucleation = 0; + vector Walls = generateWalls(parameters, topology); + ElectrolysisStatistics Statistics = {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, 0.0, 0.0, {0.0, 0.0}, {0.0, 0.0}}; + updateStatistics(Statistics, TENSION, CURRENT); + + // Output electrical (only if not called from jni) + if (outputElectrical && (jniCaller == nullptr)) { + outputMatrixI(outputFolder, "/mask.bmp", MASK, -1, 1); + outputMatrix(outputFolder, "/tension.bmp", TENSION.VALUES, TENSION.min, TENSION.max, -1.0); + outputData(outputFolder, "/tension.csv", parameters, TENSION.VALUES, TENSION.min, TENSION.max); + outputMatrix(outputFolder, "/current.bmp", CURRENT, -CURRENT.maxCoeff(), CURRENT.maxCoeff(), -1.0); + outputPoints(outputFolder, "/source.bmp", SOURCES, MAX(-SOURCES.minCoeff(), SOURCES.maxCoeff()), 7); + outputSites(outputFolder, "/sites.bmp", parameters, Sites); + } + // Bubble dynamics simulation + double t = 0.0; bool stopFlag = false; + while ((t < parameters.tf) && (!stopFlag)) { + // 1| process nucleation + nucleateBubble(Bubbles, Sites, t); + + // 2| process bubble evolution + computeBubble(parameters, Bubbles, Walls, Statistics, t, timesteps.dt, &stopFlag); + + // 3| update tension + if (electricalCoupling) if (t > TENSION.tNext) { + MatrixXi MASK_new = updateMask(parameters, MASK, Bubbles); + if (MASK != MASK_new) { + MASK = MASK_new; + MatrixTension TENSION_new = solveTension(parameters, MASK, true); + if (TENSION_new.sucess) { + TENSION = TENSION_new; + CURRENT = computeCurrent(parameters, MASK, TENSION); + SOURCES = computeSourceNormalized(parameters, MASK, TENSION); + updateNucleationSite(parameters, SOURCES, Sites); + updateStatistics(Statistics, TENSION, CURRENT); + } + TENSION.tNext += timesteps.dtElec; + } + } + + // 4| update timestep + t += timesteps.dt; timesteps.i++; + timesteps = updateTimestep(timesteps, Bubbles.maxRadius, Bubbles.maxVelocity, Sites, &stopFlag); + + // 5| output + progressRefresh1((t/parameters.tf), jniCaller); + if (jniCaller == nullptr) { + // Frame output + outputFrame(outputFolder, "/frame.blb", t, parameters, timesteps, topology, Bubbles, Statistics); + // Electrical outputs + if (outputElectrical && electricalCoupling) { + timesteps.iElec++; + outputMatrixI(outputFolder, "/mask/", timesteps.iElec, MASK, -1, 1); + outputMatrix(outputFolder, "/tension/", timesteps.iElec, TENSION.VALUES, TENSION.min, TENSION.max, -1.0); + outputMatrix(outputFolder, "/current/", timesteps.iElec, CURRENT, -CURRENT.maxCoeff(), CURRENT.maxCoeff(), -1.0); + outputPoints(outputFolder, "/source/", timesteps.iElec, SOURCES, MAX(-SOURCES.minCoeff(), SOURCES.maxCoeff()), 7); + outputSites(outputFolder, "/sites/", timesteps.iElec, parameters, Sites); + } + } + } + + // Performances + Performances performances; + if (!stopFlag) { + // Compute performances + performances = computePerformances(parameters, Statistics, Bubbles); + } else { + // Unable to compute performances + performances = zeroPerformances(); + } + + return performances; +} + diff --git a/src/Simulateur/params/consts.h b/src/Simulateur/params/consts.h new file mode 100644 index 0000000..b08fccd --- /dev/null +++ b/src/Simulateur/params/consts.h @@ -0,0 +1,92 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +/* Global const */ +double Fcst = 96485.332; // Faraday cst [A.s/mol] +double Rcst = 8.314463; // Constante gaz parfait [J/mol.K] + +/// Static pre-computation +double eThick, RSite[2], Fcapilaire[2], FcapilaireSq[2], RFritzExtra[2], Rmini[2], +nRPrefactor, nFritz[2], RFritzMax, RFritzMin, buoyancyPrefactor, QnPrefactor[2], Tprefactor[2], growingKineticFractionWait, uAscentMin, +tafelJ0A0[2], bPrefactor, +interactionLongMaximumRangeSq, interactionLubMaximumRangeSq, interactionRepulsiveRangeSq, interactionCollisionRangeSq, +interactionSafetyHalfRange, interactionSafetyHalfRangeInv, interactionSafetyHalfRangeSq, interactionSafetyRange, interactionSafetyRangeSq; +int bubbleLifeCount; + +void staticPreComputation() { + // Diametre moyen + eThick = (RFritz[0] + RFritz[1]); + for (int i = 0; i < 2; i++) { + // Taille du défaut 'r' associé au site de nucléation pour le rayon de Fritz considérée (tel que Fcap = Fbuoyancy) + RSite[i] = 2*CUBE(RFritz[i])*rho*gravity/(3*gammaLiquidGas[i]); + // Force capilaire sur un site (Fcap = 2 PI r gamma) + Fcapilaire[i] = 2.0*M_PI*RSite[i]*gammaLiquidGas[i]; + FcapilaireSq[i] = SQ(Fcapilaire[i]); + //... + RFritzExtra[i] = (extraRFritz+1)*RFritz[i]; + Rmini[i] = detachementMinimalRadiusNorm*RFritz[i]; + } + // Préfacteur pour le calcul de la quantité de matière dans une sphère (n = [4/3]PI * P0/RT * r^3) + nRPrefactor = 4.0*M_PI*P0/(3.0*Rcst*T0); + // Quantité de matière au détachement de Fritz + for (int i = 0; i < 2; i++) { + nFritz[i] = nRPrefactor*CUBE(RFritz[i]); + } + // ... + RFritzMax = mathMax(RFritz[0], RFritz[1]); + RFritzMin = mathMin(RFritz[0], RFritz[1]); + // Préfacteur pour la force de flottaison Fbuoyancy = [4/3]PI*rho*g * r^3 + buoyancyPrefactor = (4.0/3.0)*M_PI*rho*gravity; + for (int i = 0; i < 2; i++) { + // Préfacteur pour le calcul de la charge (Jelec = [Nelec*F] * Jmolaire) + QnPrefactor[i] = Nelectronic[i]*Fcst; + // Préfacteur pour le calcul de la periode de nucleation + Tprefactor[i] = Nelectronic[i]*Fcst*nFritz[i]; + } + // ... + growingKineticFractionWait = 1.0-growingKineticFraction; + // Plus petite vitesse d'ascension naturelle (pour les plus petites bulles c-a-d les bulles à RFritz) + uAscentMin = 2*rho*gravity*mathSQ(RFritzMin)/(9*mu); + + // Préfacteur de la dérivée de Tafel + for (int i = 0; i < 2; i++) { + tafelJ0A0[i] = tafelJ0[i]/tafelA0[i]; + } + // Préfacteur du terme de flotaison + bPrefactor = (2.0*rho*gravity)/(9.0*mu); // = (4/3)*PI*rho*g*R3 / 6*pi*mu*R + + // Répétition des bulles -> 2 vies + bubbleLifeCount = (repeatYBubble ? 2 : 1); + + // Version caré et inverse des parametres numériques + interactionLongMaximumRangeSq = mathSQ(interactionLongMaximumRange); + interactionLubMaximumRangeSq = mathSQ(interactionLubMaximumRange); + interactionRepulsiveRangeSq = mathSQ(interactionRepulsiveRange); + interactionCollisionRangeSq = mathSQ(interactionCollisionRange); + + interactionSafetyHalfRange = interactionSafetyD+1.0; + interactionSafetyHalfRangeInv = 1.0/interactionSafetyHalfRange; + interactionSafetyHalfRangeSq = mathSQ(interactionSafetyHalfRange); + interactionSafetyRange = 2.0+interactionSafetyD; + interactionSafetyRangeSq = mathSQ(interactionSafetyRange); +} + +#define TafelAnodeJ(deltaphi) (tafelJ0[ANODE]*exp((deltaphi)/tafelA0[ANODE])) +#define TafelCathodeJ(deltaphi) (-tafelJ0[CATHODE]*exp((deltaphi)/tafelA0[CATHODE])) + +#define dTafelAnodeJ(deltaphi) (-tafelJ0A0[ANODE]*exp((deltaphi)/tafelA0[ANODE])) +#define dTafelCathodeJ(deltaphi) (-tafelJ0A0[CATHODE]*exp((deltaphi)/tafelA0[CATHODE])) + diff --git a/src/Simulateur/params/loader.h b/src/Simulateur/params/loader.h new file mode 100644 index 0000000..6eac1d8 --- /dev/null +++ b/src/Simulateur/params/loader.h @@ -0,0 +1,127 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +double str2bool(string value) { + return value == "true"; +} +int str2int(string value) { + size_t pos; + int intValue = stoi(value, &pos); + if (pos != value.length()) throw invalid_argument("null"); + return intValue; +} +double str2double(string value) { + size_t pos; + double doubleValue = stod(value, &pos); + if (pos != value.length()) throw invalid_argument("null"); + return doubleValue; +} + +void loadParamsFromFile(Configuration& configuration) { + ifstream file("params.txt"); + if (file.is_open()) { + string line; + while (getline(file, line)) { + if (line.empty() || line[0] == '#') continue; + // Find '=' + size_t pos = line.find("="); + if (pos == string::npos) { + cerr << "Error while parsing params.txt (" << line << ")" << endl << flush; + exit(1); + } + // Get variable and value + string variable = line.substr(0, pos), value = line.substr(pos + 1); + + // Update + try { + if (variable == "dSite") dSite = str2double(value); + else if (variable == "RFritz") {RFritz[0] = str2double(value); RFritz[1] = str2double(value);} + else if (variable == "RFritz[0]") RFritz[0] = str2double(value); + else if (variable == "RFritz[1]") RFritz[1] = str2double(value); + else if (variable == "theta") {theta[0] = str2double(value); theta[1] = str2double(value);} + else if (variable == "theta[0]") theta[0] = str2double(value); + else if (variable == "theta[1]") theta[1] = str2double(value); + else if (variable == "sigma") sigma = str2double(value); + else if (variable == "electricalCoupling") electricalCoupling = str2bool(value); + else if (variable == "meshN") meshN = str2int(value); + else if (variable == "outputElectrical") outputElectrical = str2bool(value); + else if (variable == "outputFolder") outputFolder = value; + else if (variable == "useSubProgressBar") useSubProgressBar = str2bool(value); + else if (variable == "J") configuration.J = str2double(value); + else if (variable == "tf") configuration.tf = str2double(value); + else throw invalid_argument("null"); + } catch (const invalid_argument& e) { + cerr << "Error while parsing params.txt (" << line << ")" << endl << flush; + exit(1); + } + } + } +} + +// Dynamic pre-computation +struct Parameters { + double L[2]; /* = dynamic precomputation */ // [m] + int N[2]; /* = dynamic precomputation */ + double dL[2]; /* = dynamic precomputation */ // [m] + double tf, tClogg; /* = dynamic precomputation */ // [s] + double J; /* = dynamic precomputation */ // [A/m2] + double I; /* = dynamic precomputation */ // [A] +}; +Parameters loadParameters(Configuration& configuration) { + // Load params file and static precomputation + loadParamsFromFile(configuration); + staticPreComputation(); + + Parameters parameters; + // Dimensions + parameters.L[0] = configuration.L[0]; + parameters.L[1] = configuration.L[1]; + + // Mesh + if (parameters.L[0] > parameters.L[1]) { + parameters.N[0] = meshN; + parameters.N[1] = round(((double) meshN)*(parameters.L[1]/parameters.L[0])); + } else { + parameters.N[1] = meshN; + parameters.N[0] = round(((double) meshN)*(parameters.L[0]/parameters.L[1])); + } + + // Refine mesh + if (electricalCoupling) { + int refineFactor = 0; + for (int c = 0; c < 2; c++) { + int factor = (int) ceil(BUBBLE_MESH_COEFF * (parameters.L[c] / (2*Rmini[c])) / ((double) parameters.N[c])); + if (factor > 1) { + parameters.N[c] *= factor; + refineFactor = mathMax(refineFactor, factor); + } + } + } + + // Compute dL/dR + parameters.dL[0] = configuration.L[0]/(parameters.N[0]-1); + parameters.dL[1] = configuration.L[1]/(parameters.N[1]-1); + + // Current which correspond to the effective electrode surface (S = e x Ly) + parameters.J = configuration.J; + parameters.I = configuration.J*eThick*parameters.L[1]; + + // Others + parameters.tf = configuration.tf; + parameters.tClogg = cloggedTimeFactor * (parameters.L[1] / uAscentMin); + + return parameters; +} diff --git a/src/Simulateur/params/params.h b/src/Simulateur/params/params.h new file mode 100644 index 0000000..42a32c7 --- /dev/null +++ b/src/Simulateur/params/params.h @@ -0,0 +1,107 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +#define ANODE 0 +#define O2 0 +#define CATHODE 1 +#define H2 1 + + +/* --------- Material --------- */ +double dSite = 1.5/(1e-3); // Hole density [1/m] +double RFritz[2] = {60e-6, 60e-6}; // Rayon de Fritz [m] +double theta[2] = {M_PI, M_PI}; // Angle de mouillage des bulles (de l'ordre de PI) [rad] +/* -------- Electrolyte ------- */ +double mu = 1e-3; // Viscosité dynamique [Pa.s] +double gravity = 9.81; // Gravity [m/s2] +double rho = 997; // Density [kg/m3] +double P0 = 1.01e5; // Pressure [Pa] +double T0 = 273.15 + 35.0; // Temperature [K] +double gammaLiquidGas[2] = {70e-3, 70e-3}; // Surface tensions [N/m] +/* ------ Electrochemistry ---- */ +int Nelectronic[2] = {4, 2}; // Nombre d'électron nécessaire à la formation d'une molécule +double sigma = 2; // Electrolyte conductivity [S/m] +double tafelA0[2] = {(60e-3)/log(10), (30e-3)/log(10)}; // Tafel modeling : voltage slope [V] +double tafelJ0[2] = {1e-3, 1e-1}; // Tafel modeling : exchange current [A/m2] +/* ---------------------------- */ + + + + +/* ------------ Coupling ---------- */ +bool electricalCoupling = false; // Couplage faible du champ électrique +/* ------------------------------- */ + +/* ----------- Repeat y ---------- */ +bool repeatYBubble = false; // Simule la répétition du motif selon y (recopie des bulles 1 fois) +/* ------------------------------- */ + +/* -------------- Mesh ------------ */ +int meshN = 400; // Mesh size [-] +double BUBBLE_MESH_COEFF = 5.0; // Bubble meshing factor for bubble diameter [-] +/* -------------------------------- */ + +/* ------ Solver tolerance ----- */ +#define TOLERANCE_LU 1e-10 // Tolérance sur la vérification de l'inverse obtenu via la méthode LU +#define TOLERANCE_POTENTIAL 1e-6 // Critère d'arrêt sur la résolution du potentiel électrolytique +#define TOLERANCE_CURRENT 1e-3 // Critère d'arrêt sur la résolution du courant +#define DERIVATIVE_VARIATION 1e-2 // Fraction de la grandeur pour l'approximation de la dérivée +#define SOLVER_MAX_ITERATION 40 // Nombre maximum d'itération (pour les cas non solvables) +/* ----------------------------- */ + +/*---------- Time step -------- */ +double GROWING_TIME_COEF = 5.0; // Growing time step factor [-] +double BUOYANCY_TIME_COEF = 10.0; // Buoyancy time step factor [-] +double ELECTRICAL_TIME_COEF = 1.0; // Buoyancy time step factor [-] +/* ---------------------------- */ + + + +/* -- Numerical parameters -- */ +/* -------------------------- */ +double growingKineticFraction = 0.7; // Proportion de croissance/attente sur un site (croissance) +double detachementMinimalRadiusNorm = 0.8; // Fraction du rayon de détachement naturel à partir duquel le détachement dynamique est possible +/* -------------------------- */ +#define interactionLongMaximumRange 15 // Portée maximale de l'interaction champ lointain (normalisée par le rayon moyen) +#define interactionLubMaximumRange 3.5 // Portée maximale de l'interaction lubrification (normalisée par le rayon moyen) +/*--------- Repulsive ------- */ +double interactionRepulsiveRange = 2.2; // Portée du modèle de répulsion (Rref) +double interactionRepulsiveVelocity = 5e0; // Raideur du modèle de répulsion (Uref) +/* ----- Anti-collision ----- */ +#define interactionCollisionRange 3 // Portée à partir de laquelle le contact est testé (normalisée par le rayon moyen) +#define interactionSafetyD 0.1 // Distance de sécurité à partir de laquelle la vitesse normale relative est absorbée +#define extraRFritz 0.02 // Rayon supplémentaire au rayon de Fritz pour éviter le blocage numérique +/* -------------------------- */ +#define cloggedTimeFactor 1 // Temps à partir duquel une bulle est considérée coincée +/* -------------------------- */ + + + + +/* ---------- Ouptut --------- */ +/* --------------------------- */ +bool outputElectrical = true; +string outputFolder = "output"; +int outputMode = 0; +// Mode 0 : all frames +// Mode 1 : real time recording during the entire simulation +int outputFps = 60, outputTimeFactor = 1/10; +// Mode 2 : all frame between given interval +double outputStart = 0, outputStop = 1; +/* --------------------------- */ +bool useSubProgressBar = true; +/* --------------------------- */ + diff --git a/src/Simulateur/solver/Counter.h b/src/Simulateur/solver/Counter.h new file mode 100644 index 0000000..79e6b28 --- /dev/null +++ b/src/Simulateur/solver/Counter.h @@ -0,0 +1,63 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + + + +/* --- Compute the extraction proportions --- */ +/* ------------------------------------------ */ +void updateRemaining(ElectrolysisStatistics& Statistics, ListBubbles& Bubbles) { + Statistics.remaining[0] = 0; Statistics.remaining[1] = 0; + for (list::iterator it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) { + Statistics.remaining[0] += it->quantity[0]; + Statistics.remaining[1] += it->quantity[1]; + } + for (list::iterator it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + Statistics.remaining[0] += it->quantity[0]; + Statistics.remaining[1] += it->quantity[1]; + } +} +/* ------------------------------------------ */ + + +/* ---------- Compute performances ---------- */ +/* ------------------------------------------ */ +Performances zeroPerformances() { + return {{1e9, 1e9}, {0.0, 0.0}, {0.0, 0.0}, {1e9, 1e9}, {0.0, 0.0}, 1e9}; +} +Performances computePerformances(Parameters& parameters, ElectrolysisStatistics& Statistics, ListBubbles& Bubbles) { + updateRemaining(Statistics, Bubbles); + Performances performances = zeroPerformances(); + + for (int s = 0; s < 2; s++) { + performances.electrolysed[s] = (parameters.I*parameters.tf) / QnPrefactor[s]; + performances.collected[s] = Statistics.collected[s]; + performances.contamination[s] = Statistics.contamination[s]; + performances.lost[s] = Statistics.lost[s]; + performances.remaining[s] = Statistics.remaining[s]; + } + + performances.resistor = Statistics.resistor; + + return performances; +} + +void updateStatistics(ElectrolysisStatistics& Statistics, MatrixTension& TENSION, MatrixXd& CURRENT) { + Statistics.resistor = TENSION.resistor; + Statistics.overpotential = (TENSION.phi_anode-TENSION.phi_cathode); + Statistics.tension[0] = TENSION.min; Statistics.tension[1] = TENSION.max; + Statistics.current[0] = CURRENT.minCoeff(); Statistics.current[1] = CURRENT.maxCoeff(); +} +/* ------------------------------------------ */ diff --git a/src/Simulateur/solver/Current.h b/src/Simulateur/solver/Current.h new file mode 100644 index 0000000..f155000 --- /dev/null +++ b/src/Simulateur/solver/Current.h @@ -0,0 +1,106 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +// Calcul la carte des densités de courant sqrt(Jx2 + Jy2) +// J = sigma.grad(phi) +MatrixXd computeCurrent(Parameters& parameters, MatrixXi& MASK, MatrixTension& TENSION) { + int* N = parameters.N; + double dxInvSigma = sigma/parameters.dL[0], dyInvSigma = sigma/parameters.dL[1]; + MatrixXd CURRENT(N[0], N[1]); + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) { + // Electrolyte + if (MASK(i,j) == 0) { + double Jx = 0, Jy = 0; + // dx + if (i > 0) { + // electrolyte + if (MASK(i-1,j) == 0) Jx = dxInvSigma*(TENSION.VALUES(i,j) - TENSION.VALUES(i-1,j)); + // anode + else if (MASK(i-1,j) == 1) Jx = TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j)); + // cathode + else if (MASK(i-1,j) == -1) Jx = TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode); + } + // dy + int jm = (j > 0) ? j-1 : (N[1]-1); + // electrolyte + if (MASK(i,jm) == 0) Jy = dyInvSigma*(TENSION.VALUES(i,j) - TENSION.VALUES(i,jm)); + // anode + else if (MASK(i,jm) == 1) Jy = TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j)); + // cathode + else if (MASK(i,jm) == -1) Jy = TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode); + + + // norm + CURRENT(i,j) = sqrt(SQ(Jx) + SQ(Jy)); + } + // Insulation + else if (MASK(i, j) == -2) { + CURRENT(i, j) = -1; + continue; + } + // Electrode (anode/cathode) + else if ((MASK(i, j) == 1) || (MASK(i, j) == -1)) { + CURRENT(i, j) = 0; + } + } + return CURRENT; +} + +// Calcul la quantité de courant arrivant en chaque noeud frontière d'une électrode +// I = S.Jtafel +MatrixXd computeSourceNormalized(Parameters& parameters, MatrixXi& MASK, MatrixTension& TENSION) { + int* N = parameters.N; + double Sx = eThick*parameters.dL[1], Sy = eThick*parameters.dL[0]; + MatrixXd SOURCES(N[0], N[1]); + + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) { + double src = 0; + if (MASK(i,j) == 0) { + // i-1 + if (i > 0) { + // anode + if (MASK(i-1,j) == 1) src += Sx*TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j)); + // cathode + else if (MASK(i-1,j) == -1) src += Sx*TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode); + } + + // i+1 + if (i < (N[0]-1)) { + // anode + if (MASK(i+1,j) == 1) src += Sx*TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j)); + // cathode + else if (MASK(i+1,j) == -1) src += Sx*TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode); + } + + // j-1 + int jm = (j > 0) ? j-1 : (N[1]-1); + // anode + if (MASK(i,jm) == 1) src += Sy*TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j)); + // cathode + else if (MASK(i,jm) == -1) src += Sy*TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode); + + // j+1 + int jp = (j < (N[1]-1)) ? j+1 : 0; + // anode + if (MASK(i,jp) == 1) src += Sy*TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j)); + // cathode + else if (MASK(i,jp) == -1) src += Sy*TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode); + } + SOURCES(i,j) = src; + } + return SOURCES; +} + diff --git a/src/Simulateur/solver/Dynamics.h b/src/Simulateur/solver/Dynamics.h new file mode 100644 index 0000000..41b9e18 --- /dev/null +++ b/src/Simulateur/solver/Dynamics.h @@ -0,0 +1,422 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + + + +/* -------- Compute bubble evolution --------- */ +/* ------------------------------------------- */ +#include "utils/funinter.h" +#include "utils/fundyn.h" +#include "utils/funvelo.h" +void computeBubble(Parameters& parameters, ListBubbles& Bubbles, vector& Walls, ElectrolysisStatistics& Statistics, double t, double dt, bool* stopFlag) { + // 0| Prepare variables + Bubbles.maxRadius = RFritzMax; Bubbles.maxVelocity = 0; + list::iterator it; + for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) { + it->interactor.count = 0; // Reset bubble-bubble interactions + } + for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + it->interactor.count = 0; // Reset bubble-bubble interactions + } + + + // 1-2| Determine interactions and build mobility matrix + std::vector M_coefs; // Matrix coefficients for velocity system to solve + double Mij[4] { 0, 0, 0, 0 }; // Temp matrix coefficients + + int i = -1; + for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + i++; + // 1| Test which wall and bubbles are supposed to interact significantly + checkWalls(it->interactor.walls, Walls, it->x, it->y, 0.5*interactionLongMaximumRange*it->r); + checkBubbles(it, i, Bubbles, parameters.L[1], stopFlag); + + // 2| Computing mobility matrix - line i + setMobilityMatrixSphere_I(Mij); + addMobilityMatrixPlane(Mij, &(*it)); + M_coefs.push_back(Trip(2*i+0, 2*i+0, Mij[0])); M_coefs.push_back(Trip(2*i+0, 2*i+1, Mij[1])); + M_coefs.push_back(Trip(2*i+1, 2*i+0, Mij[2])); M_coefs.push_back(Trip(2*i+1, 2*i+1, Mij[3])); + + double sqI = SQ(it->r); + for (int jj = 0; jj < it->interactor.count; jj++) { + int j = it->interactor.interactionIndexs[jj]; + if (j == -1) continue; + + Bubble* jt = it->interactor.bubbles[jj]; + if (jt->r == 0) continue; + + double sqJ = SQ(jt->r); + int jMode = it->interactor.imagedModes[jj]; + + if (it->interactor.distancesSqNorm[jj] < interactionLongMaximumRangeSq) { + setMobilityMatrixSphere_J(Mij, 0.75*it->r, sqI, sqJ, it->x-jt->x, it->y-yRepeat(jt->y, parameters.L[1], jMode)); + M_coefs.push_back(Trip(2*i+0, 2*j+0, Mij[0])); M_coefs.push_back(Trip(2*i+0, 2*j+1, Mij[1])); + M_coefs.push_back(Trip(2*i+1, 2*j+0, Mij[2])); M_coefs.push_back(Trip(2*i+1, 2*j+1, Mij[3])); + } + } + } + + // 3| Invert mobility matrix + SpMat Minv; + if (Bubbles.countDetached > 0) { + SpMat M(2*Bubbles.countDetached, 2*Bubbles.countDetached); + M.setFromTriplets(M_coefs.begin(), M_coefs.end()); + + if (invertLU(Minv, M) != 0) { + // Decomposition failed + errorProgress("mobility matrix inverting error @t=" + std::to_string(t) + "s"); + *stopFlag = true; + return; + } + } + + // 4| Build resistance matrix and source term + std::vector R_coefs; // Matrix coefficients for velocity system to solve + double Rij[4] { 0, 0, 0, 0 }; // Temp matrix coefficients + Eigen::VectorXd U0(2*Bubbles.countDetached); // Source term + + i = -1; + for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + i++; + // Build source term + U0[2*i+0] = 0; + U0[2*i+1] = -bPrefactor*SQ(it->r); + + // Add repulsive term + addSphereSphereRepulsive(U0, &(*it), i); + + // Add inverted mobility matrix + Rij[0] = -Minv.coeff(2*i+0, 2*i+0); Rij[1] = -Minv.coeff(2*i+0, 2*i+1); + Rij[2] = -Minv.coeff(2*i+1, 2*i+0); Rij[3] = -Minv.coeff(2*i+1, 2*i+1); + + // Add sphere-sphere lubrification + addSphereSphereContact_I(Rij, parameters.L[1], &(*it)); + // Add sphere-plan lubrification + addSpherePlanContact_I(Rij, &(*it)); + + // Build resistance matrix ii + R_coefs.push_back(Trip(2*i+0, 2*i+0, Rij[0])); R_coefs.push_back(Trip(2*i+0, 2*i+1, Rij[1])); + R_coefs.push_back(Trip(2*i+1, 2*i+0, Rij[2])); R_coefs.push_back(Trip(2*i+1, 2*i+1, Rij[3])); + + for (int jj = 0; jj < it->interactor.count; jj++) { + int j = it->interactor.interactionIndexs[jj]; + if (j == -1) continue; + + Bubble* jt = it->interactor.bubbles[jj]; + if (jt->r == 0) continue; + + int jMode = it->interactor.imagedModes[jj]; + + if (it->interactor.distancesSqNorm[jj] < interactionLongMaximumRangeSq) { + // Add inverted mobility matrix + double ratio_ji = jt->r/it->r; + Rij[0] = -ratio_ji*Minv.coeff(2*i+0, 2*j+0); Rij[1] = -ratio_ji*Minv.coeff(2*i+0, 2*j+1); + Rij[2] = -ratio_ji*Minv.coeff(2*i+1, 2*j+0); Rij[3] = -ratio_ji*Minv.coeff(2*i+1, 2*j+1); + + // Add sphere-sphere lubrification + if (it->interactor.distancesSqNorm[jj] < interactionLubMaximumRangeSq) { + addSphereSphereContact_J(Rij, parameters.L[1], &(*it), jt, j, jMode); + } + + // Build resistance matrix ij + R_coefs.push_back(Trip(2*i+0, 2*j+0, Rij[0])); R_coefs.push_back(Trip(2*i+0, 2*j+1, Rij[1])); + R_coefs.push_back(Trip(2*i+1, 2*j+0, Rij[2])); R_coefs.push_back(Trip(2*i+1, 2*j+1, Rij[3])); + } + } + } + + // 5| Solve velocity + Eigen::VectorXd U; + if (Bubbles.countDetached > 0) { + SpMat R(2*Bubbles.countDetached, 2*Bubbles.countDetached); + R.setFromTriplets(R_coefs.begin(), R_coefs.end()); + Eigen::SparseLU solver; + solver.compute(R); + if (solver.info() != Eigen::Success) { + // Decomposition failed + errorProgress("resistance matrix solving error @t=" + std::to_string(t) + "s"); + *stopFlag = true; + return; + } + U = solver.solve(U0); + } + + // 6| Update velocity of detached bubbles + i = -1; + for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + i++; + it->ux = U(2*i+0); + it->uy = U(2*i+1); + } + + // 7| Check attached bubble detachment + double Fdyn[2]; + for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) { + // Si le détachement dynamique est possible + if (it->r != 0) { + // Calcul des actions hydrodynamiques + Fdyn[0] = 0; Fdyn[1] = 0; + for (int kk = 0; kk < it->interactor.count; kk++) { + Bubble *kt = it->interactor.bubbles[kk]; + int kMode = it->interactor.imagedModes[kk]; + if (it->interactor.distancesSqNorm[kk] < interactionCollisionRangeSq) { + double deltaX = (it->x - kt->x), deltaY = (it->y - yRepeat(kt->y, parameters.L[1], kMode)); + double norm = kt->ux*deltaX + kt->uy*deltaY; + if (norm > 0) { + // Si contact avec une autre bulle alors elle détache immédiatement + it->statLife = 0; + break; + } + } else if (it->interactor.distancesSqNorm[kk] < interactionLubMaximumRangeSq) { + computeSphereSphereForce(Fdyn, &(*it), kt, kMode, parameters.L[1]); + } + } + // Ajout de la pesanteur et retrait de la composante reprise par la paroi + Fdyn[1] += buoyancyPrefactor*CUBE(it->r); + double proj = Fdyn[0]*it->site->nx + Fdyn[1]*it->site->ny; + if (proj < 0) { + Fdyn[0] -= proj*it->site->nx; + Fdyn[1] -= proj*it->site->ny; + } + // Test de détachement hydrodynamique + if ((SQ(Fdyn[0]) + SQ(Fdyn[1]) > FcapilaireSq[it->site->type])) { + it->statLife = 0; + } + } + } + + // 8-9| Velocity correction and integration + i = -1; + for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + // 8| Velocity correction + i++; + bubbleVelocityCorrection(parameters.L[1], &(*it), stopFlag); + } + for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + // 9| Velocity integration + it->x = it->x + it->ux*dt; + it->y = it->y + it->uy*dt; + double velocity = sqrt(SQ(it->ux) + SQ(it->uy)); + if (velocity > Bubbles.maxVelocity) Bubbles.maxVelocity = velocity; + // Check for clogging + if ((it->ux == 0) && (it->uy == 0)) it->tClogg += dt; + else it->tClogg = 0; + + if (it->tClogg > parameters.tClogg) { + errorProgress("clogging detected @t=" + std::to_string(t) + "s"); + *stopFlag = true; + return; + } + // Update wall interactions for next step (detached bubble growing) + checkWalls(it->interactor.walls, Walls, it->x, it->y, 0.5*interactionLongMaximumRange*it->r); + } + + // ?| Detached bubble growing /* !!!!!!! DISABLED (\/ see end of this file \/) !!!!!!! */ + + // 10| Attached bubble growing + for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) { + // Prise en compte du dr propre au site + Site* site = it->site; + double dr = MIN(site->growth*(dt+site->tmiss), MAX(0, RFritzExtra[site->type] - it->r)); + site->tmiss = 0; + + double maxDr = INF; + double itx = it->x + interactionSafetyHalfRange*dr*site->nx, ity = it->y + interactionSafetyHalfRange*dr*site->ny; + + // Detection des parois ... + checkWalls(it->interactor.walls, Walls, itx, ity, 0.5*interactionLongMaximumRange*it->r); + for (int kk = 0; kk < wallMaxmimumIndicators; kk++) { + WallIndicator *indicator = &it->interactor.walls[kk]; + if (indicator->distance < INF) { + double deltaR = (it->r+dr); + if (SQ(indicator->distance) < interactionSafetyHalfRangeSq*SQ(deltaR)) { + double drRemaining = interactionSafetyHalfRangeInv*indicator->distance - it->r; + if (drRemaining < maxDr) { + maxDr = drRemaining; + } + } + } + } + // ... et des bulles empechant la croissance + for (int kk = 0; kk < it->interactor.count; kk++) { + Bubble* kt = it->interactor.bubbles[kk]; + int kMode = it->interactor.imagedModes[kk]; + double deltaX = (kt->x - itx), deltaY = (yRepeat(kt->y, parameters.L[1], kMode) - ity), deltaR = (it->r+dr + kt->r); + double dSQ = deltaX*deltaX + deltaY*deltaY; + if (dSQ < interactionSafetyHalfRangeSq*SQ(deltaR)) { + double drRemaining = interactionSafetyHalfRangeInv*sqrt(dSQ) - it->r - kt->r; + if (drRemaining < maxDr) { + maxDr = drRemaining; + } + } + } + + // Si le grossissement a été limitée, il faut le reporter + if (maxDr != INF) { + double drBis = MIN(MAX(0, maxDr), dr); + site->tmiss = (dr-drBis)/site->growth; + dr = drBis; + } + + // La croissance ne peut commencer qu'à partir d'un rayon mini + if (it->r == 0) { + if ((it->r+dr) < Rmini[site->type]) { + site->tmiss += dr/site->growth; + dr = 0; + } + } + + // Finalement, si grossissement il y a : + if (dr != 0) { + it->quantity[site->type] += nRPrefactor*(mathCB(it->r + dr)-mathCB(it->r)); + it->r = it->r + dr; + it->x = it->x + interactionSafetyHalfRange*dr*site->nx; + it->y = it->y + interactionSafetyHalfRange*dr*site->ny; + } + } + + // 11| Detached bubble exit + list::iterator> toKillDetached; + for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + if (it->x < 0.0) { // left exit (0) + it->statLife = 0; + // Add to stats + Statistics.collected[0] += it->quantity[0]; + Statistics.contamination[0] += it->quantity[1]; + Statistics.lost[1] += it->quantity[1]; + } else if (it->x > parameters.L[0]) { // right exit (1) + it->statLife = 0; + // Add to stats + Statistics.collected[1] += it->quantity[1]; + Statistics.contamination[1] += it->quantity[0]; + Statistics.lost[0] += it->quantity[0]; + } else if (it->y > parameters.L[1]) { // top recopy + it->y = it->y - parameters.L[1]; + it->statLife = it->statLife - 1; + // Add to stats + if (it->statLife == 0) { + Statistics.lost[0] += it->quantity[0]; + Statistics.lost[1] += it->quantity[1]; + } + } + if (it->statLife == 0) { + Bubbles.countDetached--; + Bubbles.flagDetached = true; + toKillDetached.push_front(it); + continue; + } + } + if (Bubbles.flagDetached) { // If a bubble dies, reorder the list + for (list::iterator>::iterator itKill = toKillDetached.begin(); itKill != toKillDetached.end(); itKill++) Bubbles.detached.erase(*itKill); + Bubbles.flagDetached = false; + } + + // 12| Attached bubble detachment + list::iterator> toKillAttached; + for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) { + if (it->statLife == 0) { + // Changement de statut + Bubbles.flagAttached = true; + it->statLife = bubbleLifeCount; + toKillAttached.push_front(it); + Bubbles.detached.push_back(*it); + Bubbles.countAttached--; Bubbles.countDetached++; + // Prochaine nucléation + it->site->t = t + growingKineticFractionWait*it->site->T; + it->site->occuped = false; + } + } + + // 13| Reorder bubble list + if (Bubbles.flagAttached) { // If a bubble dies, reorder the list + for (list::iterator>::iterator itKill = toKillAttached.begin(); itKill != toKillAttached.end(); itKill++) Bubbles.attached.erase(*itKill); + Bubbles.flagAttached = false; + } +} +/* ------------------------------ */ + + + + + + +// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DISABLED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +// ?| Detached bubble growing +/*double dr = dt*growingFactor; +int growingBubbleCount; Bubble* growingCollisionBubbles[2]; int growingCollisionBubblesMode[2]; double growingCollisionBubblesD[2]; +WallIndicator* growingCollisionWall = 0; int growingWallCount; +for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + growingBubbleCount = 0; growingWallCount = 0; + // Detection des bulles empechant la croissance + for (int kk = 0; kk < it->interactor.count; kk++) { + if (it->interactor.distancesSqNorm[kk] < interactionCollisionRangeSq) { + Bubble* kt = it->interactor.bubbles[kk]; + int kMode = it->interactor.imagedModes[kk]; + double deltaX = (kt->x - it->x), deltaY = (Y_repeat(kt->y, kMode) - it->y), deltaR = (it->r+dr + kt->r); + double dSQ = deltaX*deltaX + deltaY*deltaY; + if (dSQ < interactionSafetyHalfRangeSq*SQ(deltaR)) { + if (growingBubbleCount == 2) break; + growingCollisionBubbles[growingBubbleCount] = kt; + growingCollisionBubblesD[growingBubbleCount] = sqrt(dSQ); + growingCollisionBubblesMode[growingBubbleCount] = kMode; + growingBubbleCount++; + } + } + } + // Detection des parois empechant la croissance + for (int kk = 0; kk < wallMaxmimumIndicators; kk++) { + WallIndicator *indicator = &it->interactor.walls[kk]; + if (indicator->distance < (0.5*it->r*interactionCollisionRange)) { + // Recalcul de la distance exact (à cause de l'intégration) + if (indicator->distance < interactionSafetyHalfRange*(it->r+dr)) { + if (growingWallCount == 1) break; + growingCollisionWall = indicator; + growingWallCount++; + } + } + } + if ((growingBubbleCount != 0) || (growingWallCount != 0)) { + // Simple decallage si une bulle en blocage avec une autre + if ((growingBubbleCount == 1) && (growingWallCount == 0)) { + double nx = it->x - growingCollisionBubbles[0]->x, ny = it->y - Y_repeat(growingCollisionBubbles[0]->y, growingCollisionBubblesMode[0]); + double d_inv = 1.0/growingCollisionBubblesD[0]; + nx = d_inv*nx; ny = d_inv*ny; + + it->x = it->x + dr*nx; + it->y = it->y + dr*ny; + // Idem avec une paroi + } else if ((growingBubbleCount == 0) && (growingWallCount == 1)) { + it->x = it->x + dr*growingCollisionWall->n.x; + it->y = it->y + dr*growingCollisionWall->n.y; + // Résolution de la position si deux bulles + } else if ((growingBubbleCount == 2) && (growingWallCount == 0)) { + continue; + // Ou une bulle et une paroi + } else if ((growingBubbleCount == 1) && (growingWallCount == 1)) { + continue; + // Pas de croissance si plus de deux bulles en blocage ou une bulle et un mur + } else { + continue; + } + } + + // Croissance + it->r = it->r + dr; + double quantity = quantityFromDr(it->r, it->r + dr), tension = getTension(parameters, TENSION, it->x, it->y); + it->quantity[O2] += tension*quantity; it->quantity[H2] += (1-tension)*quantity; + if (it->r > Bubbles.maxRadius) + Bubbles.maxRadius = it->r; +}*/ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DISABLED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/src/Simulateur/solver/Nucleation.h b/src/Simulateur/solver/Nucleation.h new file mode 100644 index 0000000..7b43ea0 --- /dev/null +++ b/src/Simulateur/solver/Nucleation.h @@ -0,0 +1,163 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + + + +/* ------------ Boundaries ------------ */ +/* ------------------------------------ */ +void updateNucleationSite(Parameters& parameters, MatrixXd& Sources, vector& Sites) { + int* N = parameters.N; + // Create a clone of sourcesValue since it will be modified + MatrixXd SourcesCopy(N[0], N[1]); + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) { + SourcesCopy(i,j) = (double) Sources(i,j); + } + + // Setup delta step from hole density and check that it's not too small. + double delta = 1.0/dSite; + int di2 = round((delta / (2.0*parameters.dL[0])) - 0.5), dj2 = round((delta / (2.0*parameters.dL[1])) - 0.5); // -0.5 car conversion longueur case + di2 = MAX(1, di2); dj2 = MAX(1, dj2); + + // Update all sites + for (vector::iterator it = Sites.begin(); it != Sites.end(); ++it) { + Point pos = {it->x, it->y}; + int i, j; positionToIndex(N, &i, &j, pos.x/parameters.L[0], pos.y/parameters.L[1]); + double I = 0.0; int count = 0; + // Moyennage sur un rayon de 1/2d + for (int ii = -di2; ii <= +di2; ii++) for (int jj = -dj2; jj <= +dj2; jj++) { + if ((i+ii > -1) && (i+ii < N[0]) && (j+jj > -1) && (j+jj < N[1])) { + double dii = 0.0, djj = 0.0; + if (di2 > 0) dii = ((double) ii) / di2; + if (dj2 > 0) djj = ((double) jj) / dj2; + if (mathSQ(dii) + mathSQ(djj) <= 1.0) { + if (mathSign(SourcesCopy(i+ii, j+jj)) == (it->type == 0 ? +1 : -1)) { + I += ABS(SourcesCopy(i+ii, j+jj)); + SourcesCopy(i+ii, j+jj) = 0; + count++; + } + } + } + } + if (I > 0.0) { + double T = Tprefactor[it->type] / I; + it->T = T; + it->growth = RFritz[it->type]/(T*growingKineticFraction); + } else it->T = INF; + } +} + +vector buildNucleationSite(Parameters& parameters, Topology& topology, MatrixXd& Sources) { + int* N = parameters.N; + mathRandInit(); + // Create a clone of sourcesValue since it will be modified + MatrixXd SourcesCopy(N[0], N[1]); + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) { + SourcesCopy(i,j) = (double) Sources(i,j); + } + + // Setup delta step from hole density and check that it's not too small. + double delta = 1.0/dSite; + int di2 = round((delta / (2.0*parameters.dL[0])) - 0.5), dj2 = round((delta / (2.0*parameters.dL[1])) - 0.5); // -0.5 car conversion longueur case + if ((di2 == 0) || (dj2 == 0)) { + errorProgress("mesh refinement too small for nucleation."); + exit(1); + } + di2 = MAX(1, di2); dj2 = MAX(1, dj2); + list SiteList; + + // For each poly.. + for (int s = 0; s < 2; s++) for (int p = 0; p < topology.polygonCount[s]; p++) { + int n = topology.polygonIndex[s][p+1] - topology.polygonIndex[s][p]; + double pos = 0.5*delta; + // For each segment... + for (int seg = 0; seg < n; seg++) { + int index = topology.polygonIndex[s][p] + seg; + int next = topology.polygonIndex[s][p] + (seg+1)%n; + Point x0 = {topology.polygonDefinitions[index].x, topology.polygonDefinitions[index].y}; + Point x1 = {topology.polygonDefinitions[next].x, topology.polygonDefinitions[next].y}; + + // Suppression des frontière + if ( ((x0.x == 0.0) && (x1.x == 0.0)) || ((x0.x == 1.0) && (x1.x == 1.0)) + || ((x0.y == 0.0) && (x1.y == 0.0)) || ((x0.y == 1.0) && (x1.y == 1.0)) ) + continue; + + // DeNormalization + x0.x *= parameters.L[0]; x1.x *= parameters.L[0]; + x0.y *= parameters.L[1]; x1.y *= parameters.L[1]; + Vect u = {x1.x - x0.x, x1.y - x0.y}; + double length = vectLength(u); + u.x = u.x / length; u.y = u.y / length; + Vect n = {u.y, -u.x}; + + // Parcours du segmenet + while (pos <= length) { + Point current = vectForward(x0, u, pos); + int i, j; positionToIndex(N, &i, &j, current.x/parameters.L[0], current.y/parameters.L[1]); + double I = 0.0; int count = 0; + // Moyennage sur un rayon de 1/2d + for (int ii = -di2; ii <= +di2; ii++) for (int jj = -dj2; jj <= +dj2; jj++) { + if ((i+ii > -1) && (i+ii < N[0]) && (j+jj > -1) && (j+jj < N[1])) { + double dii = 0.0, djj = 0.0; + if (di2 > 0) dii = ((double) ii) / di2; + if (dj2 > 0) djj = ((double) jj) / dj2; + if (mathSQ(dii) + mathSQ(djj) <= 1.0) { + if (mathSign(SourcesCopy(i+ii, j+jj)) == (s == 0 ? +1 : -1)) { + I += ABS(SourcesCopy(i+ii, j+jj)); + SourcesCopy(i+ii, j+jj) = 0; + count++; + } + } + } + } + if (I > 0.0) { + double T = Tprefactor[s] / I; + double randomT = T*growingKineticFractionWait*mathRand(); + SiteList.push_back({current.x, current.y, T, RFritz[s]/(T*growingKineticFraction), n.x, n.y, s, randomT, 0.0, 0}); + } + pos += delta; + } + pos -= length; + } + } + vector Sites(SiteList.begin(), SiteList.end()); + return Sites; +} +/* ------------------------------------ */ + + + + + +/* ----------- Nucleation ------------ */ +/* ----------------------------------- */ +void nucleateBubble(ListBubbles& Bubbles, vector& Sites, double t) { + if (t < Bubbles.tNextNucleation) return; + Bubbles.tNextNucleation = INF; + for (vector::iterator it = Sites.begin(); it != Sites.end(); ++it) { + if (!it->occuped) if (t >= it->t) { + // Pop a new bubble ! + it->occuped = true; + it->tmiss = 0; + Bubbles.attached.push_front({it->x, it->y, 0.0, 0.0, 0.0, 1, 0, true, &(*it), {}, {0.0, 0.0}, 0.0}); + Bubbles.countAttached++; + } + // May be the next one ? + if (it->t < Bubbles.tNextNucleation) { + Bubbles.tNextNucleation = it->t; + } + } +} +/* ----------------------------------- */ diff --git a/src/Simulateur/solver/Tension.h b/src/Simulateur/solver/Tension.h new file mode 100644 index 0000000..ff706fa --- /dev/null +++ b/src/Simulateur/solver/Tension.h @@ -0,0 +1,245 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + + + +// Compute total current (at the anode) +double computeTotalCurrent(Parameters& parameters, MatrixXi& ID, vector& nodes, VectorXd& phi_n, double phi_anode) { + double I = 0; + vector::iterator it; + for (it = nodes.begin(); it != nodes.end(); it++) { + double phi_n_id = phi_n(it->id); + if (it->nx_anode != 0) I += it->nx_anode*eThick*parameters.dL[1]*TafelAnodeJ(phi_anode-phi_n_id); + if (it->ny_anode != 0) I += it->ny_anode*eThick*parameters.dL[0]*TafelAnodeJ(phi_anode-phi_n_id); + } + return I; +} + +// Non linear part of potential equation (Tafel) +VectorXd evaluateF(Parameters& parameters, MatrixXi& ID, vector& nodes, VectorXd& phi_n, double phi_anode, double phi_cathode) { + double dxSigmaInv = 1.0/(sigma*parameters.dL[0]), dySigmaInv = 1.0/(sigma*parameters.dL[1]); + + VectorXd F_n(phi_n.size()); F_n.setZero(); + vector::iterator it; + for (it = nodes.begin(); it != nodes.end(); it++) { + int id = it->id; + double phi_n_id = phi_n(it->id); + if (it->nx_anode != 0) F_n(id) -= it->nx_anode*it->denomInv*dxSigmaInv*TafelAnodeJ(phi_anode-phi_n_id); + if (it->ny_anode != 0) F_n(id) -= it->ny_anode*it->denomInv*dySigmaInv*TafelAnodeJ(phi_anode-phi_n_id); + if (it->nx_cathode != 0) F_n(id) -= it->nx_cathode*it->denomInv*dxSigmaInv*TafelCathodeJ(phi_n_id-phi_cathode); + if (it->ny_cathode != 0) F_n(id) -= it->ny_cathode*it->denomInv*dySigmaInv*TafelCathodeJ(phi_n_id-phi_cathode); + } + return F_n; +} +// Jacobian of the non linear part of potential equation (Tafel) +SpMat evaluateDF(Parameters& parameters, MatrixXi& ID, vector& nodes, VectorXd& phi_n, double phi_anode, double phi_cathode) { + double dxSigmaInv = 1.0/(sigma*parameters.dL[0]), dySigmaInv = 1.0/(sigma*parameters.dL[1]); + + SpMat dF_n(phi_n.size(), phi_n.size()); + vector coeffs; + + vector::iterator it; + for (it = nodes.begin(); it != nodes.end(); it++) { + double phi_n_id = phi_n(it->id), dF_n_id = 0; + if (it->nx_anode != 0) dF_n_id -= it->nx_anode*it->denomInv*dxSigmaInv*dTafelAnodeJ(phi_anode-phi_n_id); + if (it->ny_anode != 0) dF_n_id -= it->ny_anode*it->denomInv*dySigmaInv*dTafelAnodeJ(phi_anode-phi_n_id); + if (it->nx_cathode != 0) dF_n_id -= it->nx_cathode*it->denomInv*dxSigmaInv*dTafelCathodeJ(phi_n_id-phi_cathode); + if (it->ny_cathode != 0) dF_n_id -= it->ny_cathode*it->denomInv*dySigmaInv*dTafelCathodeJ(phi_n_id-phi_cathode); + + if (dF_n_id > 0) coeffs.push_back(Trip(it->id, it->id, dF_n_id)); + } + dF_n.setFromTriplets(coeffs.begin(), coeffs.end()); + return dF_n; +} + +// Calcul le potentiel en chaque noeuds de la solution électrolytique +// Pour la condition : anode = +phi_anode [V], cathode = -phi_cathode [V] +VectorXd subSolveTension(Parameters& parameters, SpMat& A, MatrixXi& ID, vector& nodes, double phi_anode, double phi_cathode, int* failFlag) { + VectorXd phi_n(A.cols()); phi_n.setZero(); + VectorXd residual, delta_phi; + + int iteration = 0; + double res0 = -1; + while (1) { + // Newton-Rahpson solving + VectorXd F_n = evaluateF(parameters, ID, nodes, phi_n, phi_anode, phi_cathode); + SpMat dF_n = evaluateDF(parameters, ID, nodes, phi_n, phi_anode, phi_cathode); + residual = -(A*phi_n + F_n); + + // Compute and check residual + double res = residual.norm()/phi_anode; + if (res0 == -1) res0 = res; + progressRefresh2a(); + if (res < TOLERANCE_POTENTIAL) break; + + // Start new iteration + if (iteration > SOLVER_MAX_ITERATION) { + *failFlag = -1; + return {}; + } iteration++; + + // Next step + SpMat Jacobian = A + dF_n; + *failFlag = solveLU(delta_phi, Jacobian, residual); + if (*failFlag != 0) return {}; + + phi_n = phi_n + delta_phi; + } + return phi_n; +} + +// Calcul le potentiel en chaque noeuds de la solution électrolytique +// En respectant la condition : int(J.dS) = parameters.I +MatrixTension solveTension(Parameters& parameters, MatrixXi& MASK, bool ignoreFail) { + int* N = parameters.N; + + // Parcours de cellules du mask et attribution d'id seulement dans l'électrolyte + MatrixXi ID(N[0], N[1]); ID.setZero(); + int n = 0; + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) { + if (MASK(i, j) == 0) ID(i, j) = n++; + } + double ddxInv = 1.0/mathSQ(parameters.dL[0]), ddyInv = 1.0/mathSQ(parameters.dL[1]); + + // Construction de la matrice des résistance (partie linéaire) et récupération des indices à l'interface electrolyte/electrode + vector coeffs; + vector nodes; + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) if (MASK(i,j) == 0) { + int id = ID(i,j); + ElectrodeNode node = {id, 0, 0, 0, 0, 0}; + bool electrolyteNeighboor[4] = {false, false, false, false}; + double denom = 0; + + // i-1 + if (i > 0) { + if (MASK(i-1, j) == 0) { + electrolyteNeighboor[0] = true; + denom += ddxInv; + } else if (MASK(i-1, j) == 1) node.nx_anode++; + else if (MASK(i-1, j) == -1) node.nx_cathode++; + } + + // i+1 + if (i < (N[0]-1)) { + if (MASK(i+1, j) == 0) { + electrolyteNeighboor[1] = true; + denom += ddxInv; + } else if (MASK(i+1, j) == 1) node.nx_anode++; + else if (MASK(i+1, j) == -1) node.nx_cathode++; + } + + // j-1 (with repeat condition) + int jm = (j > 0) ? j-1 : (N[1]-1); + if (MASK(i, jm) == 0) { + electrolyteNeighboor[2] = true; + denom += ddyInv; + } else if (MASK(i, jm) == 1) node.ny_anode++; + else if (MASK(i, jm) == -1) node.ny_cathode++; + + // j+1 (with repeat condition) + int jp = (j < (N[1]-1)) ? j+1 : 0; + if (MASK(i, jp) == 0) { + electrolyteNeighboor[3] = true; + denom += ddyInv; + } else if (MASK(i, jp) == 1) node.ny_anode++; + else if (MASK(i, jp) == -1) node.ny_cathode++; + + // i,j + coeffs.push_back(Trip(id, id, 1)); + double denomInv = 1.0/denom; + if (electrolyteNeighboor[0]) coeffs.push_back(Trip(id, ID(i-1,j), -ddxInv*denomInv)); + if (electrolyteNeighboor[1]) coeffs.push_back(Trip(id, ID(i+1,j), -ddxInv*denomInv)); + if (electrolyteNeighboor[2]) coeffs.push_back(Trip(id, ID(i,jm), -ddyInv*denomInv)); + if (electrolyteNeighboor[3]) coeffs.push_back(Trip(id, ID(i,jp), -ddyInv*denomInv)); + + if ((node.nx_anode != 0) || (node.ny_anode != 0) || (node.nx_cathode != 0) || (node.ny_cathode != 0)) { + // Un noeud ne peut être en contact à la fois avec l'anode et la cathode + if (((node.nx_anode != 0) || (node.ny_anode != 0)) && ((node.nx_cathode != 0) || (node.ny_cathode != 0))) { + errorProgress("electrical mesh too high (node in contact with both anode and cathode)"); + exit(1); + } + // Stockage du noeud frontière pour la partie non linéaire (Tafel) + node.denomInv = denomInv; + nodes.push_back(node); + } + } + + // Build linear matrix of resistors + SpMat A(n, n); + A.setFromTriplets(coeffs.begin(), coeffs.end()); + + // Newton-Rahpson iterative solve + int failFlag = 0; + // First voltage approximation + double phi_anode_cathode = parameters.J*parameters.L[0]/sigma + tafelA0[ANODE]*log(parameters.J/tafelJ0[ANODE]) + tafelA0[CATHODE]*log(parameters.J/tafelJ0[CATHODE]); + double d_phi_anode_cathode = DERIVATIVE_VARIATION*phi_anode_cathode; + + VectorXd phi; + double I = 0; + int iteration = 0; + double res0 = -1; + while (1) { + // Compute potential and current for delta_phi + phi = subSolveTension(parameters, A, ID, nodes, phi_anode_cathode, 0.0, &failFlag); + if (failFlag != 0) break; + I = computeTotalCurrent(parameters, ID, nodes, phi, phi_anode_cathode); + + // Check residual and post progress + double res = mathAbs(I-parameters.I)/parameters.I; + if (res0 == -1) res0 = res; + progressRefresh2b(mathMinMax(log(res0/res)/log(res0/TOLERANCE_CURRENT), 0.0, 1.0)); + if (res < TOLERANCE_CURRENT) break; + + // Start new iteration + if (iteration > SOLVER_MAX_ITERATION) { + failFlag = -2; + break; + } iteration++; + + // Compute derivative + phi = subSolveTension(parameters, A, ID, nodes, phi_anode_cathode-d_phi_anode_cathode, 0.0, &failFlag); + if (failFlag != 0) break; + double dI = I - computeTotalCurrent(parameters, ID, nodes, phi, phi_anode_cathode-d_phi_anode_cathode); + + // Next step (Newton) + phi_anode_cathode = phi_anode_cathode - d_phi_anode_cathode*(I-parameters.I)/(dI); + } + + // Return + MatrixXd TENSION(N[0], N[1]); + if (failFlag == 0) { + // Solving success : extract result + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) { + if (MASK(i,j) == 0) TENSION(i,j) = phi(ID(i,j)); + else if (MASK(i,j) == 1) TENSION(i,j) = phi.maxCoeff(); + else if (MASK(i,j) == -1) TENSION(i,j) = phi.minCoeff(); + else if (MASK(i,j) == -2) TENSION(i,j) = -1.0; + else if (MASK(i,j) == -3) TENSION(i,j) = phi.maxCoeff(); + else if (MASK(i,j) == -4) TENSION(i,j) = phi.minCoeff(); + } + // Compute resistor (for effective surface S = e x Ly) + double J = I/(eThick*parameters.L[1]); + double resistor = (phi.maxCoeff() - phi.minCoeff()) / J; + + return {TENSION, phi.minCoeff(), phi.maxCoeff(), phi_anode_cathode, 0, resistor, 0, true}; + } else if (ignoreFail) { + return {TENSION, 0, 0, 0, 0, 0, 0, false}; + } else { + // Solving fail : stop program + errorProgress("electrical solving error (code: " + std::to_string(failFlag) + ")"); + exit(1); + } +} diff --git a/src/Simulateur/solver/Timestep.h b/src/Simulateur/solver/Timestep.h new file mode 100644 index 0000000..725f171 --- /dev/null +++ b/src/Simulateur/solver/Timestep.h @@ -0,0 +1,44 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + + + +TimeSteps updateTimestep(TimeSteps& timesteps, double maxRadius, double maxVelocity, vector& Sites, bool* stopFlag) { + // Nucleation growing rate + double growingFactorMax = 0; + for (vector::iterator it = Sites.begin(); it != Sites.end(); it++) + growingFactorMax = MAX(growingFactorMax, it->growth); + double TgrowthMin = RFritzMin/growingFactorMax; + double dtMaximal1 = TgrowthMin/GROWING_TIME_COEF; + + // Equilibrium buoyancy velocity + double UElec = abs(bPrefactor*mathSQ(maxRadius)); + double dtElec = 2*maxRadius/(ELECTRICAL_TIME_COEF*UElec); + + // Check Reynolds + double U = abs(bPrefactor*mathSQ(maxRadius)); + double dtMaximal2 = 2*maxRadius/(BUOYANCY_TIME_COEF*U); + double reynolds = (rho/mu)*maxVelocity*maxRadius; + if (reynolds > 10) { + errorProgress("reynolds number too high (Re=" + std::to_string(reynolds) + ")"); + *stopFlag = true; + } + + // Return the smallest value + double dt = MIN(dtMaximal1, dtMaximal2); + + return {dt, dtElec, timesteps.i, timesteps.iElec}; +} diff --git a/src/Simulateur/solver/Topology.h b/src/Simulateur/solver/Topology.h new file mode 100644 index 0000000..53c427c --- /dev/null +++ b/src/Simulateur/solver/Topology.h @@ -0,0 +1,298 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +/* --------- Mask generation ---------- */ +/* ------------------------------------ */ +int checkTopologyPolySub(Topology& topology, double x, double y) { + Point X = {x, y}; + for (int s = 0; s < 2; s++) for (int p = 0; p < topology.polygonCount[s]; p++) { + int n = topology.polygonIndex[s][p+1] - topology.polygonIndex[s][p]; + Point *polys = new Point[n]; + for (int seg = 0; seg < n; seg++) { + int index = topology.polygonIndex[s][p] + seg; + polys[seg] = {topology.polygonDefinitions[index].x, topology.polygonDefinitions[index].y}; + } + if (isInside(polys, n, X)) { + return (s == 0) ? 1 : -1; + } + delete [] polys; + } + return 0; +} +int checkTopologyPoly(int* N, Topology& topology, double x, double y) { + double dY = 0.66666 / (N[1] - 1); + double r = 0; + if (r == 0) r = checkTopologyPolySub(topology, x, y + dY); + if (r == 0) r = checkTopologyPolySub(topology, x, y - dY); + return r; +} + +Eigen::MatrixXi generateMask(Parameters& parameters, Topology& topology) { + int* N = parameters.N; + Eigen::MatrixXi MASK(N[0], N[1]); + // Fill mask + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) { + double x, y; indexToPosition(N, i, j, &x, &y); + MASK(i, j) = checkTopologyPoly(N, topology, x, y); // anode(1) electrolyte(0) cathode(-1) insulation(-2) anode_trapped(-3) cathode_trapped(-4) + } + return MASK; +} + +void getAttachedBubbleLocation(double& x, double& y, double& r, list::iterator& it) { + // Calcul du rayon actuel + r = it->r + it->site->tmiss*it->site->growth; + double dLMouillage = r*sin(M_PI - theta[it->site->type]); + // Positionnement réel d'une bulle attachée + x = it->site->x + it->site->nx*r - it->site->nx*dLMouillage; + y = it->site->y + it->site->ny*r - it->site->ny*dLMouillage; +} + +void subupdateMask(Parameters& parameters, MatrixXi& MASK, list::iterator& it, bool attached) { + int* N = parameters.N; + double* L = parameters.L; + + double itx = it->x, ity = it->y, itr = it->r; + if (attached) getAttachedBubbleLocation(itx, ity, itr, it); + + int i, j; positionToIndex(N, &i, &j, itx/L[0], ity/L[1]); + double x, y; + + + int di = round(itr/parameters.dL[0]), dj = round(itr/parameters.dL[1]); + for (int ii = -di; ii <= di; ii++) for (int jj = -dj; jj <= dj; jj++) { + indexToPosition(N, i+ii, j+jj, &x, &y); + if ((mathSQ((L[0]*x)-itx) + mathSQ((L[1]*y)-ity)) < mathSQ(itr)) { + if (isBetween(i+ii, j+jj, N)) if (MASK(i+ii, j+jj) == 0) + MASK(i+ii, j+jj) = -2; + } + } +} + +MatrixXi updateMask(Parameters& parameters, MatrixXi MASK, ListBubbles& Bubbles) { + // Replace all insulated and others specials with electrolyte + int* N = parameters.N; + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) { + if ((MASK(i, j) == -2) || (MASK(i, j) == -3) || (MASK(i, j) == -4)) { + MASK(i, j) = 0; + } + } + // Add bubbles + list::iterator it; + for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) { + subupdateMask(parameters, MASK, it, true); + } + for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + subupdateMask(parameters, MASK, it, false); + } + + // Remove insulated areas + // Les bulles peuvent rendre certains noeuds d'électrolyte isolés (lapluspart du temps contre une électrode) + // C'est noeuds sont alors au potentiel de l'électrode en question (J=0), mais il est nécessaire de les détecter + MatrixXi VISITED(N[0], N[1]); VISITED.setZero(); + for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) { + // Vérifier qu'un noeud d'électrolyte relie bien l'anode à la cathode + if (MASK(i,j) == 0) if (VISITED(i,j) == 0) { + bool anode = false, cathode = false; + list toCheck = {{i,j}}; + VISITED(i,j) = 1; + vector checked; + while (toCheck.size() > 0) { + // Vérifie l'élement et propage + list::iterator coord = toCheck.begin(); + int ii = coord->i, jj = coord->j; + // voisin x- + if (ii > 0) { + if (MASK(ii-1, jj) == 0) if (VISITED(ii-1, jj) == 0) { + toCheck.push_back({ii-1, jj}); + VISITED(ii-1,jj) = 1; + } + if (MASK(ii-1, jj) == 1) anode = true; + if (MASK(ii-1, jj) == -1) cathode = true; + } + // voisin x+ + if (ii < (N[0]-1)) { + if (MASK(ii+1, jj) == 0) if (VISITED(ii+1, jj) == 0) { + toCheck.push_back({ii+1, jj}); + VISITED(ii+1,jj) = 1; + } + if (MASK(ii+1, jj) == 1) anode = true; + if (MASK(ii+1, jj) == -1) cathode = true; + } + // voisin y- + int jjm = (jj > 0) ? jj-1 : (N[1]-1); + if (MASK(ii, jjm) == 0) if (VISITED(ii, jjm) == 0) { + toCheck.push_back({ii, jjm}); + VISITED(ii,jjm) = 1; + } + if (MASK(ii, jjm) == 1) anode = true; + if (MASK(ii, jjm) == -1) cathode = true; + // voisin y+ + int jjp = (jj < (N[1]-1)) ? jj+1 : 0; + if (MASK(ii, jjp) == 0) if (VISITED(ii, jjp) == 0) { + toCheck.push_back({ii, jjp}); + VISITED(ii,jjp) = 1; + } + if (MASK(ii, jjp) == 1) anode = true; + if (MASK(ii, jjp) == -1) cathode = true; + + // valider l'élement testé + toCheck.erase(coord); + checked.push_back({ii, jj}); + } + // Si le groupe ne connecte pas l'anode à la caothde... + if (!(anode && cathode)) { + for (vector::iterator it = checked.begin(); it != checked.end(); it++) { + if (anode && (!cathode)) MASK(it->i, it->j) = -3; // Il sera au même potentiel que l'anode + else if ((!anode) && cathode) MASK(it->i, it->j) = -4; // ... que la cathode + else MASK(it->i, it->j) = -2; // ou completement indéterminé + } + } + } + } + return MASK; +} +/* ------------------------------------ */ + + + +/* ------ Generate wall list -------- */ +/* ---------------------------------- */ +vector generateWalls(Parameters& parameters, Topology& topology) { + list WallList; + for (int s = 0; s < 2; s++) for (int p = 0; p < topology.polygonCount[s]; p++) { + int n = topology.polygonIndex[s][p+1] - topology.polygonIndex[s][p]; + for (int seg = 0; seg < n; seg++) { + int index1 = topology.polygonIndex[s][p] + seg; + int index2 = topology.polygonIndex[s][p] + ((seg+1)%n); + Point P1 = {parameters.L[0]*topology.polygonDefinitions[index1].x, parameters.L[1]*topology.polygonDefinitions[index1].y}; + Point P2 = {parameters.L[0]*topology.polygonDefinitions[index2].x, parameters.L[1]*topology.polygonDefinitions[index2].y}; + Vect u = {P2.x-P1.x, P2.y-P1.y}; + double length = vectLength(u); + u.x = u.x / length; u.y = u.y / length; + Vect n = {u.y, -u.x}; + + WallList.push_back({{P1.x, P1.y}, {P2.x, P2.y}, u, n, length}); + WallList.push_back({{P1.x, P1.y-parameters.L[1]}, {P2.x, P2.y-parameters.L[1]}, u, n, length}); + WallList.push_back({{P1.x, P1.y+parameters.L[1]}, {P2.x, P2.y+parameters.L[1]}, u, n, length}); + } + } + vector Walls(WallList.begin(), WallList.end()); + return Walls; +} +/* ---------------------------------- */ + + + +/* --------- Wall distance ---------- */ +/* ---------------------------------- */ +WallIndicator checkWall(Wall& wall, double x, double y, double maxDistance) { + // Check if is inside + Vect P1P = {x - wall.P1.x, y - wall.P1.y}; + double pos = vectProject(P1P, wall.u); + if ((pos < -maxDistance) || (pos > wall.length+maxDistance)) + return {{0.0, 0.0}, {0.0, 0.0}, -1.0, nullptr}; + + // Test wall side + double dist; Vect n; + if (pos < 0.0) { + dist = vectLength(P1P); + if (dist > maxDistance) + return {{0.0, 0.0}, {0.0, 0.0}, -1.0, nullptr}; + n = {P1P.x/dist, P1P.y/dist}; + return {wall.P1, n, dist, nullptr}; + } else if (pos > wall.length) { + Vect P2P = {x - wall.P2.x, y - wall.P2.y}; + dist = vectLength(P2P); + if (dist > maxDistance) + return {{0.0, 0.0}, {0.0, 0.0}, -1.0, nullptr}; + n = {P2P.x/dist, P2P.y/dist}; + return {wall.P2, n, dist, nullptr}; + } else { + dist = vectProject(P1P, wall.n); + if ((dist < 0.0) || (dist > maxDistance)) + return {{0.0, 0.0}, {0.0, 0.0}, -1.0, nullptr}; + n = wall.n; + return {{wall.P1.x + pos*wall.u.x, wall.P1.y + pos*wall.u.y}, n, dist, &wall}; + } +} + + +void checkWalls(WallIndicator* indicators, vector& Walls, double x, double y, double maxDistance) { + for (int i = 0; i < wallMaxmimumIndicators; i++) { + indicators[i].P = {-1, -1}; + indicators[i].distance = INF; + }//return; + + for (vector::iterator it = Walls.begin(); it != Walls.end(); it++) { + // Test + WallIndicator indicator = checkWall(*it, x, y, maxDistance); + if (indicator.distance == -1.0) + continue; + + // Une parois est détectée comme à portée... + // Si c'est un point obtenu par linéarisation (&wall = nullptr) : + // il ne doit ni déjà être dans la liste ni faire partie des extremums d'une parois déjà ajoutée + if (indicator.wall == nullptr) { + bool removeFlag = false; + for (int i = 0; i < wallMaxmimumIndicators; i++) if (indicators[i].distance != INF) { + if (indicators[i].wall == nullptr) { + if (pointEquals(indicator.P, indicators[i].P)) { + removeFlag = true; + break; + } + } else { + if (pointEquals(indicator.P, indicators[i].wall->P1) || pointEquals(indicator.P, indicators[i].wall->P2)) { + removeFlag = true; + break; + } + } + } + if (removeFlag) continue; // Rejeté ! + // Sinon, il convient de vérifier que les points extremums de la parois à ajouter ne font pas l'objet d'un point obtenu par linéarisation + // Cependant, il peut y avoir plusieurs fois la parois (sur-contrainte logique) + } else { + for (int i = 0; i < wallMaxmimumIndicators; i++) if (indicators[i].distance != INF) { + if (indicators[i].wall == nullptr) { + if (pointEquals(indicator.wall->P1, indicators[i].P) || pointEquals(indicator.wall->P2, indicators[i].P)) { + // Retrait et décallage + for (int j = i; j < (wallMaxmimumIndicators-1); j++) { + indicators[j].P = indicators[j+1].P; + indicators[j].n = indicators[j+1].n; + indicators[j].distance = indicators[j+1].distance; + indicators[j].wall = indicators[j+1].wall; + } + indicators[wallMaxmimumIndicators-1].P = {-1, -1}; + indicators[wallMaxmimumIndicators-1].distance = INF; + i--; + } + } + } + } + + // Save the result + for (int i = 0; i < wallMaxmimumIndicators; i++) { + if (indicator.distance < indicators[i].distance) { + indicators[i].P = indicator.P; + indicators[i].distance = indicator.distance; + indicators[i].n = indicator.n; + indicators[i].wall = indicator.wall; + if (i == (wallMaxmimumIndicators-1)) maxDistance = indicator.distance; + break; + } + } + } +} +/* ---------------------------------- */ diff --git a/src/Simulateur/solver/utils/fundyn.h b/src/Simulateur/solver/utils/fundyn.h new file mode 100644 index 0000000..027339a --- /dev/null +++ b/src/Simulateur/solver/utils/fundyn.h @@ -0,0 +1,182 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +/* Dynamics functions */ + +/* + * Sphere-sphere long-range ------------------------------------ + */ + +// Mobility matrix parts for sphere-sphere longe-range interactions +void setMobilityMatrixSphere_I(double* M) { + M[0] = 1; M[1] = 0; + M[2] = 0; M[3] = 1; +} +void setMobilityMatrixSphere_J(double* M, double coeff, double sqRa, double sqRb, double rxab, double ryab) { + double r_inv, r3_inv, r5_inv, part1, part2; + r_inv = 1.0/sqrt(rxab*rxab + ryab*ryab); r3_inv = r_inv*r_inv*r_inv; r5_inv = r3_inv*r_inv*r_inv; + + double Sab2 = (sqRa + sqRb) / 3.0; + part1 = (r_inv + Sab2*r3_inv); part2 = (r3_inv - 3.0*Sab2*r5_inv); + + M[0] = coeff*(part1 + part2*rxab*rxab); M[1] = coeff*part2*rxab*ryab; + M[2] = coeff*part2*rxab*ryab; M[3] = coeff*(part1 + part2*ryab*ryab); +} + +void addMobilityMatrixPlane(double* M, Bubble* it) { + for (int kk = 0; kk < wallMaxmimumIndicators; kk++) { + if (it->interactor.walls[kk].distance < (0.5*it->r*interactionLongMaximumRange)) { + double Ra = it->r; + double da = it->interactor.walls[kk].distance; + double daInv = 1.0/da; + + double nx = it->interactor.walls[kk].n.x, ny = it->interactor.walls[kk].n.y; + + double part1 = -1.125*Ra*daInv; + double part2 = -0.5625*Ra*daInv; + double part3 = (part1 - part2); + + M[0] -= part2 + part3*nx*nx; + M[1] -= part3*nx*ny; + M[2] -= part3*nx*ny; + M[3] -= part2 + part3*ny*ny; + } + } +} + + + + + +/* + * Sphere-sphere contact ------------------------------------ + */ + +// Rlub matrix parts for sphere-sphere contact interactions +void matrixRlubS(double* Rlub, bool add, double Ra, double Rb, double rxab, double ryab) { + double Reqab = Ra*Rb/(Ra+Rb); + double rabSq = (rxab*rxab + ryab*ryab); + double epsad = MAX(SMALLEST, sqrt(rabSq) - Ra - Rb)/Reqab; + double rabSqinv = 1.0/rabSq, epsadInv = 1.0/epsad; + + double part1 = 0.25*epsadInv; + double part2 = 0.166666667*log(epsadInv); + double part3 = (part1 - part2)*rabSqinv; + + #define Bab0 part2 + part3*rxab*rxab + #define Bab1 part3*rxab*ryab + #define Bab2 part3*rxab*ryab + #define Bab3 part2 + part3*ryab*ryab + if (add) { + Rlub[0] += Bab0; + Rlub[1] += Bab1; + Rlub[2] += Bab2; + Rlub[3] += Bab3; + } else { + Rlub[0] -= Bab0; + Rlub[1] -= Bab1; + Rlub[2] -= Bab2; + Rlub[3] -= Bab3; + } +} +// Sphere-sphere contact +// i +void addSphereSphereContact_I(double* Rij, double Ly, Bubble* it) { + for (int kk = 0; kk < it->interactor.count; kk++) { + Bubble *kt = it->interactor.bubbles[kk]; + if ((it->interactor.distancesSqNorm[kk] < interactionLubMaximumRangeSq) && (kt->r != 0)) { + int kMode = it->interactor.imagedModes[kk]; + matrixRlubS(Rij, false, it->r, kt->r, kt->x-it->x, yRepeat(kt->y, Ly, kMode)-it->y); + } + } +} +// j +void addSphereSphereContact_J(double* Rij, double Ly, Bubble* it, Bubble* jt, int j, int jMode) { + matrixRlubS(Rij, true, it->r, jt->r, jt->x-it->x, yRepeat(jt->y, Ly, jMode)-it->y); +} + + + + +/* + * Sphere-plan contact ------------------------------------ + */ +// D matrix parts for sphere-plan contact +void matrixRlubW(double* Rlub, double Ra, double ra, double nx, double ny) { + double epsa = MAX(SMALLEST, ra - Ra)/Ra; + double epsaInv = 1.0/epsa; + + double part1 = -epsaInv; + double part2 = -0.53333333333*log(epsaInv); + double part3 = (part1 - part2); + + Rlub[0] += part2 + part3*nx*nx; + Rlub[1] += part3*nx*ny; + Rlub[2] += part3*nx*ny; + Rlub[3] += part2 + part3*ny*ny; +} +// Sphere-plan contact +void addSpherePlanContact_I(double* Rij, Bubble* it) { + for (int kk = 0; kk < wallMaxmimumIndicators; kk++) { + if (it->interactor.walls[kk].distance < (0.5*it->r*interactionLubMaximumRange)) { + matrixRlubW(Rij, it->r, it->interactor.walls[kk].distance, it->interactor.walls[kk].n.x, it->interactor.walls[kk].n.y); + } + } +} + + + + + +/* + * Repulsive terms ------------------------------------ + */ +// Sphere-sphere +void addSphereSphereRepulsive(VectorXd& U0, Bubble* it, int i) { + for (int jj = 0; jj < it->interactor.count; jj++) { + Bubble *jt = it->interactor.bubbles[jj]; + // Check if k is in range and valid + if ((it->interactor.distancesSqNorm[jj] < interactionRepulsiveRangeSq) && (jt->r != 0)) { + double dxij = it->x-jt->x, dyij = it->y-jt->y; + double rSq = (SQ(dxij) + SQ(dyij)); + double RrefSq = SQ(interactionRepulsiveRange); + double RrirjSq = mathSQ(it->r + jt->r); + + double r = sqrt(rSq); + dxij /= r; dyij /= r; + + double F = it->r*interactionRepulsiveVelocity*mathSQ((RrefSq - rSq) / (RrefSq - RrirjSq)); + + U0[2*i+0] -= F*dxij; + U0[2*i+1] -= F*dyij; + } + } +} + + + + + +/* + * Sphere-sphere force estimation ------------------------------------ + */ +void computeSphereSphereForce(double* F, Bubble* it, Bubble* jt, int jMode, double Ly) { + double Rlub[4] = {0, 0, 0, 0}; + matrixRlubS(Rlub, false, it->r, jt->r, jt->x-it->x, yRepeat(jt->y, Ly, jMode)-it->y); + + F[0] += Rlub[0]*jt->ux + Rlub[1]*jt->uy; + F[1] += Rlub[2]*jt->ux + Rlub[3]*jt->uy; +} diff --git a/src/Simulateur/solver/utils/funinter.h b/src/Simulateur/solver/utils/funinter.h new file mode 100644 index 0000000..42540b8 --- /dev/null +++ b/src/Simulateur/solver/utils/funinter.h @@ -0,0 +1,89 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +/* ------ Test bubble intractability ------ */ +/* ---------------------------------------- */ +void testInteractibilityMode(Bubble* it, Bubble* jt, int i, int j, int modeI, int modeJ, double deltaX, double deltaY, double RavvSq, double dMax, bool* stopFlag) { + if ((mathAbs(deltaX) < dMax) && (mathAbs(deltaY) < dMax)) { + double deltaSqNorm = (deltaX*deltaX + deltaY*deltaY) / RavvSq; // Distance inter-bulle normalisée par le rayon moyen + if (deltaSqNorm < interactionLongMaximumRangeSq) { + if (it->interactor.count == interactionMaxAmountAbsolute) {*stopFlag = true; return;} + it->interactor.interactionIndexs[it->interactor.count] = j; + it->interactor.distancesSqNorm[it->interactor.count] = deltaSqNorm; + it->interactor.bubbles[it->interactor.count] = jt; + it->interactor.imagedModes[it->interactor.count] = modeJ; + it->interactor.count++; + + if (jt->interactor.count == interactionMaxAmountAbsolute) {*stopFlag = true; return;} + jt->interactor.interactionIndexs[jt->interactor.count] = i; + jt->interactor.distancesSqNorm[jt->interactor.count] = deltaSqNorm; + jt->interactor.bubbles[jt->interactor.count] = it; + jt->interactor.imagedModes[jt->interactor.count] = modeI; + jt->interactor.count++; + } + } +} +void checkBubbles(list::iterator& it, int i, ListBubbles& Bubbles, double Ly, bool* stopFlag) { + it->statCollision = 0; + int j = i; list::iterator itbis = it; + // Bulles détachées + for (itbis++; itbis != Bubbles.detached.end(); itbis++) { + j++; + double Ravv = 0.5*(it->r + itbis->r), RavvSq = SQ(Ravv); + double dMax = Ravv * interactionLongMaximumRange; + + /* mode 0 : pas de modifications */ + double deltaX = (itbis->x - it->x), deltaY = (itbis->y - it->y); + testInteractibilityMode(&*it, &*itbis, i, j, 0, 0, deltaX, deltaY, RavvSq, dMax, stopFlag); + /* mode 1 : recopie en haut */ + if (itbis->y < dMax) { + deltaX = (itbis->x - it->x), deltaY = ((itbis->y + Ly) - it->y); + testInteractibilityMode(&*it, &*itbis, i, j, 2, 1, deltaX, deltaY, RavvSq, dMax, stopFlag); + } + /* mode 2 : recopie en bas */ + if (itbis->y > (Ly - dMax)) { + deltaX = (itbis->x - it->x), deltaY = ((itbis->y - Ly) - it->y); + testInteractibilityMode(&*it, &*itbis, i, j, 1, 2, deltaX, deltaY, RavvSq, dMax, stopFlag); + } + } + // Bulles attachées + for (itbis = Bubbles.attached.begin(); itbis != Bubbles.attached.end(); itbis++) { + double Ravv = 0.5*(it->r + itbis->r), RavvSq = SQ(Ravv); + double dMax = Ravv * interactionLongMaximumRange; + + /* mode 0 : pas de modifications */ + double deltaX = (itbis->x - it->x), deltaY = (itbis->y - it->y); + testInteractibilityMode(&*it, &*itbis, i, -1, 0, 0, deltaX, deltaY, RavvSq, dMax, stopFlag); + /* mode 1 : recopie en haut */ + if (itbis->y < dMax) { + deltaX = (itbis->x - it->x), deltaY = ((itbis->y + Ly) - it->y); + testInteractibilityMode(&*it, &*itbis, i, -1, 2, 1, deltaX, deltaY, RavvSq, dMax, stopFlag); + } + /* mode 2 : recopie en bas */ + if (itbis->y > (Ly - dMax)) { + deltaX = (itbis->x - it->x), deltaY = ((itbis->y - Ly) - it->y); + testInteractibilityMode(&*it, &*itbis, i, -1, 1, 2, deltaX, deltaY, RavvSq, dMax, stopFlag); + } + } +} + +/* Correct y for repeating symmetry */ +double yRepeat(double y, double Ly, int mode) { + if (mode == 0) return y; // mode 0 + else if (mode == 1) return y+Ly; // mode 1 + else return y-Ly; // mode 2 +} +/* ---------------------------------------- */ diff --git a/src/Simulateur/solver/utils/funvelo.h b/src/Simulateur/solver/utils/funvelo.h new file mode 100644 index 0000000..e80cac7 --- /dev/null +++ b/src/Simulateur/solver/utils/funvelo.h @@ -0,0 +1,86 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + + +/* --------- Velocity correction ---------- */ +/* ---------------------------------------- */ +void bubbleVelocityCorrection(double Ly, Bubble* it, bool* stopFlag); +bool bubbleVelocitySub(double Ly, Bubble* it, bool mode, WallIndicator** indicatorDone, Bubble** bubbleDone, bool* stopFlag) { + // avec les murs + for (int kk = 0; kk < wallMaxmimumIndicators; kk++) { + WallIndicator *indicator = &it->interactor.walls[kk]; + if (!mode) if (indicator == *indicatorDone) continue; + if (indicator->distance < (it->r*interactionSafetyHalfRange)) { + double proj0 = it->ux*indicator->n.x + it->uy*indicator->n.y; + if (proj0 < 0) { // la bulle se déplace vers le mur ? + if (mode) { + it->ux = it->ux - proj0*indicator->n.x; + it->uy = it->uy - proj0*indicator->n.y; + *indicatorDone = indicator; + return true; + } else { + // blocage définitif ! -> arrêt + it->ux = 0; it->uy = 0; + } + } + } + } + // avec les autres spheres + for (int kk = 0; kk < it->interactor.count; kk++) { + if (it->interactor.bubbles[kk]->r != 0) + if (it->interactor.distancesSqNorm[kk] < interactionSafetyRangeSq) { + Bubble* kt = it->interactor.bubbles[kk]; + int kMode = it->interactor.imagedModes[kk]; + if (!mode) if (kt == *bubbleDone) continue; + if ((kt->statCollision == 0) && (it->interactor.interactionIndexs[kk] != -1)) { + // si c'est une bulle détachée non traitée, il faut lui fixer sa vitesse avant + bubbleVelocityCorrection(Ly, kt, stopFlag); + } + double P21x = (it->x - kt->x), P21y = (it->y - yRepeat(kt->y, Ly, kMode)); + double proj1 = it->ux*P21x + it->uy*P21y; + if (proj1 < 0) { // la bulle considérée se déplace vers une autre ? + double proj2 = 0; + if (kt->statCollision == 2) proj2 = kt->ux*P21x + kt->uy*P21y; + if (proj1 < proj2) { // le mouvement relatif est-il bien un rapprochement ? + if (mode) { + double proj = (proj1-proj2)/(P21x*P21x + P21y*P21y); + it->ux = it->ux - proj*P21x; + it->uy = it->uy - proj*P21y; + *bubbleDone = kt; + return true; + } else { + double k = MAX(0, proj2 / proj1); + it->ux = k*it->ux; + it->uy = k*it->uy; + } + } + } + } + } + return false; +} +void bubbleVelocityCorrection(double Ly, Bubble* it, bool* stopFlag) { + if (it->statCollision != 0) return; + it->statCollision = 1; + + WallIndicator* indicatorDone = nullptr; Bubble* bubbleDone = nullptr; + if (bubbleVelocitySub(Ly, it, true, &indicatorDone, &bubbleDone, stopFlag)) { + bubbleVelocitySub(Ly, it, false, &indicatorDone, &bubbleDone, stopFlag); + } + + it->statCollision = 2; +} +/* ---------------------------------------- */ diff --git a/src/Simulateur/utils/loader.h b/src/Simulateur/utils/loader.h new file mode 100644 index 0000000..0af5e16 --- /dev/null +++ b/src/Simulateur/utils/loader.h @@ -0,0 +1,48 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +void readConfigurationTopology(Configuration& configuration, Topology& topology, double* args) { + int pos = 0, count; + /* -------- Parameters -------- */ + // System dimensions [m] + configuration.L[0] = args[pos++]; + configuration.L[1] = args[pos++]; + // Current density [A/m2] + configuration.J = args[pos++]; + // Simulation duration [s] + configuration.tf = args[pos++]; + /* --------- Topology --------- */ + // Polygon count + topology.polygonCount[0] = args[pos++]; + topology.polygonCount[1] = args[pos++]; + // Polygon indexs + count = topology.polygonCount[0]+1; + topology.polygonIndex[0] = new int[count]; + for (int i = 0; i < count; i++) { + topology.polygonIndex[0][i] = args[pos++]; + } + count = topology.polygonCount[1]+1; + topology.polygonIndex[1] = new int[count]; + for (int i = 0; i < count; i++) { + topology.polygonIndex[1][i] = args[pos++]; + } + // Polygon defs + count = topology.polygonIndex[1][topology.polygonCount[1]]; + topology.polygonDefinitions = new Point[count]; + for (int i = 0; i < count; i++) { + topology.polygonDefinitions[i] = {args[pos++], args[pos++]}; + } +} diff --git a/src/Simulateur/utils/math2.h b/src/Simulateur/utils/math2.h new file mode 100644 index 0000000..1e12ea0 --- /dev/null +++ b/src/Simulateur/utils/math2.h @@ -0,0 +1,174 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +#define MAX(a,b) ((a)>(b) ? (a) : (b)) +#define MIN(a,b) ((a)<(b) ? (a) : (b)) +#define ABS(x) ((x) < 0.0 ? (-(x)) : (x)) +#define SQ(x) ((x)*(x)) +#define CUBE(x) ((x)*(x)*(x)) + +#define SMALLEST 5.0e-6 +#define INF 100000.0 +#define INF_NORMALIZED 100.0 + + +/* ------ Basics ------ */ +/* -------------------- */ +double mathAbs(double x) { + return (x < 0.0) ? -x : x; +} +double mathSQ(double x) { + return x*x; +} +double mathCB(double x) { + return x*x*x; +} +int mathSign(double x) { + return (x == 0) ? 0 : ((x > 0) ? +1 : -1); +} +double mathMax(double a, double b) { + return (a > b) ? a : b; +} +double mathMin(double a, double b) { + return (a < b) ? a : b; +} +double mathMinMax(double x, double min, double max) { + return (x < min) ? min : ((x > max) ? max : x); +} +double mathMap(double x, double xMin, double xMax, double yMin, double yMax) { + return yMin + (x-xMin)*(yMax-yMin)/(xMax-xMin); +} +/* -------------------- */ + +/* ------ Random ------ */ +/* -------------------- */ +bool mathRandInitialized = false; +void mathRandInit() { + if (!mathRandInitialized) { + srand(0); // Constant seed for easy debug + mathRandInitialized = true; + } +} +double mathRand() { + return ((double) rand() / (RAND_MAX)); +} +/* -------------------- */ + + +/* ---- Vectors and co ---- */ +/* ------------------------ */ +bool pointEquals(Point& P1, Point& P2) { + return ((P1).x == (P2).x) && ((P1).y == (P2).y); +} +double vectLength(Vect& u) { + return sqrt(mathSQ(u.x) + mathSQ(u.y)); +} +double vectProject(Vect& u, Vect& v) { + return u.x*v.x + u.y*v.y; +} +Point vectForward(Point& x0, Vect& u, double pos) { + return {x0.x + (u.x*pos), x0.y + (u.y*pos)}; +} +/* ------------------------ */ + + + +/* ---- Index <-> position ---- */ +/* ---------------------------- */ +void indexToPosition(int* N, int i, int j, double *x, double *y) { + *x = ((double) i)/(N[0]-1); + *y = ((double) j)/(N[1]-1); +} +void positionToIndex(int* N, int *i, int *j, double x, double y) { + *i = round((N[0]-1)*x); + *j = round((N[1]-1)*y); +} +int indexToK(int* N, int i, int j) { + return i + (N[0]*j); +} +/* ---------------------------- */ + + + + +/* --- Linear system solver --- */ +/* ---------------------------- */ +int solveLU(VectorXd& x, SpMat& A, VectorXd& b) { + // Decompose LU + Eigen::SparseLU solver; + solver.compute(A); + + // Decomposition failed + if(solver.info() != Eigen::Success) return 1; + + // Solve + x = solver.solve(b); + + // Solving failed + if(solver.info() != Eigen::Success) return 2; + + // Check tolerance + double residual = (((A*x)-b).norm() / b.norm()); + if (residual > TOLERANCE_LU) return 3; + + return 0; +} +int invertLU(SpMat& Ainv, SpMat& A) { + // Decompose LU + Eigen::SparseLU solver; + solver.compute(A); + + // Decomposition failed + if(solver.info() != Eigen::Success) return 1; + + // Solve + SpMat I = Eigen::MatrixXd::Identity(A.rows(), A.cols()).sparseView(); + Ainv = solver.solve(I); + + // Solving failed + if(solver.info() != Eigen::Success) return 2; + + // Check tolerance + double residual = (A*Ainv-I).norm() / I.norm(); + if (residual > TOLERANCE_LU) return 3; + + return 0; +} +/* ---------------------------- */ + + + + +/* ---------- Others ---------- */ +/* ---------------------------- */ +double bounds(double x, double min, double max) { + double range = max - min; + while (x >= max) + x -= range; + while (x < min) + x += range; + return x; +} +int isApprox(double a, double b, double tol) { + return ((a > b*(1.0-tol)) && (a < b*(1.0+tol))); +} +bool isBetween(int i, int iMin, int iMax) { + return (i >= iMin) && (i <= iMax); +} +bool isBetween(int i, int j, int* N) { + return isBetween(i, 0, N[0]-1) && isBetween(j, 0, N[1]-1); +} +/* ---------------------------- */ diff --git a/src/Simulateur/utils/output.h b/src/Simulateur/utils/output.h new file mode 100644 index 0000000..bceafe2 --- /dev/null +++ b/src/Simulateur/utils/output.h @@ -0,0 +1,309 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +// Delete output content +namespace fs = std::filesystem; +void outputDelete(const char *folder) { + for (const auto& entry : fs::directory_iterator(folder)) { + if (entry.is_regular_file()) { + fs::remove(entry.path()); + } else if (entry.is_directory()) { + outputDelete(entry.path().string().c_str()); + fs::remove(entry.path()); + } + } +} +// Create folder before calling loadBMP API +unsigned int dispatch_loadbmp_encode_file(const char *filename, const unsigned char *imageData, unsigned int width, unsigned int height, unsigned int components) { + fs::path filePath(filename); + string parentFolder = filePath.parent_path().string(); + if (!fs::exists(parentFolder)) fs::create_directory(parentFolder); + + return loadbmp_encode_file(filename, imageData, width, height, components); +} + + + +// Text output +void outputData(const char *filename, Parameters& parameters, MatrixXd& MATRIX, double minVal, double maxVal) { + ofstream frameFile; + frameFile = ofstream(filename); + if (!frameFile.is_open()) { + cerr << "Unable to create outputfile." << endl; + exit(1); + } + frameFile << "Data:" << endl; + + int width = MATRIX.rows(), height = MATRIX.cols(); + double x, y; + for (int j = 0; j < (int) height; j++) for (int i = 0; i < (int) width; i++) { + int jj = height - j - 1; + indexToPosition(parameters.N, i, jj, &x, &y); + double value = MATRIX(i,j); + frameFile << parameters.L[0]*x << "," << parameters.L[1]*y << ","; + if ((value == minVal) || (value == maxVal)) frameFile << "NaN"; + else frameFile << value; + frameFile << endl; + } +} +void outputData(string folder, const char *subfilename, Parameters& parameters, MatrixXd& MATRIX, double minVal, double maxVal) { + string filename = folder + string(subfilename); + outputData(filename.c_str(), parameters, MATRIX, minVal, maxVal); +} + + +// Bitmap output +void outputMatrix(const char *filename, MatrixXd& MATRIX, double minVal, double maxVal, double nullVal) { + int width = MATRIX.rows(), height = MATRIX.cols(); + unsigned char *pixels = (unsigned char *) calloc(3*width*height, sizeof (unsigned char)); + unsigned int err; + + for (int i = 0; i < (int) width; i++) for (int j = 0; j < (int) height; j++) { + int jj = height - j - 1; + double v = (MATRIX(i,j) - minVal)/(maxVal-minVal); + + if (MATRIX(i,j) == nullVal) { + pixels[3*(i+width*jj)+0] = 34; + pixels[3*(i+width*jj)+1] = 177; + pixels[3*(i+width*jj)+2] = 76; + } else if (v <= 0.0) { + pixels[3*(i+width*jj)+0] = 255; + pixels[3*(i+width*jj)+1] = 255; + pixels[3*(i+width*jj)+2] = 255; + } else if (v >= 1.0) { + pixels[3*(i+width*jj)+0] = 255; + pixels[3*(i+width*jj)+1] = 255; + pixels[3*(i+width*jj)+2] = 255; + } else if ((v > 0) && (v < 1.0)) { + pixels[3*(i+width*jj)+0] = mathMin(mathMax(0, mathMap(v, 0.5, 1.0, 0, 255)), 255); + pixels[3*(i+width*jj)+2] = mathMin(mathMax(0, mathMap(v, 0.0, 0.5, 255, 0)), 255); + + for (int c = 0; c < 3; c++) { + pixels[3*(i+width*jj)+c] = mathMin(mathMax(0, ((int) pixels[3*(i+width*jj)+c]) + ((int) mathMin(mathMax(0, 255-mathAbs(mathMap(v, 0.25, 0.75, -255, 255))), 255))), 255); + } + } + } + + err = dispatch_loadbmp_encode_file(filename, pixels, width, height, LOADBMP_RGB); + if (err) std::cout << "FAILURE: " << EXIT_FAILURE << std::endl; + + free(pixels); +} +void outputMatrix(string folder, const char *subfilename, MatrixXd& MATRIX, double minVal, double maxVal, double nullVal) { + string filename = folder + string(subfilename); + outputMatrix(filename.c_str(), MATRIX, minVal, maxVal, nullVal); +} +void outputMatrix(string folder, const char *subfolder, int i, MatrixXd& MATRIX, double minVal, double maxVal, double nullVal) { + string filename = folder + string(subfolder) + string("frame_") + std::to_string(i) + string(".bmp"); + outputMatrix(filename.c_str(), MATRIX, minVal, maxVal, nullVal); +} + + +void outputMatrixI(const char *filename, MatrixXi& MATRIX, double minVal, double maxVal) { + int width = MATRIX.rows(), height = MATRIX.cols(); + unsigned char *pixels = (unsigned char *) calloc(3*width*height, sizeof (unsigned char)); + unsigned int err; + + for (int i = 0; i < (int) width; i++) for (int j = 0; j < (int) height; j++) { + int jj = height - j - 1; + double v = (MATRIX(i,j) - minVal)/(maxVal-minVal); + + if (v == 0.0) + pixels[3*(i+width*jj)+2] = 100; + else if (v == 1.0) + pixels[3*(i+width*jj)+0] = 100; + else if ((v > 0) && (v < 1.0)) { + pixels[3*(i+width*jj)+0] = MIN(MAX(0, 255.0*v), 255); + pixels[3*(i+width*jj)+1] = MIN(MAX(0, 255.0*v), 255); + pixels[3*(i+width*jj)+2] = MIN(MAX(0, 255.0*v), 255); + } else { + pixels[3*(i+width*jj)+0] = 255; + pixels[3*(i+width*jj)+1] = 255; + pixels[3*(i+width*jj)+2] = 255; + } + } + err = dispatch_loadbmp_encode_file(filename, pixels, width, height, LOADBMP_RGB); + if (err) std::cout << "FAILURE: " << EXIT_FAILURE << std::endl; + + free(pixels); +} +void outputMatrixI(string folder, const char *subfilename, MatrixXi& MATRIX, double minVal, double maxVal) { + string filename = folder + string(subfilename); + outputMatrixI(filename.c_str(), MATRIX, minVal, maxVal); +} + +void outputMatrixI(string folder, const char *subfolder, int i, MatrixXi& MATRIX, double minVal, double maxVal) { + string filename = folder + string(subfolder) + string("frame_") + std::to_string(i) + string(".bmp"); + outputMatrixI(filename.c_str(), MATRIX, minVal, maxVal); +} + + +void outputPoints(const char *filename, MatrixXd& POINTS, double normalizePoint, int RpxMax) { + int width = POINTS.rows(), height = POINTS.cols(); + unsigned char *pixels = (unsigned char *) calloc(3*width*height, sizeof (unsigned char)); + unsigned int err; + + // Then add the points + for (int i = 0; i < (int) width; i++) for (int j = 0; j < (int) height; j++) if (POINTS(i, j) != 0.0) { + double v = MIN(MAX(-1.0, POINTS(i,j)/normalizePoint), 1.0); + int Rpx = round(RpxMax*ABS(v)); + for (int ri = i-Rpx; ri <= i+Rpx; ri++) for (int rj = j-Rpx; rj <= j+Rpx; rj++) { + int rjj = height - rj - 1; + if ((ri > 0) && (ri < width-1) && (rj > 0) && (rj < height-1)) + if (SQ(ri-i) + SQ(rj-j) < SQ(Rpx)) { + pixels[3*(ri+width*rjj)+0] = (v > 0.0) ? 255 : 0; + pixels[3*(ri+width*rjj)+1] = 0; + pixels[3*(ri+width*rjj)+2] = (v < 0.0) ? 255 : 0; + } + } + } + + err = dispatch_loadbmp_encode_file(filename, pixels, width, height, LOADBMP_RGB); + if (err) std::cout << "FAILURE: " << EXIT_FAILURE << std::endl; + + free(pixels); +} +void outputPoints(string folder, const char *subfilename, MatrixXd& POINTS, double normalizePoint, int RpxMax) { + string filename = folder + string(subfilename); + outputPoints(filename.c_str(), POINTS, normalizePoint, RpxMax); +} +void outputPoints(string folder, const char *subfolder, int i, MatrixXd& POINTS, double normalizePoint, int RpxMax) { + string filename = folder + string(subfolder) + string("frame_") + std::to_string(i) + string(".bmp"); + outputPoints(filename.c_str(), POINTS, normalizePoint, RpxMax); +} + + +void outputPointsI(const char *filename, MatrixXi& POINTS, double normalizePoint, int RpxMax) { + int width = POINTS.rows(), height = POINTS.cols(); + unsigned char *pixels = (unsigned char *) calloc(3*width*height, sizeof (unsigned char)); + unsigned int err; + + // Then add the points + for (int i = 0; i < (int) width; i++) for (int j = 0; j < (int) height; j++) if (POINTS(i, j) != 0.0) { + double v = MIN(MAX(-1.0, POINTS(i,j)/normalizePoint), 1.0); + int Rpx = round(RpxMax*ABS(v)); + for (int ri = i-Rpx; ri <= i+Rpx; ri++) for (int rj = j-Rpx; rj <= j+Rpx; rj++) { + int rjj = height - rj - 1; + if ((ri > 0) && (ri < width-1) && (rj > 0) && (rj < height-1)) + if ((ri-i)*(ri-i) + (rj-j)*(rj-j) < Rpx*Rpx) { + pixels[3*(ri+width*rjj)+0] = (v > 0.0) ? 255 : 0; + pixels[3*(ri+width*rjj)+1] = 0; + pixels[3*(ri+width*rjj)+2] = (v < 0.0) ? 255 : 0; + } + } + } + err = dispatch_loadbmp_encode_file(filename, pixels, width, height, LOADBMP_RGB); + if (err) std::cout << "FAILURE: " << EXIT_FAILURE << std::endl; + + free(pixels); +} + + +void outputSites(const char *filename, Parameters& parameters, vector& Sites) { + MatrixXd SITES(parameters.N[0], parameters.N[1]); + for (vector::iterator it = Sites.begin(); it != Sites.end(); ++it) { + double xNorm = it->x / parameters.L[0], yNorm = it->y / parameters.L[1]; + int i = MIN(MAX(0, round(parameters.N[0]*xNorm)), parameters.N[0]-1), + j = MIN(MAX(0, round(parameters.N[1]*yNorm)), parameters.N[1]-1); + SITES(i,j) = (it->type == H2 ? -1 : 1) * (1/it->T); + } + outputPoints(filename, SITES, MAX(-SITES.minCoeff(), SITES.maxCoeff()), 14); +} +void outputSites(string folder, const char *subfilename, Parameters& parameters, vector& Sites) { + string filename = folder + string(subfilename); + outputSites(filename.c_str(), parameters, Sites); +} +void outputSites(string folder, const char *subfolder, int i, Parameters& parameters, vector& Sites) { + string filename = folder + string(subfolder) + string("frame_") + std::to_string(i) + string(".bmp"); + outputSites(filename.c_str(), parameters, Sites); +} + + + + + +/* + * Frame saving for visualization + */ + + +ofstream frameFile; +bool frameFileInit = false; +double outputDt, outputNext; +void initOutputFrame(const char *filename, Parameters& parameters, TimeSteps& timesteps, double tf, Topology& topology) { + frameFile = ofstream(filename); + if (!frameFile.is_open()) { + cerr << "Unable to create outputfile." << endl; + exit(1); + } + frameFile << parameters.L[0] << " " << parameters.L[1] << " " << eThick << " " << timesteps.dtElec << " "; + + int polyPos = 0; + for (int s = 0; s < 2; s++) { + frameFile << topology.polygonCount[s] << " "; + for (int i = 0; i <= topology.polygonCount[s]; i++) + frameFile << topology.polygonIndex[s][i] << " "; + int polygonPtCount = topology.polygonIndex[s][topology.polygonCount[s]]; + for (int i = 0; i < polygonPtCount; i++) { + frameFile << topology.polygonDefinitions[polyPos].x << ":" << topology.polygonDefinitions[polyPos].y << " "; + polyPos++; + } + } + frameFile << "!" << endl; + + if (outputMode == 1) { + outputDt = outputTimeFactor/outputFps; + outputNext = 0.0; + } + frameFileInit = true; +} +void outputFrame(string folder, const char *filename, double t, Parameters& parameters, TimeSteps& timesteps, Topology& topology, ListBubbles& Bubbles, ElectrolysisStatistics& Statistics) { + string path = folder + string(filename); + + // Check init + if (!frameFileInit) initOutputFrame(path.c_str(), parameters, timesteps, parameters.tf, topology); + + // Check if frame must be written + if (outputMode == 1) { + if (t < outputNext) return; + outputNext += outputDt; + } else if (outputMode == 2) { + if (t < outputStart) return; + if (t > outputStop) return; + } + + // Write global data + updateRemaining(Statistics, Bubbles); + frameFile << t << " " << Statistics.resistor << ":" << Statistics.overpotential << " "; + frameFile << Statistics.tension[0] << ":" << Statistics.tension[1] << ":" << Statistics.current[0] << ":" << Statistics.current[1] << " "; + frameFile << Statistics.contamination[1] << ":" << Statistics.collected[1] << ":" << Statistics.lost[1] << ":" << Statistics.remaining[1] << " "; + + // Write all bubble + list::iterator it; + for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) { + double itx, ity, itr; + double itz = eThick/2; + getAttachedBubbleLocation(itx, ity, itr, it); + frameFile << "0:" << itx << ":" << ity << ":" << itz << ":" << itr << ":" << it->quantity[0] << ":" << it->quantity[1] << " "; + } + for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) { + frameFile << "1:" << it->x << ":" << it->y << ":" << (eThick/2) << ":" << it->r << ":" << it->quantity[0] << ":" << it->quantity[1] << " "; + } + frameFile << "!" << endl; +} +void endOutputFrame(ofstream frameFile) { + frameFile.close(); +} diff --git a/src/Simulateur/utils/polygon.h b/src/Simulateur/utils/polygon.h new file mode 100644 index 0000000..5f2869f --- /dev/null +++ b/src/Simulateur/utils/polygon.h @@ -0,0 +1,110 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +// A C++ program to check if two given line segments intersect +// https://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/ + + +// Given three collinear points p, q, r, the function checks if +// point q lies on line segment 'pr' +bool onSegment(Point p, Point q, Point r) { + if (q.x <= MAX(p.x, r.x) && q.x >= MIN(p.x, r.x) && + q.y <= MAX(p.y, r.y) && q.y >= MIN(p.y, r.y)) + return true; + + return false; +} + +// To find orientation of ordered triplet (p, q, r). +// The function returns following values +// 0 --> p, q and r are collinear +// 1 --> Clockwise +// 2 --> Counterclockwise +int orientation(Point p, Point q, Point r) { + // See https://www.geeksforgeeks.org/orientation-3-ordered-points/ + // for details of below formula. + double val = (q.y - p.y) * (r.x - q.x) - + (q.x - p.x) * (r.y - q.y); + + if (val == 0) return 0; // collinear + + return (val > 0)? 1 : 2; // clock or counterclock wise +} + +// The main function that returns true if line segment 'p1q1' +// and 'p2q2' intersect. +bool doIntersect(Point p1, Point q1, Point p2, Point q2) { + // Find the four orientations needed for general and + // special cases + int o1 = orientation(p1, q1, p2); + int o2 = orientation(p1, q1, q2); + int o3 = orientation(p2, q2, p1); + int o4 = orientation(p2, q2, q1); + + // General case + if (o1 != o2 && o3 != o4) + return true; + + // Special Cases + // p1, q1 and p2 are collinear and p2 lies on segment p1q1 + if (o1 == 0 && onSegment(p1, p2, q1)) return true; + + // p1, q1 and q2 are collinear and q2 lies on segment p1q1 + if (o2 == 0 && onSegment(p1, q2, q1)) return true; + + // p2, q2 and p1 are collinear and p1 lies on segment p2q2 + if (o3 == 0 && onSegment(p2, p1, q2)) return true; + + // p2, q2 and q1 are collinear and q1 lies on segment p2q2 + if (o4 == 0 && onSegment(p2, q1, q2)) return true; + + return false; // Doesn't fall in any of the above cases +} + + + +// Returns true if the point p lies inside the polygon[] with n vertices +bool isInside(Point polygon[], int n, Point p) { + // There must be at least 3 vertices in polygon[] + if (n < 3) return false; + + // Create a point for line segment from p to infinite + Point extreme = {INF_NORMALIZED, p.y}; + + // Count intersections of the above line with sides of polygon + int count = 0, i = 0; + do + { + int next = (i+1)%n; + + // Check if the line segment from 'p' to 'extreme' intersects + // with the line segment from 'polygon[i]' to 'polygon[next]' + if (doIntersect(polygon[i], polygon[next], p, extreme)) + { + // If the point 'p' is collinear with line segment 'i-next', + // then check if it lies on segment. If it lies, return true, + // otherwise false + if (orientation(polygon[i], p, polygon[next]) == 0) + return onSegment(polygon[i], p, polygon[next]); + + count++; + } + i = next; + } while (i != 0); + + // Return true if count is odd, false otherwise + return count&1; // Same as (count%2 == 1) +} diff --git a/src/Simulateur/utils/progress.h b/src/Simulateur/utils/progress.h new file mode 100644 index 0000000..df7ee11 --- /dev/null +++ b/src/Simulateur/utils/progress.h @@ -0,0 +1,126 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +#define BAR_LENGTH_1 40 +#define BAR_LENGTH_2 5 + +int currentBar1 = -1, currentBar2a = -1, currentBar2b = 0; +char temps[4] = {'|', '/', '-', '\\'}; + +// Init progress bar +void progressInit(JNICaller* jniCaller) { +#ifdef UPDATE_PROGRESS + UPDATE_PROGRESS(*jniCaller, 10.0); +#endif + +#ifndef UPDATE_PROGRESS + string title = " Electrolysis simulation "; + int total = BAR_LENGTH_1 + 2, amount = round((total - title.length())/2); + cout << string(amount, '-') << title << string((total - amount - title.length()), '-') << endl; + + cout << "."; + for (int i = 0; i < BAR_LENGTH_1; i++) + cout << " "; + cout << "." << endl << flush; + + currentBar1 = 0; + currentBar2a = -1; + currentBar2b = 0; +#endif +} + +// Refresh progress bar +void progressPost(int bar1, int bar2a, int bar2b, JNICaller* jniCaller) { +#ifdef UPDATE_PROGRESS + if (jniCaller != nullptr) { + float percent = round(90*t/tf) + 10.0; + if (percent != *jniCaller.percent) { + UPDATE_PROGRESS(*jniCaller, percent); + *jniCaller.percent = percent; + } + } +#endif +#ifndef UPDATE_PROGRESS + // Usage of subprogress bar + if (useSubProgressBar) { + // Refresh console only if changes + if ((currentBar1 != bar1) || (currentBar2a != bar2a) || (currentBar2b != bar2b)) { + cout << '\r'; + // Bar + cout << "["; + for (int i = 0; i < BAR_LENGTH_1; i++) { + if (i < bar1) cout << "#"; + else cout << " "; + } + cout << "]"; + + // SubBar + if (bar2b != -1) { + cout << "("; + for (int i = 0; i < BAR_LENGTH_2; i++) { + if ((i < (bar2b-1)) || (bar2b == BAR_LENGTH_2)) cout << "-"; + else if (i == (bar2b-1)) cout << temps[bar2a]; + else cout << " "; + } + cout << ")"; + } else { + for (int i = 0; i < (2+BAR_LENGTH_2); i++) { + cout << " "; + } + } + cout << flush; + currentBar1 = bar1; + currentBar2a = bar2a; + currentBar2b = bar2b; + } + } else { + // Refresh console only if changes + if (currentBar1 != bar1) { + // Bar + if (currentBar1 == 0) cout << "["; + for (int i = 0; i < (bar1-currentBar1); i++) cout << "#"; + if (bar1 == BAR_LENGTH_1) cout << "]"; + cout << flush; + currentBar1 = bar1; + } + } +#endif +} +void progressRefresh1(double progress1, JNICaller* jniCaller) { + int bar1 = mathMinMax(round(progress1*BAR_LENGTH_1), 0, BAR_LENGTH_1); + progressPost(bar1, -1.0, -1.0, jniCaller); +} +void progressRefresh2b(double progress2b) { + int bar2b = round(progress2b*BAR_LENGTH_2); + progressPost(currentBar1, currentBar2a, (bar2b == 0) ? 1 : bar2b, nullptr); +} +void progressRefresh2a() { + int bar2a = (currentBar2a + 1) % 4; + if (currentBar2b == -1) progressRefresh2b(0.0); + progressPost(currentBar1, bar2a, currentBar2b, nullptr); +} + + + + + +// Error +void errorProgress(const std::string &message) { +#ifndef UPDATE_PROGRESS + cout << endl << flush; + cerr << "> Simulation stopped, caused by: " << endl << message << endl << flush; +#endif +} diff --git a/src/Simulateur/utils/typesdef.h b/src/Simulateur/utils/typesdef.h new file mode 100644 index 0000000..e692fbe --- /dev/null +++ b/src/Simulateur/utils/typesdef.h @@ -0,0 +1,123 @@ +// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters. +// +// Copyright (C) 2024 +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/ + +struct Point { + double x; + double y; +}; +struct Vect { + double x; + double y; +}; +struct Coord { + int i; + int j; +}; + +struct Topology { + int polygonCount[2]; + int* polygonIndex[2]; + Point* polygonDefinitions; +}; +struct Configuration { + double L[2]; + double J; + double tf; +}; +struct Performances { + double electrolysed[2], collected[2], contamination[2], lost[2], remaining[2]; + double resistor; +}; + +struct TimeSteps { + double dt, dtElec; + int i, iElec; +}; + +struct MatrixTension { + MatrixXd VALUES; + double min, max; + double phi_anode, phi_cathode; + double resistor; + double tNext; + bool sucess; +}; + +struct Site { + double x, y, T, growth; + double nx, ny; + int type; + double t, tmiss; + bool occuped; +}; + + +struct Wall { + Point P1, P2; + Vect u, n; + double length; +}; + +#define wallMaxmimumIndicators 2 +struct WallIndicator { + Point P; + Vect n; + double distance; + Wall* wall; +}; + +struct Bubble; +#define interactionMaxAmountAbsolute interactionLongMaximumRange*interactionLongMaximumRange +struct BubbleInteractor { + int interactionIndexs[interactionMaxAmountAbsolute]; + Bubble* bubbles[interactionMaxAmountAbsolute]; + double distancesSqNorm[interactionMaxAmountAbsolute]; + char imagedModes[interactionMaxAmountAbsolute]; + int count, countVirtual; + WallIndicator walls[wallMaxmimumIndicators]; +}; +struct Bubble { + double x, y, r; + double ux, uy; + int statLife, statCollision; + bool alive; + Site* site; + BubbleInteractor interactor; + double quantity[2]; + double tClogg; +}; + + +struct ListBubbles { + list attached; + list detached; + bool flagAttached, flagDetached; + int countAttached, countDetached; + double maxRadius, maxVelocity; + unsigned long tNextNucleation; +}; +struct ElectrolysisStatistics { + double collected[2], lost[2], contamination[2], remaining[2]; + double resistor, overpotential; + double tension[2], current[2]; +}; + +struct ElectrodeNode { + int id; + int nx_anode, ny_anode; + int nx_cathode, ny_cathode; + double denomInv; +};