Hi All,
I am trying to setup a "simple" solver for a Stokes-Cahn-Hilliard problem
but I've gotten stuck trying to implement the initial conditions for the
color/phase function. The problem looks like the following:
[image: Screenshot from 2023-11-18 11-04-57.png]
Where there's a regularized surface tension where the sharp interface would
normally be between two fluids. The separate fluids are tracked through the
function c. Using a forward Euler scheme, and introducing a dummy variable
for the curvature I end up with[image: Screenshot from 2023-11-18
11-06-46.png]
The eventual weak form isn't to important yet, since I can't even build the
initial conditions. At first I thought to extend Step-22 and use a block
system and solve "chunks" of. In particular, solving the stokes update and
then solving for the update to the phase function; but I can't figure out
how to set the initial conditions for $c^0$. Now I am wondering if the
correct approach is to use FETools and separate DoF handlers similar to
Step-35.
Any advice or direction would be greatly appreciated. I've included the
very barebones implementation I have right now.
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/c90bb603-d6bd-4f28-ba03-aa987ebe1cden%40googlegroups.com.
// Deal.II Libraries
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/types.h>
#include <deal.II/base/function.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/fe/mapping_q_generic.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparse_ilu.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools_project.h>
#include <ostream>
#include <random>
#include <string>
#include <unordered_map>
#include <fstream>
namespace stokes_cahn_hilliard {
using namespace dealii;
template<int dim>
class stokesCahnHilliard
{
public:
stokesCahnHilliard(
const std::unordered_map<std::string, double> params,
const double total_sim_time,
const unsigned int degree);
void run();
private:
void setupSystem(
const std::unordered_map<std::string, double> params
);
void addPeriodicToTriang();
void updateDoFHandler();
void reinitializeLinearSystem();
void initializeValues();
void assembleSystem();
const std::unordered_map<std::string, double> params;
const double total_sim_time;
const unsigned int degree;
const double eps;
const double eta;
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
std::vector<unsigned int> block_component;
BlockSparsityPattern block_sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
BlockVector<double> system_rhs;
BlockVector<double> solution;
double time_step;
unsigned int time_step_number;
unsigned int max_refine;
unsigned int min_refine;
};
template<int dim>
stokesCahnHilliard<dim> :: stokesCahnHilliard(
const std::unordered_map<std::string, double> params,
const double total_sim_time,
const unsigned int degree)
: params(params)
, total_sim_time(total_sim_time)
, degree(degree)
, eps(params.at("eps"))
, eta(params.at("eta"))
, fe(FE_Q<dim>(degree+1), dim,
FE_Q<dim>(degree), 1,
FE_Q<dim>(degree), 1,
FE_Q<dim>(degree), 1,
FE_Q<dim>(degree), 1)
, dof_handler(this->triangulation)
, time_step(1. / 512.)
, time_step_number(1)
, max_refine(12)
, min_refine(4)
{
std::cout << "Constructing" << std::endl;
std::cout << "Components: " << fe.n_components() << std::endl;
for(uint i = 0; i < fe.n_components(); i++){
if(i < dim) this->block_component.push_back(0);
if(i >= dim) this->block_component.push_back(i+1-dim);
}
}
template<int dim>
class InitialValuesC : public Function<dim>
{
public:
InitialValuesC(double eps);
virtual double value(
const Point<dim> &p,
const unsigned int component = 0
) const override;
private:
double eps;
};
template<int dim>
InitialValuesC<dim>::InitialValuesC(double eps)
: eps(eps)
{}
template<int dim>
double InitialValuesC<dim> :: value(
const Point<dim> &p,
const unsigned int /*component*/
) const
{
return std::tanh(
(p.norm()-0.25) / (2*std::sqrt(2)*this->eps)
);
}
template<int dim>
void stokesCahnHilliard<dim> :: addPeriodicToTriang()
{
std::cout << "Connecting nodes to neighbours due to horizontal "
<< "periodic boundary conditions"
<< std::endl;
if(dim == 2){
std::vector<GridTools::PeriodicFacePair<
typename Triangulation<dim>::cell_iterator>
> matchedPairsX;
std::vector<GridTools::PeriodicFacePair<
typename Triangulation<dim>::cell_iterator>
> matchedPairsY;
GridTools::collect_periodic_faces(this->triangulation,
0, 1, 0, matchedPairsX);
GridTools::collect_periodic_faces(this->triangulation,
2, 3, 1, matchedPairsY);
triangulation.add_periodicity(matchedPairsX);
triangulation.add_periodicity(matchedPairsY);
} else if(dim == 3){
std::vector<GridTools::PeriodicFacePair<
typename Triangulation<dim>::cell_iterator>
> matchedPairsX;
std::vector<GridTools::PeriodicFacePair<
typename Triangulation<dim>::cell_iterator>
> matchedPairsY;
std::vector<GridTools::PeriodicFacePair<
typename Triangulation<dim>::cell_iterator>
> matchedPairsZ;
GridTools::collect_periodic_faces(this->triangulation,
0, 1, 0, matchedPairsX);
GridTools::collect_periodic_faces(this->triangulation,
2, 3, 1, matchedPairsY);
GridTools::collect_periodic_faces(this->triangulation,
4, 5, 2, matchedPairsZ);
triangulation.add_periodicity(matchedPairsX);
triangulation.add_periodicity(matchedPairsY);
triangulation.add_periodicity(matchedPairsZ);
}
}
template<int dim>
void stokesCahnHilliard<dim> :: updateDoFHandler(){
std::cout << "Distributing degrees of freedom" << std::endl;
this->dof_handler.distribute_dofs(this->fe);
std::cout << "Introducing periodicity constraints and hanging node "
<< "constraints" << std::endl;
if(dim == 2){
std::vector<GridTools::PeriodicFacePair<
typename DoFHandler<dim>::cell_iterator>
> periodicity_vectorX;
std::vector<GridTools::PeriodicFacePair<
typename DoFHandler<dim>::cell_iterator>
> periodicity_vectorY;
GridTools::collect_periodic_faces(this->dof_handler,
0,1,0,periodicity_vectorX);
GridTools::collect_periodic_faces(this->dof_handler,
2,3,1,periodicity_vectorY);
DoFTools::make_periodicity_constraints<dim,dim>(periodicity_vectorX,
this->constraints);
DoFTools::make_periodicity_constraints<dim,dim>(periodicity_vectorY,
this->constraints);
DoFTools::make_hanging_node_constraints(this->dof_handler,
this->constraints);
} else if(dim == 3){
std::vector<GridTools::PeriodicFacePair<
typename DoFHandler<dim>::cell_iterator>
> periodicity_vectorX;
std::vector<GridTools::PeriodicFacePair<
typename DoFHandler<dim>::cell_iterator>
> periodicity_vectorY;
std::vector<GridTools::PeriodicFacePair<
typename DoFHandler<dim>::cell_iterator>
> periodicity_vectorZ;
GridTools::collect_periodic_faces(this->dof_handler,
0,1,0,periodicity_vectorX);
GridTools::collect_periodic_faces(this->dof_handler,
2,3,1,periodicity_vectorY);
GridTools::collect_periodic_faces(this->dof_handler,
4,5,2,periodicity_vectorZ);
DoFTools::make_periodicity_constraints<dim,dim>(periodicity_vectorX,
this->constraints);
DoFTools::make_periodicity_constraints<dim,dim>(periodicity_vectorY,
this->constraints);
DoFTools::make_periodicity_constraints<dim,dim>(periodicity_vectorZ,
this->constraints);
DoFTools::make_hanging_node_constraints(this->dof_handler,
this->constraints);
}
constraints.close();
std::cout << "Renumbering degrees of freedom" << std::endl;
DoFRenumbering::Cuthill_McKee(this->dof_handler);
DoFRenumbering::component_wise(this->dof_handler,
this->block_component);
const std::vector<types::global_dof_index> dofs_per_block =
DoFTools::count_dofs_per_fe_block(this->dof_handler,
this->block_component);
std::cout << "Total degrees of freedom: "
<< this->dof_handler.n_dofs() << std::endl;
for(uint i = 0; i < 5; i++){
std::cout << " Block " << i << ": " << dofs_per_block[i]
<< std::endl;
}
std::cout << "Building block sparsity pattern" << std::endl;
{
BlockDynamicSparsityPattern dsp(
dofs_per_block,
dofs_per_block
);
DoFTools::make_sparsity_pattern(this->dof_handler,
dsp,
this->constraints,
false);
this->block_sparsity_pattern.copy_from(dsp);
}
this->block_sparsity_pattern.compress();
std::cout << "Completed update to degrees of freedom and sparsity pattern"
<< std::endl;
}
template<int dim>
void stokesCahnHilliard<dim> :: reinitializeLinearSystem(){
std::cout << "Reinitalizing vectors and matrices" << std::endl;
const std::vector<types::global_dof_index> dofs_per_block =
DoFTools::count_dofs_per_fe_block(this->dof_handler,
this->block_component);
system_matrix.reinit(this->block_sparsity_pattern);
system_rhs.reinit(dofs_per_block);
solution.reinit(dofs_per_block);
}
template<int dim>
void stokesCahnHilliard<dim> :: initializeValues(){
std::cout << "Setting initial values for color function" << std::endl;
VectorTools::project(this->dof_handler,
this->constraints,
QGauss<dim>(2),
InitialValuesC<dim>(this->eps),
this->solution.block(3));
}
template<int dim>
void stokesCahnHilliard<dim> :: setupSystem(
const std::unordered_map<std::string, double> params){
std::cout << "Passing parameters:" << std::endl;
for(auto it=params.begin(); it!=params.end(); it++){
std::cout << " " << it->first
<< ": " << it->second
<< std::endl;
}
std::cout << "Building mesh" << std::endl;
GridGenerator::hyper_cube(this->triangulation,
-1, 1, true);
this->addPeriodicToTriang();
std::cout << "Refining grid" << std::endl;
this->triangulation.refine_global(6);
std::cout << "Mesh generated...\n"
<< " Active cells: " << triangulation.n_active_cells()
<< std::endl;
this->updateDoFHandler();
this->reinitializeLinearSystem();
this->initializeValues();
}
template<int dim>
void stokesCahnHilliard<dim> :: assembleSystem()
{
// Solve Stokes Equation
// Use the result to update Cahn-Hilliard
this->system_matrix = 0;
this->system_rhs = 0;
this->solution = 0;
}
template<int dim>
void stokesCahnHilliard<dim>::run(){
this->setupSystem(this->params);
}
}
int main(){
try{
std::cout << "Running" << std::endl << std::endl;
std::unordered_map<std::string, double> params;
params["eps"] = 1e-2;
params["eta"] = 1.;
double totalSimTime = 10;
stokes_cahn_hilliard::stokesCahnHilliard<2> stokesCahnHilliard(
params,
totalSimTime,
1
);
stokesCahnHilliard.run();
std::cout << "Completed" << std::endl;
return 0;
} catch(std::exception &exc) {
std::cerr << "--------------------------------------------------"
<< std::endl
<< "Exception while processing:" << std::endl
<< exc.what() << std::endl
<< "--------------------------------------------------"
<< std::endl;
} catch(...) {
std::cerr << "--------------------------------------------------"
<< "Unknown exception"
<< "--------------------------------------------------"
<< std::endl;
}
}
