Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/examples/reactmicp/diffusion.cpp
030defdf3929reactmicp/unsaturated

Last Change
Raw File

diffusion.cpp
View Options

#include "dfpmsolver/parabolic_driver.hpp"
#include "dfpm/meshes/uniform_mesh1d.hpp"
#include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_program.hpp"
#include "reactmicp/systems/saturated_diffusion/variables.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp"
#include "specmicp/reaction_path.hpp"
void diffusion_test()
{
specmicp::database::Database dbhandler("../data/cemdata_specmicp.js");
std::vector<std::string> to_keep = {"H2O", "H[+]", "HCO3[-]"};
dbhandler.keep_only_components(to_keep);
specmicp::RawDatabasePtr thedatabase = dbhandler.get_database();
int nb_nodes = 100+1;
double dx = 0.0025;
double radius = 3.5+dx*1;
double height = 8;
//specmicp::reactmicp::mesh::Mesh1DPtr themesh = specmicp::reactmicp::mesh::uniform_mesh1d(nb_nodes, 0.0050, 5);
specmicp::reactmicp::mesh::Mesh1DPtr themesh =
specmicp::reactmicp::mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, height);
auto parameters = std::make_shared<
specmicp::reactmicp::systems::siasaturated::SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
specmicp::reactmicp::systems::siasaturated::SIASaturatedVariables variables(thedatabase, themesh, parameters);
specmicp::reactmicp::systems::siasaturated::SaturatedDiffusionProgram the_program(themesh, thedatabase, parameters);
specmicp::Vector& concentrations = variables.total_mobile_concentrations();
concentrations.setOnes();
concentrations.block(0, 0, the_program.get_ndf(), 1).setZero();
for (int node=0; node<50+1; ++node)
{
concentrations(variables.get_dof_total_concentration(node, 1)) = 0.0;
concentrations(variables.get_dof_total_concentration(node, 2)) = 0.0;
parameters->diffusion_coefficient(node) = 1e-5;
}
the_program.number_equations({0,});
specmicp::dfpmsolver::ParabolicDriver<
specmicp::reactmicp::systems::siasaturated::SaturatedDiffusionProgram> solver(the_program);
solver.get_options().sparse_solver = specmicp::SparseSolver::GMRES;
solver.get_options().alpha = 1;
double initial_sum = 0;
for (specmicp::index_t node: themesh->range_nodes())
{
initial_sum += concentrations(variables.get_dof_total_concentration(node, 2))*themesh->get_volume_cell(node);
}
double explicit_dt = std::pow(themesh->get_dx(0),2)/parameters->diffusion_coefficient(0);
std::cout << "Explicit dt : " << explicit_dt << std::endl;
double dt= 2000; //0.1*explicit_dt;
double total = 0;
double target = 10*24*3600;
int k=0;
std::cout << 0.0 << "\t";
std::cout << initial_sum << "\t" << 0.0 << "\t" << std::endl;
while (total < target)
{
//std::cout << total << std::endl;
specmicp::dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(dt, concentrations);
//specmicp_assert(retcode == specmicp::dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
total += dt;
if (k % 10 == 0)
{
std::cout << total/24.0/3600.0 << "\t";
double sum = 0;
for (specmicp::index_t node: themesh->range_nodes())
{
sum += concentrations(variables.get_dof_total_concentration(node, 2))*themesh->get_volume_cell(node);
}
std::cout << sum << "\t" << (initial_sum - sum)/(2*3.5*M_PI*height*1e-4) << "\t" << std::endl;
}
++k;
}
}
int main()
{
diffusion_test();
}
c4science · Help