Page MenuHomec4science
Files F86263128

transport_program.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Oct 5, 09:23

transport_program.cpp
View Options

/* =============================================================================
Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
Copyright (c) 2017 F. Georget <fabien.georget@epfl.ch> EPFL
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
3. Neither the name of the copyright holder nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
============================================================================= */
#include "transport_program.hpp"
#include "variables.hpp"
#include "dfpm/meshes/mesh1d.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp_common/micpsolver/micpsolver_structs.hpp"
#include "specmicp/adimensional/adimensional_system_structs.hpp"
#include "reactmicp/solver/staggers_base/chemistry_stagger_base.hpp"
#include "specmicp_database/database_holder.hpp"
#include "specmicp_common/physics/maths.hpp"
#include "unsupported/Eigen/SparseExtra"
#include <iostream>
#include "specmicp/adimensional/adimensional_system_solver.hpp"
namespace specmicp {
namespace reactmicp {
namespace systems {
namespace chloride {
using specmicp_solver_f = std::function<ReturnCode (index_t, AdimensionalSystemSolution&)>;
ChlorideProgram::ChlorideProgram(ChlorideVariablesPtr vars,
std::vector<index_t> list_fixed_dof):
database::DatabaseHolder(vars->get_database()),
m_vars(vars),
m_mesh(vars->get_mesh())
{
m_scaling.setOnes(vars->ndof());
number_equations(list_fixed_dof);
}
void ChlorideProgram::number_equations(std::vector<index_t>& list_fixed_dof)
{
const auto tot_dof = m_vars->tot_dof();
m_ideq.resize(tot_dof);
m_ideq.setZero();
// boundary conditions is provided as a set of fixed dofs
//
// Only zero-flux bcs possible for now
for (auto ind: list_fixed_dof)
{
m_ideq(ind) = no_equation;
}
// standard neq counting
index_t neq = 0;
for (auto dof: RangeIterator<index_t>(tot_dof))
{
if (m_ideq(dof) == no_equation) continue;
m_ideq(dof) = neq;
++neq;
}
m_neq=neq;
}
void ChlorideProgram::compute_residuals(
const Vector& displacement,
const Vector& velocity,
Vector& residual,
bool use_chemistry_rates,
bool use_grad_psi
)
{
const index_t ndof = m_vars->ndof();
residual = Vector::Zero(m_neq);
m_current.setZero(m_mesh->nb_elements());
#ifdef SPECMICP_HAVE_OPENMP
#pragma omp parallel for schedule(dynamic, 5)
#endif
for (index_t e=0; e<m_mesh->nb_elements(); ++e)
{
Vector eresidual = Vector::Zero(2*m_vars->ndof());
const auto node_0 = m_mesh->get_node(e, 0);
const auto node_1 = m_mesh->get_node(e, 1);
auto evel = m_vars->get_element_dofs(e, velocity);
auto edisp = m_vars->get_element_dofs(e, displacement);
compute_element_residuals(e, edisp, evel,
m_vars->adim_solution(node_0),
m_vars->adim_solution(node_1),
eresidual, use_chemistry_rates,
use_grad_psi);
for (auto enode: RangeIterator<index_t>(2) )
{
index_t node = m_mesh->get_node(e, enode);
for (auto nodal_dof: RangeIterator<index_t>(ndof))
{
const auto ideq = m_ideq(m_vars->ndof_to_dof(node, nodal_dof));
if (ideq != no_equation) {
#ifdef SPECMICP_HAVE_OPENMP
#pragma omp atomic
#endif
residual(ideq) += eresidual(m_vars->ndof_to_edof(enode, nodal_dof));
}
}
}
}
}
void ChlorideProgram::compute_element_residuals(
index_t element,
const Vector& edisplacement,
const Vector& evelocity,
const AdimensionalSystemSolution& adim_sol_0,
const AdimensionalSystemSolution& adim_sol_1,
Vector& eresidual,
bool use_chemistry_rate,
bool use_grad_psi
)
{
// general variables that will be used
const index_t node_0 = m_mesh->get_node(element, 0);
const index_t node_1 = m_mesh->get_node(element, 1);
const scalar_t mass_coeff_0 = m_mesh->get_volume_cell_element(element, 0);
const scalar_t mass_coeff_1 = m_mesh->get_volume_cell_element(element, 1);
const scalar_t mu_face = average<Average::harmonic>(m_vars->resistance_factor(node_0), m_vars->resistance_factor(node_1));
const scalar_t flux_geom_coeff = -(m_mesh->get_face_area(element) / m_mesh->get_dx(element));
scalar_t courant_face = 0;
AdimensionalSystemSolutionExtractor extr_0(adim_sol_0, m_vars->get_database(), m_vars->get_units());
AdimensionalSystemSolutionExtractor extr_1(adim_sol_1, m_vars->get_database(), m_vars->get_units());
// aqueous components
for (auto i: m_data->range_aqueous_component())
{
const auto dof_0 = m_vars->dof_component(node_0, i);
const auto dof_1 = m_vars->dof_component(node_1, i);
const auto edof_0 = m_vars->edof_component(0, i);
const auto edof_1 = m_vars->edof_component(1, i);
const scalar_t velocity_0 = mass_coeff_0*(evelocity(edof_0)*m_vars->porosity(node_0)
+m_vars->velocity_porosity(node_0)*edisplacement(edof_0));
const scalar_t velocity_1 = mass_coeff_1*(evelocity(edof_1)*m_vars->porosity(node_1)
+ m_vars->velocity_porosity(node_1)*edisplacement(edof_1));
scalar_t transport_0 = 0.0;
scalar_t transport_1 = 0.0;
// advection
if (m_vars->fluid_velocity() != 0.0)
{
const scalar_t vel = m_vars->fluid_velocity();
scalar_t flux_advection = m_mesh->get_face_area(element)*vel;
if (vel > 0)
{
flux_advection *= (edisplacement(edof_0) - edisplacement(edof_1));
transport_1 += flux_advection;
}
else
{
flux_advection *= (edisplacement(edof_1) - edisplacement(edof_0));
transport_0 -= flux_advection;
}
}
// diffusion and nernst planck
{
scalar_t diff_0 = 0;
scalar_t diff_1 = 0;
scalar_t np_plus_0 = 0;
scalar_t np_minus_0 = 0;
scalar_t np_plus_1 = 0;
scalar_t np_minus_1 = 0;
const scalar_t mola_0 = extr_0.molality_component(i);
const scalar_t mola_1 = extr_1.molality_component(i);
const scalar_t gamma_0 = extr_0.activity_coefficient_component(i);
const scalar_t gamma_1 = extr_1.activity_coefficient_component(i);
const scalar_t mola_face = average<Average::arithmetic>(mola_0, mola_1);
const auto diff_c = m_vars->intrinsic_diff_coeff_component(i)*mu_face*flux_geom_coeff*(mola_0 - mola_1);
const auto diff_g = m_vars->intrinsic_diff_coeff_component(i)*mu_face*flux_geom_coeff*mola_face*(gamma_0 - gamma_1);
diff_0 = + diff_c + diff_g;
diff_1 = - diff_c - diff_g;
const scalar_t charge = m_data->charge_component(i);
if (charge != 0.0) {
const scalar_t coeff_np = m_vars->intrinsic_diff_coeff_component(i)*mu_face*charge;
if (charge > 0) {
np_plus_0 += coeff_np*mola_0;
np_plus_1 += coeff_np*mola_1;
} else {
np_minus_0 += coeff_np*mola_0;
np_minus_1 += coeff_np*mola_1;
}
}
for (auto j: m_data->range_aqueous())
{
const auto nu_ji = m_data->nu_aqueous(j, i);
if (nu_ji == 0.0) continue;
const scalar_t mola_0 = extr_0.molality_aqueous(j);
const scalar_t mola_1 = extr_1.molality_aqueous(j);
const scalar_t gamma_0 = extr_0.activity_coefficient_aqueous(j);
const scalar_t gamma_1 = extr_1.activity_coefficient_aqueous(j);
const scalar_t mola_face = average<Average::arithmetic>(mola_0, mola_1);
const auto diffj_c = m_vars->intrinsic_diff_coeff_aqueous(j)*mu_face*flux_geom_coeff*(mola_0 - mola_1);
const auto diffj_g = m_vars->intrinsic_diff_coeff_aqueous(j)*mu_face*flux_geom_coeff*mola_face*(gamma_0 - gamma_1);
diff_0 += nu_ji*(+ diffj_c + diffj_g);
diff_1 += nu_ji*(- diffj_c - diffj_g);
const scalar_t charge = m_data->charge_aqueous(j);
if (charge != 0.0) {
const scalar_t coeff_np = nu_ji*m_vars->intrinsic_diff_coeff_aqueous(j)*mu_face*charge;
if (charge > 0) {
np_plus_0 += coeff_np*mola_0;
np_plus_1 += coeff_np*mola_1;
} else {
np_minus_0 += coeff_np*mola_0;
np_minus_1 += coeff_np*mola_1;
}
}
}
// grad Psi
if (use_grad_psi) {
const scalar_t f_o_rt = m_vars->faraday()/m_vars->rt();
const scalar_t psi_0 = edisplacement(m_vars->edof_potential(0));
const scalar_t psi_1 = edisplacement(m_vars->edof_potential(1));
const scalar_t grad_psi =( psi_0 - psi_1 )*flux_geom_coeff;
if (grad_psi > 0) {
diff_0 += -(-np_minus_0-np_plus_1)*grad_psi*f_o_rt;
diff_1 += -(+np_minus_0+np_plus_1)*grad_psi*f_o_rt;
} else if (grad_psi < 0) {
diff_0 += (+np_minus_1+np_plus_0)*grad_psi*f_o_rt;
diff_1 += (-np_minus_1-np_plus_0)*grad_psi*f_o_rt;
}
}
transport_0 += m_vars->rho_w()*diff_0;
transport_1 += m_vars->rho_w()*diff_1;
}
courant_face += m_data->charge_component(i)*m_vars->faraday()*(transport_0 - transport_1);
m_current(element) = courant_face;
eresidual(edof_0) = (velocity_0 - transport_0)/m_scaling(edof_0);
eresidual(edof_1) = (velocity_1 - transport_1)/m_scaling(edof_0);
if (use_chemistry_rate) {
eresidual(edof_0) -= mass_coeff_0*m_vars->chemistry_rate(dof_0)/m_scaling(edof_0);
eresidual(edof_1) -= mass_coeff_1*m_vars->chemistry_rate(dof_1)/m_scaling(edof_0); //scaling vector is ndofx1
}
}
// potential
const auto edof_0 = m_vars->edof_potential(0);
const auto edof_1 = m_vars->edof_potential(1);
eresidual(edof_0) = - courant_face/m_scaling(edof_0);
eresidual(edof_1) = + courant_face/m_scaling(edof_0);
}
void ChlorideProgram::compute_jacobian(
Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian,
scalar_t alphadt
)
{
dfpm::list_triplet_t t_jacob;
const index_t estimation = m_mesh->nb_nodes()*(m_vars->ndof()*m_vars->ndof());
t_jacob.reserve(estimation);
#ifdef SPECMICP_HAVE_OPENMP
#pragma omp parallel for schedule(dynamic, 5)
#endif
for (index_t element=0; element<m_mesh->nb_elements(); ++element)
{
auto evel = m_vars->get_element_dofs(element, velocity);
auto edisp = m_vars->get_element_dofs(element, displacement);
compute_element_jacobian(element, edisp, evel, t_jacob, alphadt);
}
jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
jacobian.setFromTriplets(t_jacob.begin(), t_jacob.end());
//Eigen::saveMarket(jacobian, "jacob.mkt");
}
AdimensionalSystemSolution ChlorideProgram::perturb_adim_solution(
index_t element, index_t enode, const Vector& edisplacement)
{
AdimensionalSystemSolution sol;
specmicp_assert(m_chem_stagger);
auto node = m_mesh->get_node(element, enode);
ReturnCode retcode = m_chem_stagger->perturb_adim_solution(
node, m_vars.get(),
edisplacement.segment(enode*m_vars->ndof(), m_vars->ndof()),
m_vars->adim_solution(node), sol
);
if (retcode < ReturnCode::Success) {
throw std::runtime_error("Failed to perturb chemistry solution - Jacobian");
}
return sol;
}
AdimensionalSystemSolution ChlorideProgram::perturb_adim_solution(
index_t node, const Eigen::Ref<const Vector>& displacement)
{
AdimensionalSystemSolution sol;
specmicp_assert(m_chem_stagger);
ReturnCode retcode = m_chem_stagger->perturb_adim_solution(
node, m_vars.get(), displacement,
m_vars->adim_solution(node), sol
);
if (retcode < ReturnCode::Success) {
throw std::runtime_error("Failed to perturb chemistry solution - Residual");
}
return sol;
}
void ChlorideProgram::compute_element_jacobian(
index_t element,
Vector& edisplacement,
Vector& evelocity,
dfpm::list_triplet_t& jacobian,
scalar_t alphadt
)
{
Vector i_res(m_vars->edof());
Vector p_res(m_vars->edof());
std::array<AdimensionalSystemSolution, 2> adim_sols = {
m_vars->adim_solution(m_mesh->get_node(element, 0)),
m_vars->adim_solution(m_mesh->get_node(element, 1))
};
compute_element_residuals(element, edisplacement, evelocity,
adim_sols[0],
adim_sols[1],
i_res, NO_CHEMISTRY_RATE, true);
for (auto enode: m_mesh->range_nen())
{
const auto node = m_mesh->get_node(element, enode);
// concentration
for (auto component: m_data->range_aqueous_component())
{
const auto dof = m_vars->dof_component(node, component);
const auto ideq = m_ideq(dof);
if (ideq == no_equation) continue;
const auto edof = m_vars->edof_component(enode, component);
auto tmp = evelocity(edof);
auto tmpd = edisplacement(edof);
auto dh = eps_jacobian*std::abs(tmp);
if (dh < eps_jacobian) dh = eps_jacobian;
evelocity(edof) = tmp + dh;
dh = evelocity(edof) - tmp;
edisplacement(edof) = m_vars->predictors()(dof) + alphadt*evelocity(edof);
adim_sols[enode] = perturb_adim_solution(element, enode, edisplacement);
compute_element_residuals(element, edisplacement, evelocity,
adim_sols[0],
adim_sols[1],
p_res, NO_CHEMISTRY_RATE, true);
evelocity(edof) = tmp;
edisplacement(edof) = tmpd;
adim_sols[enode] = m_vars->adim_solution(node);
for (auto ernode: m_mesh->range_nen())
{
for (auto rndof: RangeIterator<index_t>(m_vars->ndof()))
{
const auto rdof = m_vars->ndof_to_dof(m_mesh->get_node(element, ernode), rndof);
const auto rideq = m_ideq(rdof);
if (rideq == no_equation) continue;
const auto redof = m_vars->ndof_to_edof(ernode, rndof);
#ifdef SPECMICP_HAVE_OPENMP
#pragma omp critical
{
#endif
jacobian.push_back(dfpm::triplet_t(rideq, ideq, (p_res(redof) - i_res(redof))/dh));
#ifdef SPECMICP_HAVE_OPENMP
}
#endif
}
}
}
// potential
{
const auto dof = m_vars->dof_potential(node);
const auto ideq = m_ideq(dof);
if (ideq == no_equation) continue;
const auto edof = m_vars->edof_potential(enode);
scalar_t tmp = edisplacement(edof);
scalar_t h = eps_jacobian*std::abs(tmp);
if (h < 1e-4*eps_jacobian) h = 1e-4*eps_jacobian;
edisplacement(edof) = tmp + h;
compute_element_residuals(element, edisplacement, evelocity,
adim_sols[0],
adim_sols[1],
p_res, NO_CHEMISTRY_RATE, true);
edisplacement(edof) = tmp;
for (auto enode: m_mesh->range_nen())
{
for (auto rndof: RangeIterator<index_t>(m_vars->ndof()))
{
const auto rdof = m_vars->ndof_to_dof(m_mesh->get_node(element, enode), rndof);
const auto rideq = m_ideq(rdof);
if (rideq == no_equation) continue;
const auto redof = m_vars->ndof_to_edof(enode, rndof);
auto value = (p_res(redof) - i_res(redof))/(h);
#ifdef SPECMICP_HAVE_OPENMP
#pragma omp critical
{
#endif
jacobian.push_back(dfpm::triplet_t(rideq, ideq, value));
#ifdef SPECMICP_HAVE_OPENMP
}
#endif
}
}
}
}
}
void ChlorideProgram::update_solution(
const Vector& update,
scalar_t lambda,
scalar_t alpha_dt,
Vector& predictor,
Vector& displacement,
Vector& velocity)
{
for (index_t node: m_mesh->range_nodes())
{
for (index_t i: m_data->range_aqueous_component())
{
const index_t dof = m_vars->dof_component(node, i);
const index_t id = id_equation(dof);
if (id == no_equation) continue;
velocity(dof) += lambda*update(id);
}
{
const index_t dofp = m_vars->dof_potential(node);
const index_t idp = id_equation(dofp);
if (idp != no_equation) {
velocity(dofp) += lambda*update(idp)/alpha_dt;
}
}
}
m_vars->velocities() = velocity;
displacement = predictor + alpha_dt*velocity;
m_vars->main_variables() = displacement;
#ifdef SPECMICP_HAVE_OPENMP
#pragma omp parallel for schedule(dynamic, 5)
#endif
for (index_t node=0; node<m_mesh->nb_nodes(); ++node) {
m_vars->adim_solution(node) =
perturb_adim_solution(node, m_vars->get_nodal_dofs(node, displacement));
}
}
void ChlorideProgram::compute_transport_rate(scalar_t dt, const Vector &displacement)
{
for (index_t node: m_mesh->range_nodes())
{
for (index_t i: m_data->range_aqueous_component())
{
const auto dof = m_vars->dof_component(node, i);
const auto ideq = id_equation(dof);
if (ideq == no_equation) continue;
const scalar_t transient = (
((m_vars->porosity(node)*displacement(dof)) -
(m_vars->predictor_porosity(node)*m_vars->predictors()(dof)))/dt
);
const scalar_t chem_rates = m_vars->chemistry_rate(dof);
m_vars->transport_rate(dof) = transient - chem_rates;
}
}
}
void ChlorideProgram::register_chem_stagger(
std::shared_ptr<ChlorideChemistryStaggerBase> chem_stagger)
{
m_chem_stagger = chem_stagger;
}
} // namespace chloride
} // namespace systems
} // namespace reactmicp
} // namespace specmicp

Event Timeline

c4science · Help