Applications/ApplicationsLib/ProjectData.cpp

@@ -85,6 +85,9 @@
 #ifdef OGS_BUILD_PROCESS_PHASEFIELD
 #include "ProcessLib/PhaseField/CreatePhaseFieldProcess.h"
 #endif
+#ifdef OGS_BUILD_PROCESS_PHASEFIELDACID
+#include "ProcessLib/PhaseFieldAcid/CreatePhaseFieldAcidProcess.h"
+#endif
 #ifdef OGS_BUILD_PROCESS_RICHARDSCOMPONENTTRANSPORT
 #include "ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.h"
 #endif
@@ -818,6 +821,16 @@ void ProjectData::parseProcesses(
         }
         else
 #endif
+#ifdef OGS_BUILD_PROCESS_PHASEFIELDACID
+            if (type == "PHASE_FIELD_ACID")
+        {
+            process = ProcessLib::PhaseFieldAcid::createPhaseFieldAcidProcess(
+                name, *_mesh_vec[0], std::move(jacobian_assembler),
+                _process_variables, _parameters, _local_coordinate_system,
+                integration_order, process_config);
+        }
+        else
+#endif
 #ifdef OGS_BUILD_PROCESS_RICHARDSCOMPONENTTRANSPORT
             if (type == "RichardsComponentTransport")
         {

NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h

@@ -8,6 +8,8 @@
  *
  */
 
+#pragma once
+
 #include <cassert>
 
 #include "MathLib/Point3d.h"

NumLib/Fem/InitShapeMatrices.h

@@ -11,6 +11,7 @@
 
 #include <vector>
 
+#include "CoordinatesMapping/NaturalNodeCoordinates.h"
 #include "FiniteElement/TemplateIsoparametric.h"
 #include "MeshLib/Elements/Element.h"
 
@@ -67,6 +68,29 @@ initShapeMatrices(MeshLib::Element const& e, bool const is_axially_symmetric,
         e, is_axially_symmetric, points);
 }
 
+
+template <typename ShapeFunction, typename ShapeMatricesType, int GlobalDim,
+          ShapeMatrixType SelectedShapeMatrixType = ShapeMatrixType::ALL>
+std::vector<typename ShapeMatricesType::ShapeMatrices,
+            Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
+initShapeMatricesInNodes(MeshLib::Element const& e,
+                         bool const is_axially_symmetric)
+{
+    std::vector<MathLib::Point3d> points;
+
+    unsigned const n_nodes = e.getNumberOfNodes();
+    for (unsigned i = 0; i < n_nodes; ++i)
+    {
+        points.emplace_back(
+            NumLib::NaturalCoordinates<
+                typename ShapeFunction::MeshElement>::coordinates[i]);
+    }
+
+    return computeShapeMatrices<ShapeFunction, ShapeMatricesType, GlobalDim,
+                                SelectedShapeMatrixType>(
+        e, is_axially_symmetric, points);
+}
+
 template <typename ShapeFunction, typename ShapeMatricesType>
 double interpolateXCoordinate(
     MeshLib::Element const& e,

ProcessLib/PhaseFieldAcid/CMakeLists.txt

+get_source_files(SOURCES)
+
+ogs_add_library(PhaseFieldAcid ${SOURCES})
+target_link_libraries(PhaseFieldAcid PUBLIC ProcessLib PRIVATE ParameterLib)
+
+if(OGS_BUILD_TESTING)
+    include(Tests.cmake)
+endif()

ProcessLib/PhaseFieldAcid/CreatePhaseFieldAcidProcess.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreatePhaseFieldAcidProcess.h"
+
+#include <cassert>
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "ParameterLib/Parameter.h"
+#include "ParameterLib/Utils.h"
+#include "PhaseFieldAcidProcess.h"
+#include "PhaseFieldAcidProcessData.h"
+#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+namespace PhaseFieldAcid
+{
+std::unique_ptr<Process> createPhaseFieldAcidProcess(
+    std::string name,
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    int const integration_order,
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{prj__processes__process__type}
+    config.checkConfigParameter("type", "PHASE_FIELD_ACID");
+    DBUG("Create PhaseFieldAcidProcess.");
+
+    auto const coupling_scheme =
+        //! \ogs_file_param{prj__processes__process__PHASE_FIELD_ACID__coupling_scheme}
+        config.getConfigParameterOptional<std::string>("coupling_scheme");
+    const bool use_monolithic_scheme =
+        !(coupling_scheme && (*coupling_scheme == "staggered"));
+
+    // Process variable.
+
+    //! \ogs_file_param{prj__processes__process__PHASE_FIELD_ACID__process_variables}
+    auto const pv_config = config.getConfigSubtree("process_variables");
+
+    const int _concentration_process_id = 0;
+    const int _phasefield_process_id = 1;
+    ProcessVariable* variable_c;
+    ProcessVariable* variable_ph;
+
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+        process_variables;
+    if (use_monolithic_scheme)  // monolithic scheme.
+    {
+        OGS_FATAL("Monolithic implementation is not available.");
+    }
+    else  // staggered scheme.
+    {
+        using namespace std::string_literals;
+        for (
+            auto const& variable_name :
+            {//! \ogs_file_param_special{prj__processes__process__PHASE_FIELD_ACID__process_variables__concentration}
+             "concentration"s,
+             //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD_ACID__process_variables__phasefield}
+             "phasefield"s})
+        {
+            auto per_process_variables =
+                findProcessVariables(variables, pv_config, {variable_name});
+            process_variables.push_back(std::move(per_process_variables));
+        }
+        variable_c = &process_variables[0][0].get();
+        variable_ph = &process_variables[1][0].get();
+    }
+
+    DBUG("Associate concentration with process variable '{:s}'.",
+         variable_c->getName());
+
+    if (variable_c->getNumberOfGlobalComponents() != 1)
+    {
+        OGS_FATAL(
+            "Phasefield process variable '{}' is not a scalar variable but has "
+            "{} components.",
+            variable_c->getName(),
+            variable_c->getNumberOfGlobalComponents());
+    }
+
+    DBUG("Associate phase field with process variable '{:s}'.",
+         variable_ph->getName());
+    if (variable_ph->getNumberOfGlobalComponents() != 1)
+    {
+        OGS_FATAL(
+            "Phasefield process variable '{:s}' is not a scalar variable but "
+            "has {:d} components.",
+            variable_ph->getName(),
+            variable_ph->getNumberOfGlobalComponents());
+    }
+    auto const phasefieldacid_parameters_config =
+        //! \ogs_file_param{prj__processes__process__PHASE_FIELD_ACID__phasefieldacid_parameters}
+        config.getConfigSubtree("phasefieldacid_parameters");
+
+    // Chemical diffusivity
+    auto& chemical_diffusivity = ParameterLib::findParameter<double>(
+        phasefieldacid_parameters_config,
+        //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD_ACID__phasefieldacid_parameters__checmial_diffusivity}
+        "chemical_diffusivity", parameters, 1);
+    DBUG("Use '{}' as chemical diffusivity.", chemical_diffusivity.name);
+
+    // alpha
+    auto& alpha = ParameterLib::findParameter<double>(
+        phasefieldacid_parameters_config,
+        //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD_ACID__phasefieldacid_parameters__alpha}
+        "alpha", parameters, 1);
+    DBUG("Use '{}' as alpha.", alpha.name);
+
+    // rrc(reaction rate coefficient, k)
+    auto& rrc = ParameterLib::findParameter<double>(
+        phasefieldacid_parameters_config,
+        //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD_ACID__phasefieldacid_parameters__rrc}
+        "rrc", parameters, 1);
+    DBUG("Use '{}' as rrc.", rrc.name);
+
+    // epsilon
+    auto& epsilon = ParameterLib::findParameter<double>(
+        phasefieldacid_parameters_config,
+        //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD_ACID__phasefieldacid_parameters__epsilon}
+        "epsilon", parameters, 1);
+    DBUG("Use '{}' as epsi.", epsilon.name);
+
+    // tau
+    auto& tau = ParameterLib::findParameter<double>(
+        phasefieldacid_parameters_config,
+        //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD_ACID__phasefieldacid_parameters__tau}
+        "tau", parameters, 1);
+    DBUG("Use '{}' as tau.", tau.name);
+
+    // grad_phi_cutoff
+    auto const grad_phi_cutoff =
+        //! \ogs_file_param_special{prj__processes__process__PHASE_FIELD_ACID__phasefieldacid_parameters__grad_phi_cutoff}
+        phasefieldacid_parameters_config.getConfigParameter<double>(
+            "grad_phi_cutoff");
+
+    // Specific body force parameter.
+    Eigen::VectorXd specific_body_force;
+    std::vector<double> const b =
+        //! \ogs_file_param{prj__processes__process__PHASE_FIELD_ACID__specific_body_force}
+        config.getConfigParameter<std::vector<double>>("specific_body_force");
+    if (b.size() != mesh.getDimension())
+    {
+        OGS_FATAL(
+            "specific body force (gravity vector) has {} components, mesh "
+            "dimension is {}",
+            b.size(), mesh.getDimension());
+    }
+    specific_body_force.resize(b.size());
+    std::copy_n(b.data(), b.size(), specific_body_force.data());
+
+    PhaseFieldAcidProcessData process_data{materialIDs(mesh),
+                                           specific_body_force,
+                                           chemical_diffusivity,
+                                           alpha,
+                                           rrc,
+                                           epsilon,
+                                           tau,
+                                           grad_phi_cutoff};
+
+    SecondaryVariableCollection secondary_variables;
+
+    ProcessLib::createSecondaryVariables(config, secondary_variables);
+
+    return std::make_unique<PhaseFieldAcidProcess>(
+        std::move(name), mesh, std::move(jacobian_assembler), parameters,
+        integration_order, std::move(process_variables),
+        std::move(process_data), std::move(secondary_variables),
+        use_monolithic_scheme, _concentration_process_id,
+        _phasefield_process_id);
+}
+
+}  // namespace PhaseFieldAcid
+}  // namespace ProcessLib

ProcessLib/PhaseFieldAcid/CreatePhaseFieldAcidProcess.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace MeshLib
+{
+class Mesh;
+}
+namespace ParameterLib
+{
+struct CoordinateSystem;
+struct ParameterBase;
+}  // namespace ParameterLib
+namespace ProcessLib
+{
+class AbstractJacobianAssembler;
+class Process;
+class ProcessVariable;
+}  // namespace ProcessLib
+
+namespace ProcessLib
+{
+namespace PhaseFieldAcid
+{
+std::unique_ptr<Process> createPhaseFieldAcidProcess(
+    std::string name,
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    int const integration_order,
+    BaseLib::ConfigTree const& config);
+
+}  // namespace PhaseFieldAcid
+}  // namespace ProcessLib

ProcessLib/PhaseFieldAcid/IntegrationPointData.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+#include <Eigen/Core>
+
+namespace ProcessLib
+{
+namespace PhaseFieldAcid
+{
+template <typename ShapeMatrixType>
+struct IntegrationPointData final
+{
+    IntegrationPointData(
+        typename ShapeMatrixType::NodalRowVectorType const N_,
+        typename ShapeMatrixType::GlobalDimNodalMatrixType const dNdx_,
+        double const integration_weight_)
+        : N(N_),
+          dNdx(dNdx_),
+          integration_weight(integration_weight_),
+          mass_operator{N.transpose() * N * integration_weight}
+    {
+    }
+
+    typename ShapeMatrixType::NodalRowVectorType const N;
+    typename ShapeMatrixType::GlobalDimNodalMatrixType const dNdx;
+    double const integration_weight;
+    typename ShapeMatrixType::NodalMatrixType const mass_operator;
+
+    double dummy;
+    double dummy_prev;
+    double kappa;
+
+    void pushBackState() { dummy_prev = dummy; }
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+}  // namespace PhaseFieldAcid
+}  // namespace ProcessLib

ProcessLib/PhaseFieldAcid/LocalAssemblerInterface.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <vector>
+
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+
+namespace ProcessLib
+{
+namespace PhaseFieldAcid
+{
+struct PhaseFieldAcidLocalAssemblerInterface
+    : public ProcessLib::LocalAssemblerInterface,
+      public NumLib::ExtrapolatableElement
+{
+};
+
+}  // namespace PhaseFieldAcid
+}  // namespace ProcessLib

ProcessLib/PhaseFieldAcid/PhaseFieldAcidProcess.cpp

+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "PhaseFieldAcidProcess.h"
+
+#include <cassert>
+
+#include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "PhaseFieldAcidFEM.h"
+#include "ProcessLib/Process.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+namespace PhaseFieldAcid
+{
+PhaseFieldAcidProcess::PhaseFieldAcidProcess(
+    std::string name, MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+        process_variables,
+    PhaseFieldAcidProcessData&& process_data,
+    SecondaryVariableCollection&& secondary_variables,
+    bool const use_monolithic_scheme, const int concentration_process_id,
+    const int phasefield_process_id)
+    : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
+              integration_order, std::move(process_variables),
+              std::move(secondary_variables), use_monolithic_scheme),
+      _process_data(std::move(process_data)),
+      _concentration_process_id(concentration_process_id),
+      _phasefield_process_id(phasefield_process_id)
+{
+    if (use_monolithic_scheme)
+    {
+        OGS_FATAL(
+            "Monolithic scheme is not implemented for the PhaseField process.");
+    }
+}
+
+bool PhaseFieldAcidProcess::isLinear() const
+{
+    return false;
+}
+
+MathLib::MatrixSpecifications PhaseFieldAcidProcess::getMatrixSpecifications(
+    const int process_id) const
+{
+    auto const& l = *_local_to_global_index_map_single_component;
+    return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+            &l.getGhostIndices(), &_sparsity_pattern_with_single_component};
+}
+
+NumLib::LocalToGlobalIndexMap const& PhaseFieldAcidProcess::getDOFTable(
+    const int process_id) const
+{
+    // For the equation of phasefield
+    return *_local_to_global_index_map_single_component;
+}
+
+void PhaseFieldAcidProcess::constructDofTable()
+{
+    const int concentration_process_id = 0;
+    constructDofTableOfSpecifiedProcessStaggeredScheme(
+        concentration_process_id);
+
+    // TODO move the two data members somewhere else.
+    // for extrapolation of secondary variables of stress or strain
+    std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
+        *_mesh_subset_all_nodes};
+    _local_to_global_index_map_single_component =
+        std::make_unique<NumLib::LocalToGlobalIndexMap>(
+            std::move(all_mesh_subsets_single_component),
+            // by location order is needed for output
+            NumLib::ComponentOrder::BY_LOCATION);
+
+    assert(_local_to_global_index_map_single_component);
+
+    // For phase field equation.
+    _sparsity_pattern_with_single_component = NumLib::computeSparsityPattern(
+        *_local_to_global_index_map_single_component, _mesh);
+}
+
+void PhaseFieldAcidProcess::initializeConcreteProcess(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    MeshLib::Mesh const& mesh,
+    unsigned const integration_order)
+{
+    const int process_id = 0;
+    ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+    ProcessLib::createLocalAssemblers<PhaseFieldAcidLocalAssembler>(
+        mesh.getDimension(), mesh.getElements(), dof_table,
+        pv.getShapeFunctionOrder(), _local_assemblers,
+        mesh.isAxiallySymmetric(), integration_order, _process_data,
+        _concentration_process_id, _phasefield_process_id);
+
+    _process_data.element_kappa = MeshLib::getOrCreateMeshProperty<double>(
+        const_cast<MeshLib::Mesh&>(mesh), "kappa_avg",
+        MeshLib::MeshItemType::Cell, 1);
+    _process_data.element_h = MeshLib::getOrCreateMeshProperty<double>(
+        const_cast<MeshLib::Mesh&>(mesh), "h_avg", MeshLib::MeshItemType::Cell,
+        1);
+    _process_data.element_g = MeshLib::getOrCreateMeshProperty<double>(
+        const_cast<MeshLib::Mesh&>(mesh), "g_avg", MeshLib::MeshItemType::Cell,
+        1);
+    _process_data.element_f = MeshLib::getOrCreateMeshProperty<double>(
+        const_cast<MeshLib::Mesh&>(mesh), "f_avg", MeshLib::MeshItemType::Cell,
+        1);
+
+    // Initialize local assemblers after all variables have been set.
+    GlobalExecutor::executeMemberOnDereferenced(
+        &LocalAssemblerInterface::initialize, _local_assemblers,
+        *_local_to_global_index_map);
+}
+
+void PhaseFieldAcidProcess::initializeBoundaryConditions()
+{
+    // Staggered scheme:
+    // for the equations of deformation.
+    const int concentration_process_id = 0;
+    initializeProcessBoundaryConditionsAndSourceTerms(
+        *_local_to_global_index_map, concentration_process_id);
+    // for the phase field
+    const int phasefield_process_id = 1;
+    initializeProcessBoundaryConditionsAndSourceTerms(
+        *_local_to_global_index_map_single_component, phasefield_process_id);
+}
+
+void PhaseFieldAcidProcess::assembleConcreteProcess(
+    const double t, double const dt, std::vector<GlobalVector*> const& x,
+    std::vector<GlobalVector*> const& xdot, int const process_id,
+    GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+{
+    DBUG("Assemble PhaseFieldAcidProcess.");
+
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
+    switch (process_id)
+    {
+        case 1:
+        {
+            DBUG(
+                "Assemble the equations of phase-field in "
+                "PhaseFieldAcidProcess for the staggered scheme.");
+            break;
+        }
+        case 0:
+        {
+            DBUG(
+                "Assemble the equations of concentration in "
+                "PhaseFieldAcidProcess for the staggered scheme.");
+            break;
+        }
+    }
+
+    setCoupledSolutionsOfPreviousTimeStep();
+    dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+    dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+
+    ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeSelectedMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+        pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
+        b);
+}
+
+void PhaseFieldAcidProcess::assembleWithJacobianConcreteProcess(
+    const double t, double const dt, std::vector<GlobalVector*> const& x,
+    std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
+    const double dx_dx, int const process_id, GlobalMatrix& M, GlobalMatrix& K,
+    GlobalVector& b, GlobalMatrix& Jac)
+{
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
+
+    // For the staggered scheme
+    if (process_id == _concentration_process_id)
+    {
+        DBUG(
+            "Assemble the Jacobian equations of concentration in "
+            "PhaseFieldAcidProcess for the staggered scheme.");
+    }
+    else if (process_id == _phasefield_process_id)
+    {
+        DBUG(
+            "Assemble the Jacobian equations of phase-field in "
+            "PhaseFieldAcidProcess for the staggered scheme.");
+    }
+    else
+    {
+        OGS_FATAL("Unknown process id {}.");
+    }
+    dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+    dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+
+    ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeSelectedMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+        _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
+}
+
+void PhaseFieldAcidProcess::preTimestepConcreteProcess(
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
+    const int process_id)
+{
+    DBUG("PreTimestep PhaseFieldAcidProcess {}.", process_id);
+
+    /*  TODO (yoshioka) Maybe needed for xdot.
+    _x_previous_timestep =
+        MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+        */
+    assert(process_id < 2);
+
+    if (_use_monolithic_scheme)
+    {
+        return;
+    }
+
+    if (!_xs_previous_timestep[process_id])
+    {
+        _xs_previous_timestep[process_id] =
+            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+                *x[process_id]);
+    }
+    else
+    {
+        auto& x0 = *_xs_previous_timestep[process_id];
+        MathLib::LinAlg::copy(*x[process_id], x0);
+    }
+
+    auto& x0 = *_xs_previous_timestep[process_id];
+    MathLib::LinAlg::setLocalAccessibleVector(x0);
+    ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+
+    GlobalExecutor::executeSelectedMemberOnDereferenced(
+        &LocalAssemblerInterface::preTimestep, _local_assemblers,
+        pv.getActiveElementIDs(), getDOFTable(process_id), *x[process_id], t,
+        dt);
+}
+
+void PhaseFieldAcidProcess::postTimestepConcreteProcess(
+    std::vector<GlobalVector*> const& x, const double t, const double dt,
+    int const process_id)
+{
+    if (isPhaseFieldProcess(process_id))
+    {
+        DBUG("PostTimestep PhaseFieldProcess.");
+
+        std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
+        dof_tables.reserve(x.size());
+
+        std::generate_n(
+            std::back_inserter(dof_tables), x.size(),
+            [&]()
+            { return _local_to_global_index_map_single_component.get(); });
+
+        ProcessLib::ProcessVariable const& pv =
+            getProcessVariables(process_id)[0];
+
+        GlobalExecutor::executeSelectedMemberOnDereferenced(
+            &LocalAssemblerInterface::postTimestep, _local_assemblers,
+            pv.getActiveElementIDs(), dof_tables, x, t, dt);
+    }
+}
+
+void PhaseFieldAcidProcess::setCoupledSolutionsOfPreviousTimeStepPerProcess(
+    const int process_id)
+{
+    const auto& x_t0 = _xs_previous_timestep[process_id];
+    if (x_t0 == nullptr)
+    {
+        OGS_FATAL(
+            "Memory is not allocated for the global vector of the solution of "
+            "the previous time step for the staggered scheme.\n It can be done "
+            "by overriding Process::preTimestepConcreteProcess (ref. "
+            "HTProcess::preTimestepConcreteProcess) ");
+    }
+
+    _coupled_solutions->coupled_xs_t0[process_id] = x_t0.get();
+}
+
+void PhaseFieldAcidProcess::setCoupledSolutionsOfPreviousTimeStep()
+{
+    _coupled_solutions->coupled_xs_t0.resize(2);
+    setCoupledSolutionsOfPreviousTimeStepPerProcess(_concentration_process_id);
+    setCoupledSolutionsOfPreviousTimeStepPerProcess(_phasefield_process_id);
+}
+
+void PhaseFieldAcidProcess::postNonLinearSolverConcreteProcess(
+    GlobalVector const& x, GlobalVector const& /*xdot*/, const double t,
+    double const dt, const int process_id)
+{
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
+
+    dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+    dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+
+    ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+}
+
+void PhaseFieldAcidProcess::updateConstraints(GlobalVector& lower,
+                                              GlobalVector& upper,
+                                              int const process_id)
+{
+    if (process_id == _phasefield_process_id)
+    {
+        MathLib::LinAlg::set(lower, -1.0);
+        MathLib::LinAlg::set(upper, 1.0);
+        return;
+    }
+    if (process_id == _concentration_process_id)
+    {
+        MathLib::LinAlg::set(lower, 0);
+        MathLib::LinAlg::set(upper, 0.1);
+        return;
+    }
+    OGS_FATAL(
+        "Unknown process id {}. Must be {} for phasefield process or {} for "
+        "concentration process.",
+        process_id, _phasefield_process_id, _concentration_process_id);
+}
+
+constexpr bool PhaseFieldAcidProcess::isPhaseFieldProcess(
+    int const process_id) const
+{
+    return process_id == _phasefield_process_id;
+}
+
+}  // namespace PhaseFieldAcid
+}  // namespace ProcessLib

ProcessLib/PhaseFieldAcid/PhaseFieldAcidProcess.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "LocalAssemblerInterface.h"
+#include "PhaseFieldAcidProcessData.h"
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace PhaseFieldAcid
+{
+class PhaseFieldAcidProcess final : public Process
+{
+public:
+    PhaseFieldAcidProcess(
+        std::string name, MeshLib::Mesh& mesh,
+        std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+            jacobian_assembler,
+        std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+            parameters,
+        unsigned const integration_order,
+        std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+            process_variables,
+        PhaseFieldAcidProcessData&& process_data,
+        SecondaryVariableCollection&& secondary_variables,
+        bool const use_monolithic_scheme, const int concentration_process_id,
+        const int phasefield_process_id);
+
+    //! \name ODESystem interface
+    //! @{
+    bool isLinear() const override;
+    //! @}
+
+    MathLib::MatrixSpecifications getMatrixSpecifications(
+        const int process_id) const override;
+
+    NumLib::LocalToGlobalIndexMap const& getDOFTable(
+        const int process_id) const override;
+
+private:
+    using LocalAssemblerInterface = PhaseFieldAcidLocalAssemblerInterface;
+
+    void constructDofTable() override;
+
+    void initializeBoundaryConditions() override;
+
+    void initializeConcreteProcess(
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh,
+        unsigned const integration_order) override;
+
+    void assembleConcreteProcess(const double t, double const dt,
+                                 std::vector<GlobalVector*> const& x,
+                                 std::vector<GlobalVector*> const& xdot,
+                                 int const process_id, GlobalMatrix& M,
+                                 GlobalMatrix& K, GlobalVector& b) override;
+
+    void assembleWithJacobianConcreteProcess(
+        const double t, double const dt, std::vector<GlobalVector*> const& x,
+        std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
+        const double dx_dx, int const process_id, GlobalMatrix& M,
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
+                                    const int process_id) override;
+
+    void setCoupledSolutionsOfPreviousTimeStepPerProcess(const int process_id);
+
+    /// Set the solutions of the previous time step to the coupled term.
+    /// It is only for the staggered scheme, and it must be called within
+    /// the coupling loop because that the coupling term is only created there.
+    void setCoupledSolutionsOfPreviousTimeStep();
+
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                     const double t, const double delta_t,
+                                     int const process_id) override;
+
+    void postNonLinearSolverConcreteProcess(GlobalVector const& x,
+                                            GlobalVector const& xdot,
+                                            const double t, double const dt,
+                                            int const process_id) override;
+
+    void updateConstraints(GlobalVector& lower, GlobalVector& upper,
+                           int process_id) override;
+
+private:
+    PhaseFieldAcidProcessData _process_data;
+
+    std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap>
+        _local_to_global_index_map_single_component;
+
+    /// Sparsity pattern for the phase field equation, and it is initialized
+    ///  only if the staggered scheme is used.
+    GlobalSparsityPattern _sparsity_pattern_with_single_component;
+
+    /// Check whether the process represented by \c process_id is/has
+    /// mechanical process. In the present implementation, the mechanical
+    /// process has process_id == 0 in the staggered scheme.
+    constexpr bool isPhaseFieldProcess(int const process_id) const;
+
+    /// Solutions of the previous time step
+    std::array<std::unique_ptr<GlobalVector>, 2> _xs_previous_timestep;
+    const int _concentration_process_id;
+    const int _phasefield_process_id;
+};
+
+}  // namespace PhaseFieldAcid
+}  // namespace ProcessLib

ProcessLib/PhaseFieldAcid/PhaseFieldAcidProcessData.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+#pragma once
+
+#include <Eigen/Eigen>
+#include <memory>
+#include <utility>
+
+namespace MeshLib
+{
+template <typename T>
+class PropertyVector;
+}
+namespace ParameterLib
+{
+template <typename T>
+struct Parameter;
+}
+namespace ProcessLib::PhaseFieldAcid
+{
+struct PhaseFieldAcidProcessData
+{
+    MeshLib::PropertyVector<int> const* const material_ids = nullptr;
+    Eigen::VectorXd const specific_body_force;
+
+    ParameterLib::Parameter<double> const& chemical_diffusivity;
+    ParameterLib::Parameter<double> const& alpha;
+    ParameterLib::Parameter<double> const& rrc;
+    ParameterLib::Parameter<double> const& epsilon;
+    ParameterLib::Parameter<double> const& tau;
+    double const grad_phi_cutoff;
+
+    MeshLib::PropertyVector<double>* element_kappa = nullptr;
+    MeshLib::PropertyVector<double>* element_g = nullptr;
+    MeshLib::PropertyVector<double>* element_f = nullptr;
+    MeshLib::PropertyVector<double>* element_h = nullptr;
+};
+
+}  // namespace ProcessLib::PhaseFieldAcid

ProcessLib/PhaseFieldAcid/SecondaryData.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+#include <vector>
+
+namespace ProcessLib
+{
+namespace PhaseFieldAcid
+{
+/// Used for the extrapolation of the integration point values. It is
+/// ordered (and stored) by integration points.
+template <typename ShapeMatrixType>
+struct SecondaryData
+{
+    std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
+};
+
+}  // namespace PhaseFieldAcid
+}  // namespace ProcessLib

scripts/cmake/ProcessesSetup.cmake

@@ -11,6 +11,7 @@ set(_processes_list
     LIE
     ThermoRichardsMechanics
     PhaseField
+    PhaseFieldAcid
     RichardsComponentTransport
     RichardsFlow
     RichardsMechanics