From 3fd2c17c73584cc0ec1d9fcee94d2a3bc3196390 Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Mon, 29 Jun 2026 20:11:52 +0100 Subject: [PATCH 1/4] Construct electrode reaction current only once for consistency --- .../ElectrodeReactionCurrentCoefficient.hpp | 24 ++++++++--------- operators/EChemOperator.cpp | 26 ++++++++++++++----- operators/EChemOperator.hpp | 2 ++ 3 files changed, 33 insertions(+), 19 deletions(-) diff --git a/coefficients/ElectrodeReactionCurrentCoefficient.hpp b/coefficients/ElectrodeReactionCurrentCoefficient.hpp index 832331a..1c04abe 100644 --- a/coefficients/ElectrodeReactionCurrentCoefficient.hpp +++ b/coefficients/ElectrodeReactionCurrentCoefficient.hpp @@ -7,26 +7,26 @@ class ElectrodeReactionCurrentCoefficient : public Coefficient ExchangeCurrentCoefficient * _jex = nullptr; OverPotentialCoefficient * _op = nullptr; - const real_t _a; - const Region _r; - const int _sign; + real_t _a; + Region _r; + int _sign; public: /// P2D - ElectrodeReactionCurrentCoefficient(const Region & r, - const int & sign, - ExchangeCurrentCoefficient & jex, + ElectrodeReactionCurrentCoefficient(ExchangeCurrentCoefficient & jex, OverPotentialCoefficient & op) : _jex(&jex), - _op(&op), - _a(r == NE ? AN * LNE / NNE * NX - : r == PE ? AP * LPE / NPE * NX - : 0), - _r(r), - _sign(sign) + _op(&op) { } + void SetRegionSign(const Region &r, const int & sign) + { + _a = r == NE ? AN * LNE / NNE * NX : r == PE ? AP * LPE / NPE * NX : 0; + _r = r; + _sign = sign; + } + /// P2D virtual real_t Eval(ElementTransformation & T, const IntegrationPoint & ip) override { diff --git a/operators/EChemOperator.cpp b/operators/EChemOperator.cpp index 77e9a4c..ef967b4 100644 --- a/operators/EChemOperator.cpp +++ b/operators/EChemOperator.cpp @@ -84,6 +84,12 @@ EChemOperator::EChemOperator(ParFiniteElementSpace *& x_h1space, _ec_gf = CE0; _sc_gf = 0.; + ConstructExchangeCurrent(); + ConstructOpenCircuitPotential(); + ConstructOverPotential(); + ConstructReactionCurrent(); + ConstructElectrodeReactionCurrent(); + // Construct equation ojects, first the 3 macro equations, then the NPAR micro eqs if (P2D) { @@ -117,11 +123,6 @@ EChemOperator::EChemOperator(ParFiniteElementSpace *& x_h1space, } } - ConstructExchangeCurrent(); - ConstructOpenCircuitPotential(); - ConstructOverPotential(); - ConstructReactionCurrent(); - if (P2D) { // Construct space for discontinuous functions like the reaction current j @@ -318,9 +319,9 @@ EChemOperator::GetElectrodeReactionCurrent(const Region & r, const int & sign) "Cannot get partial electrode reaction current, only negative (NE) and positive " "electrodes (PE) are supported."); - ElectrodeReactionCurrentCoefficient j(r, sign, *_jex, *_op); + _je->SetRegionSign(r, sign); QuadratureSpace x_qspace(_x_h1space->GetParMesh(), 2 * _x_h1space->FEColl()->GetOrder()); - return x_qspace.Integrate(j); + return x_qspace.Integrate(*_je); } // @@ -399,6 +400,17 @@ EChemOperator::ConstructReactionCurrent() _j = new ReactionCurrentCoefficient(T, *_jex, *_op); } +// +// Electrode Reaction Current +// + +void +EChemOperator::ConstructElectrodeReactionCurrent() +{ + if (P2D) + _je = new ElectrodeReactionCurrentCoefficient(*_jex, *_op); +} + // // Exchange Current // diff --git a/operators/EChemOperator.hpp b/operators/EChemOperator.hpp index c059415..803d8ba 100644 --- a/operators/EChemOperator.hpp +++ b/operators/EChemOperator.hpp @@ -45,6 +45,7 @@ class EChemOperator : public TimeDependentOperator /// Coefficients for derived, i.e. not solved for, quantities ReactionCurrentCoefficient * _j; + ElectrodeReactionCurrentCoefficient * _je; ExchangeCurrentCoefficient * _jex; OpenCircuitPotentialCoefficient * _ocp; OverPotentialCoefficient * _op; @@ -156,6 +157,7 @@ class EChemOperator : public TimeDependentOperator delete _sc[p]; delete _j; + delete _je; delete _jex; delete _ocp; delete _op; From 541f8db3f1eaf002954d3d2bc08c5898c9d74c5c Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Mon, 29 Jun 2026 20:59:22 +0100 Subject: [PATCH 2/4] Unify coefficient construction --- operators/EChemOperator.cpp | 91 ++++++++++++++----------------------- operators/EChemOperator.hpp | 5 +- 2 files changed, 34 insertions(+), 62 deletions(-) diff --git a/operators/EChemOperator.cpp b/operators/EChemOperator.cpp index ef967b4..aa05f3d 100644 --- a/operators/EChemOperator.cpp +++ b/operators/EChemOperator.cpp @@ -84,11 +84,7 @@ EChemOperator::EChemOperator(ParFiniteElementSpace *& x_h1space, _ec_gf = CE0; _sc_gf = 0.; - ConstructExchangeCurrent(); - ConstructOpenCircuitPotential(); - ConstructOverPotential(); - ConstructReactionCurrent(); - ConstructElectrodeReactionCurrent(); + ConstructCoefficients(); // Construct equation ojects, first the 3 macro equations, then the NPAR micro eqs if (P2D) @@ -376,6 +372,38 @@ EChemOperator::GetParticleReactionCurrent() return j; } +void +EChemOperator::ConstructCoefficients() +{ + if (SPM) + { + _jex = new ExchangeCurrentCoefficient( + KN, KP, GetSurfaceConcentration(NE), GetSurfaceConcentration(PE), CE0); + _ocp = new OpenCircuitPotentialCoefficient( + UN, UP, GetSurfaceConcentration(NE), GetSurfaceConcentration(PE)); + _op = new OverPotentialCoefficient(T, *_jex); + _j = new ReactionCurrentCoefficient(); + } + else if (SPMe) + { + _jex = new ExchangeCurrentCoefficient( + KN, KP, GetSurfaceConcentration(NE), GetSurfaceConcentration(PE), _ec_gfc); + _ocp = new OpenCircuitPotentialCoefficient( + UN, UP, GetSurfaceConcentration(NE), GetSurfaceConcentration(PE)); + _op = new OverPotentialCoefficient(T, *_jex); + _j = new ReactionCurrentCoefficient(); + } + else if (P2D) + { + _jex = new ExchangeCurrentCoefficient(KN, KP, _sc_gfc, _ec_gfc); + _ocp = new OpenCircuitPotentialCoefficient(UN, UP, _sc_gfc); + _op = new OverPotentialCoefficient( + GetReferencePotential(E), GetReferencePotential(PE), _sp_gfc, _ep_gfc, *_ocp); + _j = new ReactionCurrentCoefficient(T, *_jex, *_op); + _je = new ElectrodeReactionCurrentCoefficient(*_jex, *_op); + } +} + // // Reaction Current // @@ -391,26 +419,6 @@ EChemOperator::GetReactionCurrent(const Region & r) return _j->Eval()(r); } -void -EChemOperator::ConstructReactionCurrent() -{ - if (SPM || SPMe) - _j = new ReactionCurrentCoefficient(); - else if (P2D) - _j = new ReactionCurrentCoefficient(T, *_jex, *_op); -} - -// -// Electrode Reaction Current -// - -void -EChemOperator::ConstructElectrodeReactionCurrent() -{ - if (P2D) - _je = new ElectrodeReactionCurrentCoefficient(*_jex, *_op); -} - // // Exchange Current // @@ -426,19 +434,6 @@ EChemOperator::GetExchangeCurrent(const Region & r) return _jex->Eval()(r); } -void -EChemOperator::ConstructExchangeCurrent() -{ - if (SPM) - _jex = new ExchangeCurrentCoefficient( - KN, KP, GetSurfaceConcentration(NE), GetSurfaceConcentration(PE), CE0); - else if (SPMe) - _jex = new ExchangeCurrentCoefficient( - KN, KP, GetSurfaceConcentration(NE), GetSurfaceConcentration(PE), _ec_gfc); - else if (P2D) - _jex = new ExchangeCurrentCoefficient(KN, KP, _sc_gfc, _ec_gfc); -} - // // Open Circuit Potential // @@ -454,16 +449,6 @@ EChemOperator::GetOpenCircuitPotential(const Region & r) return _ocp->Eval()(r); } -void -EChemOperator::ConstructOpenCircuitPotential() -{ - if (SPM || SPMe) - _ocp = new OpenCircuitPotentialCoefficient( - UN, UP, GetSurfaceConcentration(NE), GetSurfaceConcentration(PE)); - else if (P2D) - _ocp = new OpenCircuitPotentialCoefficient(UN, UP, _sc_gfc); -} - // // OverPotential // @@ -478,16 +463,6 @@ EChemOperator::GetOverPotential(const Region & r) return _op->Eval()(r); } -void -EChemOperator::ConstructOverPotential() -{ - if (SPM || SPMe) - _op = new OverPotentialCoefficient(T, *_jex); - else if (P2D) - _op = new OverPotentialCoefficient( - GetReferencePotential(E), GetReferencePotential(PE), _sp_gfc, _ep_gfc, *_ocp); -} - // // Voltage // diff --git a/operators/EChemOperator.hpp b/operators/EChemOperator.hpp index 803d8ba..8b2a2e8 100644 --- a/operators/EChemOperator.hpp +++ b/operators/EChemOperator.hpp @@ -120,10 +120,7 @@ class EChemOperator : public TimeDependentOperator Array GetParticleReactionCurrent(); /// Construct coefficients for derived quantities - void ConstructReactionCurrent(); - void ConstructExchangeCurrent(); - void ConstructOpenCircuitPotential(); - void ConstructOverPotential(); + void ConstructCoefficients(); /// Helpers for quantities which are constant within a region const real_t & GetSurfaceConcentration(const Region & r); From bf8dd3a03a9a5eb02336a70ea2ad78fc48078409 Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Mon, 29 Jun 2026 21:40:42 +0100 Subject: [PATCH 3/4] Move Update() input parameters to eqs. constructors --- equations/ElectrolyteConcentration.cpp | 6 +++--- equations/ElectrolyteConcentration.hpp | 8 ++++++-- equations/ElectrolytePotential.cpp | 8 ++++---- equations/ElectrolytePotential.hpp | 6 ++++-- equations/Equation.hpp | 3 --- equations/SolidConcentration.cpp | 4 ++-- equations/SolidConcentration.hpp | 2 +- equations/SolidPotential.cpp | 4 ++-- equations/SolidPotential.hpp | 7 +++++-- operators/EChemOperator.cpp | 14 +++++++------- 10 files changed, 34 insertions(+), 28 deletions(-) diff --git a/equations/ElectrolyteConcentration.cpp b/equations/ElectrolyteConcentration.cpp index c6e27c0..92eb226 100644 --- a/equations/ElectrolyteConcentration.cpp +++ b/equations/ElectrolyteConcentration.cpp @@ -2,7 +2,7 @@ #include "equations/ElectrolyteConcentration.hpp" void -ElectrolyteConcentration::Update(const GridFunctionCoefficient & ec_gfc, const Coefficient & j) +ElectrolyteConcentration::Update() { // Mass coefficient. Vector mass_vec({/* NE */ EPS_N /* length scaling */ * (LNE / NNE * NX), @@ -17,14 +17,14 @@ ElectrolyteConcentration::Update(const GridFunctionCoefficient & ec_gfc, const C /* PE */ (1 - TPLUS) * AP /* length scaling */ * (LPE / NPE * NX)}); PWConstCoefficient source_part(source_vec); - ProductCoefficient source(source_part, const_cast(j)); + ProductCoefficient source(source_part, j); // Diffusion coefficient. Vector D_scale_vec({/* NE */ BNE /* length scaling */ / (LNE / NNE * NX), /* SEP */ BSEP /* length scaling */ / (LSEP / NSEP * NX), /* PE */ BPE /* length scaling */ / (LPE / NPE * NX)}); - TransformedCoefficient D_coeff(&const_cast(ec_gfc), DE); + TransformedCoefficient D_coeff(&ec_gfc, DE); PWConstCoefficient D_scale_coeff(D_scale_vec); ProductCoefficient D(D_scale_coeff, D_coeff); diff --git a/equations/ElectrolyteConcentration.hpp b/equations/ElectrolyteConcentration.hpp index 26a7f51..64cdf5f 100644 --- a/equations/ElectrolyteConcentration.hpp +++ b/equations/ElectrolyteConcentration.hpp @@ -7,7 +7,11 @@ using namespace mfem; class ElectrolyteConcentration : public Equation { +protected: + Coefficient & j; + GridFunctionCoefficient & ec_gfc; + public: - using Equation::Equation; - virtual void Update(const GridFunctionCoefficient & ec_gfc, const Coefficient & j) override; + ElectrolyteConcentration(ParFiniteElementSpace & f, Coefficient & j, GridFunctionCoefficient & ec_gfc) : Equation(f), j(j), ec_gfc(ec_gfc) {}; + virtual void Update(); }; diff --git a/equations/ElectrolytePotential.cpp b/equations/ElectrolytePotential.cpp index 4515e79..a1b8c40 100644 --- a/equations/ElectrolytePotential.cpp +++ b/equations/ElectrolytePotential.cpp @@ -2,7 +2,7 @@ #include "equations/ElectrolytePotential.hpp" void -ElectrolytePotential::Update(const GridFunctionCoefficient & ec_gfc, const Coefficient & j) +ElectrolytePotential::Update() { // Source term. Vector source_vec({/* NE */ AN /* length scaling */ * (LNE / NNE * NX), @@ -10,7 +10,7 @@ ElectrolytePotential::Update(const GridFunctionCoefficient & ec_gfc, const Coeff /* PE */ AP /* length scaling */ * (LPE / NPE * NX)}); PWConstCoefficient source_part(source_vec); - ProductCoefficient source(source_part, const_cast(j)); + ProductCoefficient source(source_part, j); Vector b_vec({/* NE */ BNE /* length scaling */ / (LNE / NNE * NX), /* SEP */ BSEP /* length scaling */ / (LSEP / NSEP * NX), @@ -18,11 +18,11 @@ ElectrolytePotential::Update(const GridFunctionCoefficient & ec_gfc, const Coeff PWConstCoefficient b_part(b_vec); - TransformedCoefficient kappa(&const_cast(ec_gfc), Kappa); + TransformedCoefficient kappa(&ec_gfc, Kappa); ProductCoefficient kappa_eff(b_part, kappa); GradientGridFunctionCoefficient grad_ec(ec_gfc.GetGridFunction()); - RatioCoefficient ec_inv(2 * T * (1 - TPLUS), const_cast(ec_gfc)); + RatioCoefficient ec_inv(2 * T * (1 - TPLUS), ec_gfc); ScalarVectorProductCoefficient grad_ln_ec(ec_inv, grad_ec); ScalarVectorProductCoefficient grad_ln_ec_kappad(kappa_eff, grad_ln_ec); diff --git a/equations/ElectrolytePotential.hpp b/equations/ElectrolytePotential.hpp index 257c9e8..0da3673 100644 --- a/equations/ElectrolytePotential.hpp +++ b/equations/ElectrolytePotential.hpp @@ -8,15 +8,17 @@ using namespace mfem; class ElectrolytePotential : public Equation { protected: + Coefficient & j; + GridFunctionCoefficient & ec_gfc; ParLinearForm * Qc = nullptr; Vector bc; public: - ElectrolytePotential(ParFiniteElementSpace & f) : Equation(f), bc(f.GetTrueVSize()) + ElectrolytePotential(ParFiniteElementSpace & f, Coefficient & j, GridFunctionCoefficient & ec_gfc) : Equation(f), j(j), ec_gfc(ec_gfc), bc(f.GetTrueVSize()) { f.GetEssentialTrueDofs(Array({1, 0}), ess_tdof_list); } - virtual void Update(const GridFunctionCoefficient & ec_gfc, const Coefficient & j) override; + virtual void Update(); void Reset() { delete K; diff --git a/equations/Equation.hpp b/equations/Equation.hpp index 0ebcd62..aa0bb99 100644 --- a/equations/Equation.hpp +++ b/equations/Equation.hpp @@ -31,9 +31,6 @@ class Equation const HypreParMatrix & GetK() const { return Kmat; }; const Vector & GetZ() const { return b; }; - virtual void Update(const Coefficient & j) {} - virtual void Update(const GridFunctionCoefficient & u, const Coefficient & j) {} - virtual ~Equation() { delete M; diff --git a/equations/SolidConcentration.cpp b/equations/SolidConcentration.cpp index 901a1aa..40364b7 100644 --- a/equations/SolidConcentration.cpp +++ b/equations/SolidConcentration.cpp @@ -1,7 +1,7 @@ #include "equations/SolidConcentration.hpp" void -SolidConcentration::Update(const Coefficient & j) +SolidConcentration::Update(const real_t & j) { MFEM_ASSERT(particle_region == NE || particle_region == PE, "Particle not in electrode!"); @@ -11,7 +11,7 @@ SolidConcentration::Update(const Coefficient & j) FunctionCoefficient r2([](const Vector & r) { return r(0) * r(0); }); ProductCoefficient dr2(D / R / R, r2); - ProductCoefficient jjr2(const_cast(j), r2); + ProductCoefficient jjr2(j, r2); ProductCoefficient jr2(-1. / R / t_scale, jjr2); if (!M) diff --git a/equations/SolidConcentration.hpp b/equations/SolidConcentration.hpp index 57fa198..bd98cf9 100644 --- a/equations/SolidConcentration.hpp +++ b/equations/SolidConcentration.hpp @@ -37,7 +37,7 @@ class SolidConcentration : public Equation { } - virtual void Update(const Coefficient & j) override; + virtual void Update(const real_t & j); real_t SurfaceConcentration(const BlockVector & x); int GetParticleRank() { return particle_rank; } int GetSurfaceRank() { return surface_rank; } diff --git a/equations/SolidPotential.cpp b/equations/SolidPotential.cpp index f5dd227..323af9b 100644 --- a/equations/SolidPotential.cpp +++ b/equations/SolidPotential.cpp @@ -2,7 +2,7 @@ #include "equations/SolidPotential.hpp" void -SolidPotential::Update(const Coefficient & j) +SolidPotential::Update() { // Source term. Vector source_vec({/* NE */ -AN /* length scaling */ * (LNE / NNE * NX), @@ -10,7 +10,7 @@ SolidPotential::Update(const Coefficient & j) /* PE */ -AP /* length scaling */ * (LPE / NPE * NX)}); PWConstCoefficient source_part(source_vec); - ProductCoefficient source(source_part, const_cast(j)); + ProductCoefficient source(source_part, j); // Effective conductivity (does not account for electrode filler). Vector sigma_vec({/* NE */ (1 - EPS_N) * SIGN /* length scaling */ / (LNE / NNE * NX), diff --git a/equations/SolidPotential.hpp b/equations/SolidPotential.hpp index 860fb0d..f3bf6da 100644 --- a/equations/SolidPotential.hpp +++ b/equations/SolidPotential.hpp @@ -7,10 +7,13 @@ using namespace mfem; class SolidPotential : public Equation { +protected: + Coefficient &j; + public: - SolidPotential(ParFiniteElementSpace & f) : Equation(f) + SolidPotential(ParFiniteElementSpace & f, Coefficient & j) : Equation(f), j(j) { f.GetEssentialTrueDofs(Array({1, 1}), ess_tdof_list); } - virtual void Update(const Coefficient & j) override; + virtual void Update(); }; diff --git a/operators/EChemOperator.cpp b/operators/EChemOperator.cpp index aa05f3d..6e431c1 100644 --- a/operators/EChemOperator.cpp +++ b/operators/EChemOperator.cpp @@ -89,10 +89,10 @@ EChemOperator::EChemOperator(ParFiniteElementSpace *& x_h1space, // Construct equation ojects, first the 3 macro equations, then the NPAR micro eqs if (P2D) { - _ep = new ElectrolytePotential(*_x_h1space); - _sp = new SolidPotential(*_x_h1space); + _ep = new ElectrolytePotential(*_x_h1space, *_j, _ec_gfc); + _sp = new SolidPotential(*_x_h1space, *_j); } - _ec = new ElectrolyteConcentration(*_x_h1space); + _ec = new ElectrolyteConcentration(*_x_h1space, *_j, _ec_gfc); if (SPM || SPMe) { @@ -217,17 +217,17 @@ EChemOperator::SetGridFunctionsFromTrueVectors() void EChemOperator::UpdatePotentialEquations() { - _ep->Update(_ec_gfc, *_j); - _sp->Update(*_j); + _ep->Update(); + _sp->Update(); } void EChemOperator::UpdateConcentrationEquations() { - _ec->Update(_ec_gfc, *_j); + _ec->Update(); const Array & j = GetParticleReactionCurrent(); for (unsigned p = 0; p < NPAR; p++) - _sc[p]->Update(ConstantCoefficient(j[p])); + _sc[p]->Update(j[p]); } // From 70f0bcdf138e8bf279013ea8fc585df974381985 Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Mon, 29 Jun 2026 21:43:25 +0100 Subject: [PATCH 4/4] Add grad-log coefficient for ElectrolytePotential --- coefficients/EPGradientLogCoefficient.hpp | 37 +++++++++++++++++++++++ equations/ElectrolytePotential.cpp | 7 +---- equations/ElectrolytePotential.hpp | 4 ++- makefile | 3 +- 4 files changed, 43 insertions(+), 8 deletions(-) create mode 100644 coefficients/EPGradientLogCoefficient.hpp diff --git a/coefficients/EPGradientLogCoefficient.hpp b/coefficients/EPGradientLogCoefficient.hpp new file mode 100644 index 0000000..f5ecca4 --- /dev/null +++ b/coefficients/EPGradientLogCoefficient.hpp @@ -0,0 +1,37 @@ +#include "mfem.hpp" +using namespace mfem; + +class EPGradientLogCoefficient : public VectorCoefficient +{ +private: + GridFunctionCoefficient * _electrolyte_concentration_gfc = nullptr; + GradientGridFunctionCoefficient _electrolyte_concentration_ggfc; + +public: + EPGradientLogCoefficient(GridFunctionCoefficient & ec) + : VectorCoefficient(1), + _electrolyte_concentration_gfc(&ec), + _electrolyte_concentration_ggfc(ec.GetGridFunction()) + { + } + + virtual void Eval(Vector &V, ElementTransformation & T, const IntegrationPoint & ip) override + { + _electrolyte_concentration_ggfc.Eval(V, T, ip); + const real_t ec = _electrolyte_concentration_gfc->Eval(T, ip); + V[0] *= 2 * constants::T * (1 - TPLUS) * Kappa(ec) / ec; + + switch (T.Attribute) + { + case NE: + V[0] *= BNE / (LNE / NNE * NX); + break; + case SEP: + V[0] *= BSEP / (LSEP / NSEP * NX); + break; + case PE: + V[0] *= BPE / (LPE / NPE * NX); + break; + } + } +}; diff --git a/equations/ElectrolytePotential.cpp b/equations/ElectrolytePotential.cpp index a1b8c40..8233a81 100644 --- a/equations/ElectrolytePotential.cpp +++ b/equations/ElectrolytePotential.cpp @@ -21,11 +21,6 @@ ElectrolytePotential::Update() TransformedCoefficient kappa(&ec_gfc, Kappa); ProductCoefficient kappa_eff(b_part, kappa); - GradientGridFunctionCoefficient grad_ec(ec_gfc.GetGridFunction()); - RatioCoefficient ec_inv(2 * T * (1 - TPLUS), ec_gfc); - ScalarVectorProductCoefficient grad_ln_ec(ec_inv, grad_ec); - ScalarVectorProductCoefficient grad_ln_ec_kappad(kappa_eff, grad_ln_ec); - if (!K) { K = new ParBilinearForm(&fespace); @@ -37,7 +32,7 @@ ElectrolytePotential::Update() if (!Qc) { Qc = new ParLinearForm(&fespace); - Qc->AddDomainIntegrator(new DomainLFGradIntegrator(grad_ln_ec_kappad)); + Qc->AddDomainIntegrator(new DomainLFGradIntegrator(gl)); Qc->Assemble(); Qc->ParallelAssemble(bc); } diff --git a/equations/ElectrolytePotential.hpp b/equations/ElectrolytePotential.hpp index 0da3673..92419d6 100644 --- a/equations/ElectrolytePotential.hpp +++ b/equations/ElectrolytePotential.hpp @@ -2,6 +2,7 @@ #include "mfem.hpp" #include "equations/Equation.hpp" +#include "coefficients/EPGradientLogCoefficient.hpp" using namespace mfem; @@ -10,11 +11,12 @@ class ElectrolytePotential : public Equation protected: Coefficient & j; GridFunctionCoefficient & ec_gfc; + EPGradientLogCoefficient gl; ParLinearForm * Qc = nullptr; Vector bc; public: - ElectrolytePotential(ParFiniteElementSpace & f, Coefficient & j, GridFunctionCoefficient & ec_gfc) : Equation(f), j(j), ec_gfc(ec_gfc), bc(f.GetTrueVSize()) + ElectrolytePotential(ParFiniteElementSpace & f, Coefficient & j, GridFunctionCoefficient & ec_gfc) : Equation(f), j(j), ec_gfc(ec_gfc), gl(ec_gfc), bc(f.GetTrueVSize()) { f.GetEssentialTrueDofs(Array({1, 0}), ess_tdof_list); } diff --git a/makefile b/makefile index 04037df..fc9f595 100644 --- a/makefile +++ b/makefile @@ -24,7 +24,8 @@ CF_INC_FILES = ExchangeCurrentCoefficient.hpp \ ReactionCurrentCoefficient.hpp \ ElectrodeReactionCurrentCoefficient.hpp \ OpenCircuitPotentialCoefficient.hpp \ - OverPotentialCoefficient.hpp + OverPotentialCoefficient.hpp \ + EPGradientLogCoefficient.hpp CL_INC_FILES = LGM50.hpp PR_SRC_FILES = settings.cpp PR_INC_FILES = settings.hpp \