Dear Wolfgang:
Thanks for your replying.I have tried some times to run my code on
the HPC. It seems that every application runs for a long time, this
error(fatal error in PMPI_Waitwall:See the MPI_ERROR field in MPI_Status
for the error code ) wil happen near "setup matrix-free （CPU/wall）" in the
file dealii-main.cc line 420.I think the MPI function MPI_BARRIER may help.
But I am confused that
how to use this function in my code. And this is the code How I realize my
PDEs in phase field of solidification. Another question is that when I try
to run my code in dimension3, it seems that this code run much slower than
dimension2. I didn‘ change other things just chang the dimension, is there
something went wrong
during this code.Thanks a lot!
Best Wishes!
On Saturday, April 1, 2023 at 12:44:06 AM UTC+8 Wolfgang Bangerth wrote:
> On 3/29/23 07:02, Tom Li wrote:
> > I am using matrix-free technology to solve a phase field problem.When I
> > compile my program everythiong went well, but something went wrong when
> the
> > program runs for a long time. It seems that this error has some
> relationships
> > with MPI;
>
> Tom:
> it's impossible to say from just the error message. It could be that your
> job
> ran out of time in the queue you submitted it to. It could be a hardware
> failure. It could be a software bug.
>
> The first step is to find out whether the problem is reproducible. If you
> run
> the exact same job, does it error out in the exact same place?
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>
--
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/9f382dc5-80e9-4e1c-9c93-42cd7accf3adn%40googlegroups.com.
##
# CMake script for the step-37 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "pftest6-10")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 3.3.0)
FIND_PACKAGE(deal.II 9.4.0
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
#
# Are all dependencies fulfilled?
#
IF(NOT DEAL_II_WITH_LAPACK) # keep in one line
MESSAGE(FATAL_ERROR "
Error! This tutorial requires a deal.II library that was configured with the
following options:
DEAL_II_WITH_LAPACK = ON
However, the deal.II library found at ${DEAL_II_PATH} was configured with these
options:
DEAL_II_WITH_LAPACK = ${DEAL_II_WITH_LAPACK}
This conflicts with the requirements."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
// alloy.h
#ifndef _parameters_h
#define _parameters_h
const double epsilon_4 = 0.0106;//0.0106; // anisotropy of interface energy
/** Parameters of Phase-field Model **/
const double a1 = 0.8839;
const double a2 = 0.6267;
const double D_l = 2.42e-05;
const double D_s = 0.0 * 1.15e-08;
const double pc = 0.15;
const double c_l0 = 4.0;
const double delta_T = 1.0; //过冷度
const double cooling_rate = 0.5/60; //冷却速率
const double T_m = 934.0; //纯铝的熔点
const double T_0 = 920.0; //初始的温度
const double m_l = (T_0 - T_m) / c_l0; //液相线的斜率
const double omega = 0.55;
const double NP_pf_upper = 10.26;
const double NP_pf_lower = -10.26;
const double pf_upper = tanh(NP_pf_upper/sqrt(2.0));
const double pf_lower = tanh(NP_pf_lower/sqrt(2.0));
const double Gamma = 2.41e-05;//[cm*K]吉布斯-汤姆森系数
const double delta_T_0 = m_l*(pc-1.0)*c_l0;//[K]freezing range
const double d0 = Gamma/delta_T_0;//[cm] chemical capillart length
const double xi = 0.2410/(d0*1.0e4);//ratio of interface thicknss to capillary length
const double lambda = a1*xi; //
const double w0 = d0*xi; // [cm] interface width
const double tau0 = a2*lambda*w0*w0/D_l;//[s] relaxation time
const double D = a1*a2*xi; // dimensionless solute diffusion coefficient in liquid (rescaled by w0^2/tau0)
const double CR=cooling_rate*tau0;//cooling rate rescaled by [1K/tau0]
const double amp_pf=0.15;//amplitude of noise for phase field
const double amp_c=0.00;//amplitude of noise for concentration field
//////////////////////////////////////////////////
const unsigned int degree_finite_element = 2;
const unsigned int dimension = 2;
const double _time_step = 0.2;
const double pi_2 = 1.57079632679489661923;
double tt =0;
const double length = 4150.0;
double undercooling = 0.0;//实时过冷度
using namespace dealii;
/*interface energy anisotropy*/
template <int dim, typename Number, typename VectorizedArrayType>
VectorizedArrayType
get_interface_isotropy(const Tensor<1,dim, VectorizedArrayType> &grad, VectorizedArrayType& square_grad)
{
VectorizedArrayType as_n = 0.0;
Number a_si = 0.0;
Number n_fourth_power = 0.0; // [(grad_x)^4 + (grad_y)^4 + (grad_z)^4]/(square_grad * square_grad)
//Point<dim> n_two_power_grad(0,0); // [(grad_x)^2 /(square_grad )
for (unsigned int v =0; v < VectorizedArray<Number>::size(); ++v)
{
n_fourth_power = 0.0;
if (square_grad[v] >= 1.0e-16)
{for (unsigned int d = 0; d < dim; ++ d)
{ n_fourth_power += pow(grad[d][v],4)/(square_grad[v]*square_grad[v]);
}
}
a_si=1.0-3.0*epsilon_4+4.0*epsilon_4*n_fourth_power;
as_n[v]=a_si*a_si;
}
return as_n;
}
template <int dim, typename Number, typename VectorizedArrayType>
VectorizedArrayType get_isotropy(const Tensor<1,dim, VectorizedArrayType> &grad)
{
VectorizedArrayType as_n = 0.0;
VectorizedArrayType square_grad = 0.0;
bool interface = false;
for (unsigned int d = 0; d < dim; ++ d)
square_grad += grad[d]*grad[d];
for (unsigned int v =0; v < VectorizedArray<Number>::size(); ++v)
if(square_grad[v]>1.0e-16)
{
interface=true;
break;
}
if(interface)
as_n = get_interface_isotropy<dim,Number, VectorizedArray<Number> >(grad, square_grad);
else
{ as_n = (1.0-3.0*epsilon_4)*(1.0-3.0*epsilon_4);
}
// if(tt >10.) std::cout<<"grad="<<grad<<":" <<as_n<<std::endl;
return as_n;
}
template <int dim, typename Number, typename VectorizedArrayType>
std::pair<VectorizedArrayType, Tensor<1,dim, VectorizedArrayType>> get_anisotropy(const Tensor<1,dim, VectorizedArrayType> &grad)
{
Tensor<1,dim, VectorizedArrayType> anisotropy;
VectorizedArrayType as_n = 0.0;
Number square_grad = 0.0;
Number a_si = 0.0;
Number n_fourth_power = 0.0; // [(grad_x)^4 + (grad_y)^4 + (grad_z)^4]/(square_grad * square_grad)
for (unsigned int v =0; v < VectorizedArray<Number>::size(); ++v)
{
square_grad = 0.0;
n_fourth_power = 0.0;
for (unsigned int d = 0; d < dim; ++d)
square_grad+= grad[d][v]*grad[d][v]; // square of the phase-field gradient
if (square_grad < 1.0e-16) n_fourth_power = 0.0;
else
{
for (unsigned int d = 0; d < dim; ++ d)
n_fourth_power += pow(grad[d][v],4)/(square_grad*square_grad);
}
a_si=1.0-3.0*epsilon_4+4.0*epsilon_4*n_fourth_power;
// std::cout<<"a_s1="<<a_si<<std::endl;
as_n[v]=a_si*a_si;
if (square_grad < 1.0e-16)
{
for (unsigned int d = 0; d < dim; ++ d)
{
anisotropy[d][v] = as_n[v];
}
}
else
{
for (unsigned int d = 0; d < dim; ++ d)
{
Number skew_anisotropy = 16.0*a_si*epsilon_4*(pow(grad[d][v],2)/square_grad
- n_fourth_power);
// std::cout<<"skew1="<<skew_anisotropy<<std::endl;
anisotropy[d][v] = skew_anisotropy;
}
}
}
return std::make_pair(as_n, anisotropy);
}
template <int dim, typename Number, typename VectorizedArrayType>
std::pair<VectorizedArrayType, Tensor<1,dim, VectorizedArrayType>>
get_interface_anisotropy(const Tensor<1,dim, VectorizedArrayType> &grad, VectorizedArrayType& square_grad)
{
Tensor<1,dim, VectorizedArrayType> anisotropy;
VectorizedArrayType as_n = 0.0;
Number a_si = 0.0;
Number n_fourth_power = 0.0; // [(grad_x)^4 + (grad_y)^4 + (grad_z)^4]/(square_grad * square_grad)
//Point<dim> n_two_power_grad(0,0); // [(grad_x)^2 /(square_grad )
for (unsigned int v =0; v < VectorizedArray<Number>::size(); ++v)
{
n_fourth_power = 0.0;
Point<dim> n_two_power_grad(0,0); // [(grad_x)^2 /(square_grad )
if (square_grad[v] >= 1.0e-16)
{for (unsigned int d = 0; d < dim; ++ d)
{ n_fourth_power += pow(grad[d][v],4)/(square_grad[v]*square_grad[v]);
n_two_power_grad[d] = pow(grad[d][v],2)/square_grad[v];
}
}
a_si=1.0-3.0*epsilon_4+4.0*epsilon_4*n_fourth_power;
as_n[v]=a_si*a_si;
for (unsigned int d = 0; d < dim; ++ d)
anisotropy[d][v] = as_n[v]+ 16.0*a_si*epsilon_4*(n_two_power_grad[d]
- n_fourth_power);
}
return std::make_pair(as_n, anisotropy);
}
template <int dim, typename Number, typename VectorizedArrayType>
std::pair<VectorizedArrayType, Tensor<1,dim, VectorizedArrayType>> get_anisotropy_vectorized(const Tensor<1,dim, VectorizedArrayType> &grad)
{
Tensor<1,dim, VectorizedArrayType> anisotropy;
VectorizedArrayType as_n = 0.0;
VectorizedArrayType square_grad = 0.0;
bool interface = false;
for (unsigned int d = 0; d < dim; ++ d)
square_grad += grad[d]*grad[d];
for (unsigned int v =0; v < VectorizedArray<Number>::size(); ++v)
if(square_grad[v]>1.0e-15)
{
interface=true;
break;
}
if(interface) return get_interface_anisotropy<dim,Number, VectorizedArray<Number> >(grad, square_grad);
else
{ as_n = (1.0-3.0*epsilon_4)*(1.0-3.0*epsilon_4);
for (unsigned int d = 0; d < dim; ++ d)
anisotropy[d]= as_n;
}
return std::make_pair(as_n, anisotropy);
}
//////////////////////////////////////////////////
#endif
// end of file
/* ---------------------------------------------------------------------
*
* Copyright (C) 2009 - 2022 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Authors: Katharina Kormann, Martin Kronbichler, Uppsala University,
* 2009-2012, updated to MPI version with parallel vectors in 2016
*/
/// 基于bt4修改，修改FEEaluation类，测试matrixfree 在左端项计算中考虑各向异性。
// First include the necessary files from the deal.II library.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_lib.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_matrix.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/error_estimator.h>
// This includes the data structures for the efficient implementation of
// matrix-free methods or more generic finite element operators with the class
// MatrixFree.
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/operators.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/distributed/solution_transfer.h>
// The following classes are used in parallel distributed computations and
// have all already been introduced in step-40:
#include <deal.II/base/index_set.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <iostream>
#include <fstream>
#include <time.h>
#include <memory>
#include "parameters.h"
#include "PFOperator.h"
#include "soluteOperator.h"
/*build phasefield constraints*/
#define DEAL_II_WITH_P4EST
namespace PhaseField
{
using namespace dealii;
/*initialzie phase field values*/
template <int dim>
class PFInitialValues : public Function<dim>
{
public:
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const override;
};
template <int dim>
double PFInitialValues<dim>::value(const Point<dim> &p,
const unsigned int component) const
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
double pf = -1.0;
double NP_pf = 0.0;
double d = 0, r0 = 2.1;
Point<dim> p0(length / 2.00, length / 2.0);
d = p.distance(p0);
pf = -tanh((d - r0) / sqrt(2.0));
if (pf > pf_upper)
{
NP_pf = NP_pf_upper;
}
else if (pf < pf_lower)
{
NP_pf = NP_pf_lower;
}
else
NP_pf = atanh(pf) * sqrt(2.0);
// if(d<r0)
// {NP_pf=r0-d;}
// if(d>=r0)
// {NP_pf=r0-d;}
return NP_pf;
}
/*initialize solute values*/
template <int dim>
class SoluteInitialValues : public Function<dim>
{
public:
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const override;
};
template <int dim>
double SoluteInitialValues<dim>::value(const Point<dim> &p,
const unsigned int component) const
{
// (void)component;
// Assert(component == 0, ExcIndexRange(component, 0, 1));
// double pf = -1.0;
// double d=0, r0=6.094;
// Point<dim> p0(length/2.00,length/2.0);
// d = p.distance(p0);
// pf = -tanh((d - r0)/sqrt(2.0));
// double u0 = 1.0-omega*(1.0-pc);
// double cl = u0*c_l0;
// double cs = pc*c_l0;
// double u = log(u0);//log((2.0*((1.0-pf)/2.0*cl+(1.0+pf)*cs/2.0)/c_l0)/(1.0+pc-(1.0-pc)*pf));
return 0.0;
}
////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
template <int dim>
class PFProblem
{
public:
PFProblem();
void run();
void setup_dof_system();
void setup_free_matrix();
void assemble_pf_rhs();
void assemble_solute_rhs();
void refine_mesh(const unsigned int min_grid_level,
const unsigned int max_grid_level);
// void numerical_truncation();
void solve_pf();
void solve_solute();
void output_results(const std::string &directory) const;
parallel::distributed::Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
MappingQ1<dim> mapping;
AffineConstraints<double> constraints;
// AffineConstraints<double> constraintsT;
using LevelMatrixType = PFOperator<dim, degree_finite_element, float>;
PFOperator<dim, degree_finite_element, double> pf_matrix;
soluteOperator<dim, degree_finite_element, double> solute_matrix;
MGConstrainedDoFs mg_constrained_dofs;
MGLevelObject<PFOperator<dim, degree_finite_element, float>> pf_mg_matrices;
MGLevelObject<soluteOperator<dim, degree_finite_element, float>> solute_mg_matrices;
MGTransferMatrixFree<dim, float> mg_transfer;
IndexSet locally_relevant_dofs;
LinearAlgebra::distributed::Vector<double> pf_rhs;
LinearAlgebra::distributed::Vector<double> solute_rhs;
LinearAlgebra::distributed::Vector<double> pf_solution;
LinearAlgebra::distributed::Vector<double> solute_solution;
LinearAlgebra::distributed::Vector<double> last_pf_solution;
LinearAlgebra::distributed::Vector<double> last_solute_solution;
MGLevelObject<LinearAlgebra::distributed::Vector<float>> pf_mg_solution;
MGLevelObject<LinearAlgebra::distributed::Vector<float>> solute_mg_solution;
typename MatrixFree<dim, double>::AdditionalData additional_data;
std::shared_ptr<MatrixFree<dim, double>> system_mf_storage;
std::shared_ptr<MatrixFree<dim, double>> system_mf_storage_solute;
typename MatrixFree<dim, float>::AdditionalData additional_data_mg;
double setup_time, cal_time, time_step;
unsigned int timestep_number;
bool setup_matrix_flag;
ConditionalOStream pcout;
ConditionalOStream time_details;
};
template <int dim>
PFProblem<dim>::PFProblem()
: triangulation(MPI_COMM_WORLD,
typename Triangulation<dim>::MeshSmoothing(
Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening),
// Triangulation<dim>::limit_level_difference_at_vertices,
typename parallel::distributed::Triangulation<
dim>::Settings(parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy)),
fe(degree_finite_element), dof_handler(triangulation), system_mf_storage(new MatrixFree<dim, double>()), system_mf_storage_solute(new MatrixFree<dim, double>()), setup_time(0.), cal_time(0.), timestep_number(0), time_step(_time_step), setup_matrix_flag(false), pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0), time_details(std::cout,
true &&
Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
{
}
template <int dim>
void PFProblem<dim>::setup_dof_system()
{
Timer time;
setup_time = 0;
dof_handler.distribute_dofs(fe);
dof_handler.distribute_mg_dofs();
pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
const std::set<types::boundary_id> dirichlet_boundary_ids = {0};
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, dirichlet_boundary_ids);
locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
constraints.close();
{
additional_data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none; // TasksParallelScheme::partition_color; //TasksParallelScheme::partition_partition;
additional_data.mapping_update_flags = (update_values | update_gradients | update_JxW_values | update_quadrature_points);
additional_data_mg.tasks_parallel_scheme = MatrixFree<dim, float>::AdditionalData::none; // TasksParallelScheme::partition_color;///TasksParallelScheme::partition_partition;
additional_data_mg.mapping_update_flags = (update_values | update_gradients | update_JxW_values | update_quadrature_points);
system_mf_storage = std::make_shared<MatrixFree<dim, double>>();
system_mf_storage->reinit(mapping,
dof_handler,
constraints,
QGauss<1>(fe.degree + 1),
additional_data);
/* system_mf_storage_solute= std::make_shared<MatrixFree<dim, double>>() ;
system_mf_storage_solute->reinit(mapping,
dof_handler,
constraints,
QGauss<1>(fe.degree + 1),
additional_data); */
mg_transfer.clear();
mg_transfer.initialize_constraints(mg_constrained_dofs);
mg_transfer.build(dof_handler);
}
///////////////////////////////
pcout << "initialize the mg_matrices..." << std::endl;
{
pf_mg_matrices.clear_elements();
solute_mg_matrices.clear_elements();
const unsigned int nlevels = triangulation.n_global_levels();
pf_mg_matrices.resize(0, nlevels - 1);
pf_mg_solution.resize(0, nlevels - 1);
solute_mg_matrices.resize(0, nlevels - 1);
solute_mg_solution.resize(0, nlevels - 1);
for (unsigned int level = 0; level < nlevels; ++level)
{
const IndexSet relevant_dofs = DoFTools::extract_locally_relevant_level_dofs(dof_handler, level);
AffineConstraints<double> level_constraints;
level_constraints.reinit(relevant_dofs);
level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level));
level_constraints.close();
additional_data_mg.mg_level = level;
auto mg_mf_storage_level = std::make_shared<MatrixFree<dim, float>>();
mg_mf_storage_level->reinit(mapping,
dof_handler,
level_constraints,
QGauss<1>(fe.degree + 1),
additional_data_mg);
pf_mg_matrices[level].initialize(mg_mf_storage_level,
mg_constrained_dofs,
level);
solute_mg_matrices[level].initialize(mg_mf_storage_level,
mg_constrained_dofs,
level);
}
}
pcout << "finishing set up dof system.." << std::endl;
setup_free_matrix();
pf_matrix.initialize_dof_vector(pf_solution);
pf_matrix.initialize_dof_vector(pf_rhs);
pf_matrix.initialize_dof_vector(last_pf_solution);
solute_matrix.initialize_dof_vector(solute_solution);
solute_matrix.initialize_dof_vector(last_solute_solution);
solute_matrix.initialize_dof_vector(solute_rhs);
//////////////////////////////
}
template <int dim>
void PFProblem<dim>::setup_free_matrix()
{
pcout << "steup free matrix ...." << std::endl;
Timer time;
setup_time = 0;
pf_matrix.clear();
solute_matrix.clear();
///////////////////////////////////////////
{
pf_matrix.initialize(system_mf_storage);
solute_matrix.initialize(system_mf_storage);
}
/////////////////////////////////////////////////
const unsigned int nlevels = triangulation.n_global_levels();
for (unsigned int level = 0; level < nlevels; ++level)
{
pf_mg_matrices[level].initialize_dof_vector(pf_mg_solution[level]);
solute_mg_matrices[level].initialize_dof_vector(solute_mg_solution[level]);
}
//////////////////////////////////////////
setup_matrix_flag = true;
setup_time += time.wall_time();
time_details << "Setup matrix-free (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << 's' << std::endl;
}
template <int dim>
void PFProblem<dim>::assemble_pf_rhs()
{
pcout << "assembling pf rhs..." << std::endl;
pf_rhs = 0;
FEEvaluation<dim, degree_finite_element> phi(*pf_matrix.get_matrix_free());
FEEvaluation<dim, degree_finite_element> phi_solute(*solute_matrix.get_matrix_free());
for (unsigned int cell = 0; cell < pf_matrix.get_matrix_free()->n_cell_batches(); ++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(last_pf_solution);
phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
phi_solute.reinit(cell);
phi_solute.read_dof_values_plain(last_solute_solution);
phi_solute.evaluate(EvaluationFlags::values);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
auto NP_pf = phi.get_value(q);
auto pf = NP_pf;
for (unsigned int v = 0; v < NP_pf.size(); ++v)
{
pf[v] = tanh(NP_pf[v] / sqrt(2.0));
if (pf[v] > pf_upper)
{
pf[v] = 1.0;
}
if (pf[v] < pf_lower)
{
pf[v] = -1.0;
}
}
auto c_val = phi_solute.get_value(q);
auto as_n = get_isotropy<dim, double, VectorizedArray<double>>(phi.get_gradient(q));
auto NP_pf_grad = phi.get_gradient(q);
// auto rhs_val= (tau0 - undercooling / ((1 - pc) * c_l0 * m_l))*as_n*pf+ time_step*(pf-lambda*(c_val + undercooling / ((1 - pc) * c_l0 * m_l))*(1.0-pf*pf))*(1.0-pf*pf);
double random = double(rand()) / double(RAND_MAX) * 2.0 - 1.0;
auto noise = (1.0 - undercooling / (c_l0 * m_l)) * as_n * amp_pf * random;
auto rhs_val = noise + (1.0 - undercooling / (c_l0 * m_l)) * as_n * NP_pf/time_step - as_n * sqrt(2.0) * pf * (NP_pf_grad * NP_pf_grad) + sqrt(2.0) * pf - lambda * (1.0 - pf * pf) * sqrt(2.0) * (c_val + undercooling / ((1.0 - pc) * c_l0 * m_l));
phi.submit_value(rhs_val, q);
}
phi.integrate_scatter(EvaluationFlags::values, pf_rhs);
}
pf_rhs.compress(VectorOperation::add);
pcout << "finishing assembling the pf_rhs" << std::endl;
// MPI_Barrier(MPI_COMM_WORLD) ;
}
template <int dim>
void PFProblem<dim>::assemble_solute_rhs()
{
pcout << "assembling solute rhs..." << std::endl;
solute_rhs = 0;
FEEvaluation<dim, degree_finite_element> phi(*pf_matrix.get_matrix_free());
FEEvaluation<dim, degree_finite_element> phi_last(*pf_matrix.get_matrix_free());
FEEvaluation<dim, degree_finite_element> phi_solute(*solute_matrix.get_matrix_free());
for (unsigned int cell = 0; cell < solute_matrix.get_matrix_free()->n_cell_batches(); ++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(pf_solution);
phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
phi_last.reinit(cell);
phi_last.read_dof_values_plain(last_pf_solution);
phi_last.evaluate(EvaluationFlags::values);
phi_solute.reinit(cell);
phi_solute.read_dof_values_plain(last_solute_solution);
phi_solute.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
auto NP_pf = phi.get_value(q);
auto pf = NP_pf;
auto last_NP_pf = phi_last.get_value(q);
auto c_val = phi_solute.get_value(q);
auto NP_pf_grad = phi.get_gradient(q);
// auto c_grad = phi_solute.get_gradient(q);
for (unsigned int v = 0; v < NP_pf.size(); ++v)
{
pf[v] = tanh(NP_pf[v] / sqrt(2.0));
if (pf[v] > pf_upper)
{
pf[v] = 1.0;
}
if (pf[v] < pf_lower)
{
pf[v] = -1.0;
}
}
auto square_grad = NP_pf_grad * NP_pf_grad;
auto coe_solute = 0.5 * ((1.0 + pc) - (1.0 - pc) * pf); // 计算0.5*[(1+k) - (1-k)*phi]
auto coe_pf = 0.5 * (1.0 + (1.0 - pc) * c_val); // 计算0.5*[1 + (1-k)*U]
auto inverse_grad = square_grad;
for (unsigned int v = 0; v < square_grad.size(); ++v)
{
if (square_grad[v] < 1.0e-14)
{
inverse_grad[v] = 0.0;
}
else
{
inverse_grad[v] = 1.0 / sqrt(square_grad[v]);
}
}
double random = double(rand()) / double(RAND_MAX) * 2.0 - 1.0;
auto noise = 0.5 * (1.0 - pf) * random * amp_c;
auto rhs_val = coe_solute * c_val/_time_step + coe_pf * (1.0 - pf * pf) / sqrt(2.0) * (NP_pf - last_NP_pf)/time_step + noise;
// 0.5*[(1+k) - (1-k)*phi]*U + 0.5*[1 + (1 - k)*U] * (phi- phi_1)
phi_solute.submit_value(rhs_val, q);
// -1 / sqrt(2) * 0.5*[1 + (1-k)*U] * (phi - phi_1) * grad_phi / |grad_phi|
phi_solute.submit_gradient(-0.5 * (1.0 - pf * pf) * coe_pf * (NP_pf - last_NP_pf) * inverse_grad * NP_pf_grad/_time_step, q);
}
phi_solute.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, solute_rhs);
// phi_solute.integrate_scatter(EvaluationFlags::values, solute_rhs);
}
solute_rhs.compress(VectorOperation::add);
}
template <int dim>
void PFProblem<dim>::refine_mesh(const unsigned int min_grid_level,
const unsigned int max_grid_level)
{
pcout << " ---------- Refining mesh" << std::endl;
Timer time;
setup_time += time.wall_time();
QGauss<dim> quadrature_formula(degree_finite_element + 1);
FEValues<dim> phasefield_fe_values(fe,
quadrature_formula,
update_values);
const unsigned int n_q_points = quadrature_formula.size();
std::vector<double> phasefield_values(n_q_points);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell != endc; ++cell)
{
if (cell->is_locally_owned())
{
phasefield_fe_values.reinit(cell);
phasefield_fe_values.get_function_values(pf_solution, phasefield_values);
unsigned int cell_level = cell->level();
bool flag_interface = false;
for (unsigned int q = 0; q < n_q_points; ++q)
{
double NP_pf = phasefield_values[q];
double pf = tanh(NP_pf / sqrt(2.0));
if ((!flag_interface) && (pf > -0.98) && (pf < 0.98))
{
flag_interface = true;
}
}
if (flag_interface)
{
if (cell_level < max_grid_level)
cell->set_refine_flag();
}
else if (cell_level > min_grid_level)
cell->set_coarsen_flag();
if (cell_level >= max_grid_level)
cell->clear_refine_flag();
if ((cell_level <= min_grid_level) && (cell->coarsen_flag_set()))
cell->clear_coarsen_flag();
}
}
pcout << " transferimg the data..." << std::endl;
parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>> solution_trans_pf(dof_handler);
parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>> solution_trans_solute(dof_handler);
triangulation.prepare_coarsening_and_refinement();
solution_trans_pf.prepare_for_coarsening_and_refinement(pf_solution);
solution_trans_solute.prepare_for_coarsening_and_refinement(solute_solution);
triangulation.execute_coarsening_and_refinement();
setup_dof_system();
MPI_Barrier(MPI_COMM_WORLD);
pf_solution.zero_out_ghost_values();
solute_solution.zero_out_ghost_values();
solution_trans_pf.interpolate(pf_solution);
solution_trans_solute.interpolate(solute_solution);
MPI_Barrier(MPI_COMM_WORLD);
pf_solution.update_ghost_values();
solute_solution.update_ghost_values();
pcout << " refining mesh is finished..." << std::endl;
time_details << "refining mesh time (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << "s\n";
}
// template<int dim>
// void PFProblem<dim>::numerical_truncation()
// {
// pcout << " ---------- numerical_truncation" << std::endl;
// Timer time;
// setup_time += time.wall_time();
// QGauss<dim> quadrature_formula(degree_finite_element +1);
// FEValues<dim> phasefield_fe_values (fe,
// quadrature_formula,
// update_values);
// const unsigned int n_q_points = quadrature_formula.size();
// std::vector<double> phasefield_values(n_q_points);
// typename DoFHandler<dim>::active_cell_iterator
// cell = dof_handler.begin_active(),
// endc = dof_handler.end();
// for (; cell!=endc; ++cell)
// {
// if (cell->is_locally_owned())
// {
// phasefield_fe_values.reinit (cell);
// phasefield_fe_values.get_function_values(pf_solution, phasefield_values);
// unsigned int cell_level = cell->level();
// for (unsigned int q = 0; q < n_q_points; ++q)
// {
// double pf = phasefield_values[q];
// if(pf>10.26)
// {
// cell->set_boundary_id(10);
// }
// if(pf<-10.26)
// {
// cell->set_boundary_id(11);
// }
// }
// }
// }
// }
template <int dim>
void PFProblem<dim>::solve_pf()
{
// pcout<<"====preparing solving the pf system...."<<std::endl;
Timer time;
last_pf_solution.update_ghost_values();
pf_matrix.evaluate_coefficient(last_pf_solution);
mg_transfer.interpolate_to_mg(dof_handler, pf_mg_solution, last_pf_solution);
using SmootherType =
PreconditionChebyshev<LevelMatrixType,
LinearAlgebra::distributed::Vector<float>>;
mg::SmootherRelaxation<SmootherType,
LinearAlgebra::distributed::Vector<float>>
mg_smoother;
MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
smoother_data.resize(0, triangulation.n_global_levels() - 1);
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level)
{
if (level > 0)
{
smoother_data[level].smoothing_range = 15.;
smoother_data[level].degree = 5;
smoother_data[level].eig_cg_n_iterations = 10;
}
else
{
smoother_data[0].smoothing_range = 1e-3;
smoother_data[0].degree = numbers::invalid_unsigned_int;
smoother_data[0].eig_cg_n_iterations = pf_mg_matrices[0].m();
}
pf_mg_matrices[level].evaluate_coefficient(pf_mg_solution[level]);
pf_mg_matrices[level].compute_diagonal();
smoother_data[level].preconditioner =
pf_mg_matrices[level].get_matrix_diagonal_inverse();
}
mg_smoother.initialize(pf_mg_matrices, smoother_data);
MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<float>> mg_coarse;
mg_coarse.initialize(mg_smoother);
mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(pf_mg_matrices);
MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<LevelMatrixType>> mg_interface_matrices;
mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1);
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level)
mg_interface_matrices[level].initialize(pf_mg_matrices[level]);
mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_interface(
mg_interface_matrices);
Multigrid<LinearAlgebra::distributed::Vector<float>> mg(
mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
mg.set_edge_matrices(mg_interface, mg_interface);
PreconditionMG<dim,
LinearAlgebra::distributed::Vector<float>,
MGTransferMatrixFree<dim, float>>
preconditioner(dof_handler, mg, mg_transfer);
SolverControl solver_control(pf_matrix.m(), 1e-8 * pf_rhs.l2_norm());
SolverCG<LinearAlgebra::distributed::Vector<double>> cg(solver_control);
setup_time += time.wall_time();
time_details << "MG build smoother time (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s\n";
pcout << "Total setup time (wall) " << setup_time << "s\n";
MPI_Barrier(MPI_COMM_WORLD);
time.reset();
time.start();
pf_solution = 0;
pf_solution.zero_out_ghost_values();
// PreconditionSSOR<SparseMatrix<double>> preconditioner2;
// preconditioner2.initialize(system_matrix, 1.2);
pcout << "solving the pf system....." << std::endl;
cg.solve(pf_matrix, pf_solution, pf_rhs, preconditioner);
constraints.distribute(pf_solution);
// constraints.distribute(pf_solution);
// locally_relevant_solution.copy_locally_owned_data_from(pf_solution);
// locally_relevant_solution.update_ghost_values();
double pf_val[2] = {std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()};
// pf_solution.zero_out_ghost_values();
// this_pf_solution.zero_out_ghost_values();
for (unsigned int i = 0;
i < pf_solution.size();
++i)
{
if (pf_solution.in_local_range(i))
{
double NP_pf = pf_solution[i];
pf_val[0] = std::min<double>(pf_val[0], NP_pf);
pf_val[1] = std::max<double>(pf_val[1], NP_pf);
if (NP_pf > NP_pf_upper)
{
NP_pf = NP_pf_upper;
}
else if (NP_pf < NP_pf_lower)
{
NP_pf = NP_pf_lower;
}
// if (pf > 0.99999) pf = 1.0;
// else if (pf< -0.99999) pf = -1.0;
pf_solution[i] = NP_pf;
// this_pf_solution[i] = pf;
}
}
pf_solution.update_ghost_values();
double global_pf[2];
pf_val[0] *= -1.0;
Utilities::MPI::max(pf_val, MPI_COMM_WORLD, global_pf);
global_pf[0] *= -1.0;
pcout << "Time solve (" << solver_control.last_step() << " iterations)"
<< (solver_control.last_step() < 10 ? " " : " ") << "(CPU/wall) "
<< time.cpu_time() << "s/" << time.wall_time() << "s\n";
pcout << "residul=" << solver_control.last_value() << std::endl
<< " --------------------------------------------------\n"
<< " Phase-field Range : "
<< global_pf[0] << ' ' << global_pf[1] << "\n"
<< " --------------------------------------------------"
<< std::endl;
}
template <int dim>
void PFProblem<dim>::solve_solute()
{
// pcout<<"====preparing solving the solute system...."<<std::endl;
Timer time;
pf_solution.update_ghost_values();
solute_matrix.evaluate_coefficient(pf_solution);
mg_transfer.interpolate_to_mg(dof_handler, solute_mg_solution, pf_solution);
time.restart();
using SmootherType =
PreconditionChebyshev<soluteOperator<dim, degree_finite_element, float>,
LinearAlgebra::distributed::Vector<float>>;
mg::SmootherRelaxation<SmootherType,
LinearAlgebra::distributed::Vector<float>>
mg_smoother;
MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
smoother_data.resize(0, triangulation.n_global_levels() - 1);
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level)
{
if (level > 0)
{
smoother_data[level].smoothing_range = 15.;
smoother_data[level].degree = 5;
smoother_data[level].eig_cg_n_iterations = 10;
}
else
{
smoother_data[0].smoothing_range = 1e-3;
smoother_data[0].degree = numbers::invalid_unsigned_int;
smoother_data[0].eig_cg_n_iterations = solute_mg_matrices[0].m();
}
// solute_mg_matrices[level].evaluate_coefficient_lhs(solute_mg_solution[level]);
solute_mg_matrices[level].evaluate_coefficient(solute_mg_solution[level]);
solute_mg_matrices[level].compute_diagonal();
smoother_data[level].preconditioner =
solute_mg_matrices[level].get_matrix_diagonal_inverse();
}
mg_smoother.initialize(solute_mg_matrices, smoother_data);
MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<float>>
mg_coarse;
mg_coarse.initialize(mg_smoother);
mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(
solute_mg_matrices);
MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<soluteOperator<dim, degree_finite_element, float>>>
mg_interface_matrices;
mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1);
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level)
mg_interface_matrices[level].initialize(solute_mg_matrices[level]);
mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_interface(
mg_interface_matrices);
Multigrid<LinearAlgebra::distributed::Vector<float>> mg(
mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
mg.set_edge_matrices(mg_interface, mg_interface);
PreconditionMG<dim,
LinearAlgebra::distributed::Vector<float>,
MGTransferMatrixFree<dim, float>>
preconditioner(dof_handler, mg, mg_transfer);
pcout << "solving the solute system....." << std::endl;
// SolverControl solver_control(600, 1e-12 * solute_rhs.l2_norm());
SolverControl solver_control(solute_matrix.m(), 1e-8 * solute_rhs.l2_norm());
SolverCG<LinearAlgebra::distributed::Vector<double>> cg(solver_control);
setup_time += time.wall_time();
time_details << "MG build smoother time (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s\n";
pcout << "Total setup time (wall) " << setup_time << "s\n";
time.reset();
time.start();
solute_solution = 0;
solute_solution.zero_out_ghost_values();
cg.solve(solute_matrix, solute_solution, solute_rhs, preconditioner);
constraints.distribute(solute_solution);
double solute_val[2] = {std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()};
for (unsigned int i = 0;
i < solute_solution.size();
++i)
{
if (solute_solution.in_local_range(i))
{
double solute = solute_solution[i];
solute_val[0] = std::min<double>(solute_val[0], solute);
solute_val[1] = std::max<double>(solute_val[1], solute);
// double pf = pf_solution[i];
// if (pf > 0.999999) solute_solution[i] =-0.001;
}
}
solute_solution.update_ghost_values();
double global_solute[2];
solute_val[0] *= -1.0;
Utilities::MPI::max(solute_val, MPI_COMM_WORLD, global_solute);
global_solute[0] *= -1.0;
pcout << "Time solve (" << solver_control.last_step() << " iterations)"
<< (solver_control.last_step() < 10 ? " " : " ") << "(CPU/wall) "
<< time.cpu_time() << "s/" << time.wall_time() << "s\n";
pcout << "residul=" << solver_control.last_value() << std::endl
<< " --------------------------------------------------\n"
<< " solute Range : "
<< global_solute[0] << ' ' << global_solute[1] << "\n"
<< " --------------------------------------------------"
<< std::endl;
}
template <int dim>
void PFProblem<dim>::output_results(const std::string &directory) const
{
pcout << "outputing data ....." << std::endl;
Timer time;
if (triangulation.n_global_active_cells() > 100000000)
return;
DataOut<dim> data_out;
LinearAlgebra::distributed::Vector<double> output_solute(solute_solution);
for (unsigned int i = 0; i < solute_solution.size(); ++i)
{
if (solute_solution.in_local_range(i))
{
double U_val = solute_solution[i];
double pf = pf_solution[i];
output_solute[i] = (((1 - pc) * U_val + 1) * (1 + pc - (1 - pc) * pf) / 2.0 * c_l0);
}
}
output_solute.update_ghost_values();
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(pf_solution, "pf");
data_out.add_data_vector(solute_solution, "solute");
// data_out.add_data_vector(output_solute, "solute");
data_out.build_patches(mapping);
DataOutBase::VtkFlags flags;
flags.compression_level = DataOutBase::VtkFlags::best_speed;
data_out.set_flags(flags);
data_out.write_vtu_with_pvtu_record(directory, "solution",
timestep_number, MPI_COMM_WORLD, 6);
time_details << "Time write output (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << "s, out dir is " << directory << "\n";
pf_solution.zero_out_ghost_values();
solute_solution.zero_out_ghost_values();
}
template <int dim>
void PFProblem<dim>::run()
{
{
const unsigned int n_vect_doubles = VectorizedArray<double>::size();
const unsigned int n_vect_bits = 8 * sizeof(double) * n_vect_doubles;
pcout << "VecSize= " << VectorizedArray<double>::size();
pcout << "Vectorization over " << n_vect_doubles
<< " doubles = " << n_vect_bits << " bits ("
<< Utilities::System::get_current_vectorization_level() << ')'
<< std::endl;
pcout <<"tau0=" <<tau0<<std::endl;
pcout <<"w0=" <<w0<<std::endl;
}
const unsigned int initial_global_refinement = 5;
const unsigned int maximum_level = 6, mini_level = 2;
const unsigned int n_adaptive_pre_refinement_steps = maximum_level - initial_global_refinement;
// const unsigned int adapt_interval =8;
double interval = 0.1 / time_step;
// const unsigned int save_interval =100*int(interval);
const unsigned int save_interval = 24;
const unsigned int adapt_interval = 100;
Point<dim> p1, p2;
for (int l = 0; l < dim; ++l)
{
p1[l] = 0.0;
p2[l] = length;
}
std::vector<unsigned int> rep(dim);
for (int k = 0; k < dim; k++)
rep[k] = p2[1] / 50.0;
GridGenerator::subdivided_hyper_rectangle(triangulation, rep, p1, p2);
triangulation.refine_global(initial_global_refinement);
setup_dof_system();
VectorTools::interpolate(dof_handler,
PFInitialValues<dim>(),
pf_solution);
pf_solution.update_ghost_values();
solute_solution.update_ghost_values();
// numerical_truncation();/*set pf_value>10.26 id = 10 pf_value<-10.26 id=11*/
// constraints.clear();
// locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
// DoFTools::make_hanging_node_constraints(dof_handler,constraints);
// VectorTools::interpolate_boundary_values(dof_handler,10,ConstantFunction<dim>(NP_pf_upper),constraints);
// VectorTools::interpolate_boundary_values(dof_handler,11,ConstantFunction<dim>(NP_pf_lower),constraints);
// constraints.close();
// pf_solution.update_ghost_values();
// solute_solution.update_ghost_values();
for (unsigned i = 0; i < initial_global_refinement; ++i)
{
refine_mesh(mini_level, maximum_level);
pf_solution = 0.0;
VectorTools::interpolate(dof_handler, PFInitialValues<dim>(), pf_solution);
pf_solution.update_ghost_values();
solute_solution.update_ghost_values();
// numerical_truncation();
// constraints.clear();
// locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
// DoFTools::make_hanging_node_constraints(dof_handler,constraints);
// VectorTools::interpolate_boundary_values(dof_handler,10,ConstantFunction<dim>(NP_pf_upper),constraints);
// VectorTools::interpolate_boundary_values(dof_handler,11,ConstantFunction<dim>(NP_pf_lower),constraints);
// constraints.close();
// pf_solution.update_ghost_values();
// solute_solution.update_ghost_values();
}
solute_solution = 0.0;
VectorTools::interpolate(dof_handler, SoluteInitialValues<dim>(), solute_solution);
solute_solution.update_ghost_values();
output_results("result1/");
time_t first, second;
first = time(NULL);
Timer time1;
while (undercooling <= 2.0)
{
cal_time += time_step;
++timestep_number;
undercooling = delta_T + cal_time * CR; // 实时的过冷度
pcout << "Time step " << timestep_number << " at t=" << cal_time
<< std::endl;
pcout << "undercooling = " << undercooling << "K" << std::endl;
last_pf_solution.zero_out_ghost_values();
last_solute_solution.zero_out_ghost_values();
for (unsigned int i = 0; i < pf_solution.size(); ++i)
{
if (pf_solution.in_local_range(i))
{
last_pf_solution[i] = pf_solution[i];
last_solute_solution[i] = solute_solution[i];
}
}
last_solute_solution.update_ghost_values();
last_pf_solution.update_ghost_values();
// numerical_truncation();
// constraints.clear();
// locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
// DoFTools::make_hanging_node_constraints(dof_handler,constraints);
// VectorTools::interpolate_boundary_values(dof_handler,10,ConstantFunction<dim>(NP_pf_upper),constraints);
// VectorTools::interpolate_boundary_values(dof_handler,11,ConstantFunction<dim>(NP_pf_lower),constraints);
// constraints.close();
// pf_solution.update_ghost_values();
// solute_solution.update_ghost_values();
if (!setup_matrix_flag)
{
setup_free_matrix();
}
assemble_pf_rhs();
solve_pf();
assemble_solute_rhs();
solve_solute();
setup_matrix_flag = false;
if (timestep_number % adapt_interval == 0)
{
refine_mesh(mini_level, maximum_level);
}
if (timestep_number % save_interval == 0)
{
output_results("result1/");
}
}
pcout << "The computing time is (CPU/wall) " << time1.cpu_time() << "s/" << time1.wall_time() << 's' << std::endl;
second = time(NULL);
pcout << "the used time is :" << difftime(second, first) << std::endl;
}
}
int main(int argc, char *argv[])
{
try
{
using namespace PhaseField;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
PFProblem<dimension> pf_problem;
pf_problem.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
////PFOperator.h
#ifndef PFOperator_h
#define PFOperator_h
#include"parameters.h"
using namespace dealii;
namespace PhaseField
{
template <int dim, int fe_degree, typename number>
class PFOperator
: public MatrixFreeOperators::
Base<dim, LinearAlgebra::distributed::Vector<number>>
{
public:
using value_type = number;
using FECellIntegrator =
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number>;
PFOperator();
void clear() override;
void evaluate_coefficient(
const LinearAlgebra::distributed::Vector<number> &solution);
virtual void compute_diagonal() override;
private:
virtual void apply_add(
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src) const override;
void
local_apply(const MatrixFree<dim, number> & data,
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;
void local_compute_diagonal(FECellIntegrator &integrator) const;
Table<2, VectorizedArray<number> > coefficient;
Table<2, std::pair<VectorizedArray<number>,
Tensor<1,dim, VectorizedArray<number> > > > anisotropy;
double kn;
};
template <int dim, int fe_degree, typename number>
PFOperator<dim, fe_degree, number>::PFOperator()
: MatrixFreeOperators::Base<dim,
LinearAlgebra::distributed::Vector<number>>(),
kn(_time_step)
{}
template <int dim, int fe_degree, typename number>
void PFOperator<dim, fe_degree, number>::clear()
{
coefficient.reinit(0, 1.0);
anisotropy.reinit(0, 0.0);
MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number>>::
clear();
}
template <int dim, int fe_degree, typename number>
void PFOperator<dim, fe_degree, number>::evaluate_coefficient(const LinearAlgebra::distributed::Vector<number> &solution)
{
const unsigned int n_cells = this->data->n_cell_batches();
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(*this->data);
anisotropy.reinit(n_cells, phi.n_q_points);
for (unsigned int cell = 0; cell < n_cells; ++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(solution);
phi.evaluate( EvaluationFlags::gradients);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
anisotropy(cell, q) =
get_anisotropy_vectorized<dim,number, VectorizedArray<number> >(phi.get_gradient(q));
}
}
}
template <int dim, int fe_degree, typename number>
void PFOperator<dim, fe_degree, number>::local_apply(
const MatrixFree<dim, number> & data,
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> & cell_range) const
{
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(data);
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
{
phi.reinit(cell);
phi.read_dof_values(src);
phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
auto bas_grad = phi.get_gradient(q);
for (unsigned int d =0; d<dim; ++d)
{
bas_grad[d] = anisotropy(cell,q).second[d]*bas_grad[d];
//anisotropy_grad[d] = anisotropy(cell, q).second[d];
}
//phi.submit_value(anisotropy(cell, q).first * phi.get_value(q), q);
phi.submit_value(anisotropy(cell, q).first * phi.get_value(q) * (1.0-undercooling/(c_l0*m_l))/kn, q);
phi.submit_gradient(bas_grad, q);
}
phi.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
phi.distribute_local_to_global(dst);
//phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
}
}
template <int dim, int fe_degree, typename number>
void PFOperator<dim, fe_degree, number>::apply_add(
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src) const
{
this->data->cell_loop(&PFOperator::local_apply, this, dst, src);
}
////////////////////////////////////////////////////////////////////
template <int dim, int fe_degree, typename number>
void PFOperator<dim, fe_degree, number>::local_compute_diagonal(
FECellIntegrator &phi) const
{
const unsigned int cell = phi.get_current_cell_index();
phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
auto bas_grad = phi.get_gradient(q);
// auto anisotropy_grad = phi.get_gradient(q);
for (unsigned int d =0; d<dim; ++d)
{
bas_grad[d] = anisotropy(cell, q).second[d]*bas_grad[d];
//anisotropy_grad[d] = anisotropy(cell, q).second[d] * anisotropy_grad[d];
}
phi.submit_value(anisotropy(cell, q).first * phi.get_value(q)* (tau0 - undercooling / ((1 - pc) * c_l0 * m_l))/kn, q);
phi.submit_gradient(bas_grad, q);
}
phi.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
}
template <int dim, int fe_degree, typename number>
void PFOperator<dim, fe_degree, number>::compute_diagonal()
{
this->inverse_diagonal_entries.reset(
new DiagonalMatrix<LinearAlgebra::distributed::Vector<number>>());
LinearAlgebra::distributed::Vector<number> &inverse_diagonal = this->inverse_diagonal_entries->get_vector();
this->data->initialize_dof_vector(inverse_diagonal);
MatrixFreeTools::compute_diagonal(*this->data,
inverse_diagonal,
&PFOperator::local_compute_diagonal,
this);
for (auto &diagonal_element : inverse_diagonal)
{
diagonal_element = (std::abs(diagonal_element) > 1.0e-10) ?
(1.0 / diagonal_element) :
1.0;
}
}
}
namespace ChangeVectorTypes
{
template <typename number>
void copy(TrilinosWrappers::MPI::Vector & out,
const LinearAlgebra::distributed::Vector<number> &in)
{
LinearAlgebra::ReadWriteVector<double> rwv(
out.locally_owned_elements());
in.update_ghost_values();
rwv.import(in, VectorOperation::insert);
out.import(rwv, VectorOperation::insert);
out.compress(VectorOperation::insert);
}
template <typename number>
void copy(LinearAlgebra::distributed::Vector<number> &out,
const TrilinosWrappers::MPI::Vector & in)
{
LinearAlgebra::ReadWriteVector<double> rwv;
rwv.reinit(in);
out.import(rwv, VectorOperation::insert);
// out.update_ghost_values();
}
} // namespace ChangeVectorTypes
#endif
////soluteOperator.h
#ifndef soluteOperator_h
#define soluteOperator_h
#include"parameters.h"
using namespace dealii;
namespace PhaseField
{
template <int dim, int fe_degree, typename number>
class soluteOperator:public MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number>>
{
public:
using value_type = number;
using FECellIntegrator = FEEvaluation<dim, fe_degree, fe_degree+1, 1, number>;
soluteOperator();
void clear() override;
virtual void compute_diagonal() override;
// void evaluate_coefficient_lhs(const LinearAlgebra::distributed::Vector<number> &phasefield_solution); //计算0.5 * [(1+k) -(1-k) * phi]
void evaluate_coefficient(const LinearAlgebra::distributed::Vector<number> &phasefield_solution);//计算0.5*(1 - pf);
// void evaluate_second(const LinearAlgebra::distributed::Vector<number> &phasefield_solution,
// const LinearAlgebra::distributed::Vector<number> &last_phasefield_solution);
private:
virtual void apply_add(
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src
)const override;
double kn;
void local_apply(const MatrixFree<dim, number> &data,
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;
void local_compute_diagonal(FECellIntegrator &integrator) const;
Table<2, std::pair<VectorizedArray<number>, VectorizedArray<number>>> coefficient;
//Table<2, VectorizedArray<number>> coefficient;
};
template<int dim, int fe_degree, typename number>
soluteOperator<dim, fe_degree, number>::soluteOperator()
:MatrixFreeOperators::Base<dim,
LinearAlgebra::distributed::Vector<number>>(),
kn(_time_step)//时间步长
{}
template<int dim, int fe_degree, typename number>
void soluteOperator<dim, fe_degree, number>::clear()
{
coefficient.reinit(0, 0.0);
MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number>>::clear();
}
template<int dim, int fe_degree, typename number>
void soluteOperator<dim, fe_degree, number>::evaluate_coefficient(const LinearAlgebra::distributed::Vector<number> &phasefield_solution)
{
const unsigned int n_cells = this->data->n_cell_batches();
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(*this->data);
coefficient.reinit(n_cells, phi.n_q_points);
for(unsigned int cell = 0; cell<n_cells;++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(phasefield_solution);
phi.evaluate(EvaluationFlags::values);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
auto NP_pf = phi.get_value(q);
auto pf = NP_pf;
for(unsigned int v=0;v<NP_pf.size();++v)
{
pf[v] = tanh(NP_pf[v]/sqrt(2.0));
if(pf[v]>pf_upper)
{pf[v]=1.0;}
if(pf[v]<pf_lower)
{pf[v]=-1.0;}
}
// for (unsigned int v =0; v < VectorizedArray<number>::size(); ++v)
// {
// //if(!(pf[v]<1.01 && pf[v]>=-1.01)) std::cout<<" wrong ....."<<pf[v]<<std::endl;
// if(pf[v] < pf_lower) pf[v] =-1.00;
// else if(pf[v] > pf_upper) pf[v] = 1.00;
// }
auto coe1 = 0.5 * ((1.0 + pc) - (1.0 - pc) * pf);
auto coe2 = 0.5 * D * (1.0 - pf) *kn;
coefficient(cell,q) = std::make_pair(coe1,coe2);
}
}
}
template <int dim, int fe_degree, typename number>
void soluteOperator<dim, fe_degree, number>::local_apply(
const MatrixFree<dim, number> &data,
LinearAlgebra::distributed::Vector<number> &dst,
const LinearAlgebra::distributed::Vector<number> &src,
const std::pair<unsigned int, unsigned int> &cell_range) const
{
FEEvaluation<dim, fe_degree, fe_degree+1, 1, number> phi(data);
for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
{
phi.reinit(cell);
phi.read_dof_values(src);
phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
for(unsigned int q=0;q<phi.n_q_points;q++)
{
phi.submit_value(coefficient(cell,q).first * phi.get_value(q)/kn, q);
phi.submit_gradient(coefficient(cell, q).second * phi.get_gradient(q)/kn, q);
//phi.submit_value(coefficient(cell,q) * phi.get_value(q), q);
}
phi.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
//phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
//phi.integrate_scatter(EvaluationFlags::values, dst);
phi.distribute_local_to_global(dst);
}
}
template <int dim, int fe_degree, typename number>
void soluteOperator<dim, fe_degree, number>::apply_add(
LinearAlgebra::distributed::Vector<number> & dst,
const LinearAlgebra::distributed::Vector<number> &src) const
{
this->data->cell_loop(&soluteOperator::local_apply, this, dst, src);
}
template<int dim, int fe_degree, typename number>
void soluteOperator<dim, fe_degree, number>::local_compute_diagonal(FECellIntegrator &phi) const
{
const unsigned int cell = phi.get_current_cell_index();
phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
// phi.evaluate(EvaluationFlags::values);
for(unsigned int q=0;q<phi.n_q_points;q++)
{
phi.submit_value(coefficient(cell,q).first * phi.get_value(q)/kn, q);
phi.submit_gradient(coefficient(cell,q).second * phi.get_gradient(q)/kn , q);
//phi.submit_value(coefficient(cell,q) * phi.get_value(q), q);
}
phi.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
//phi.integrate(EvaluationFlags::values);
}
template<int dim, int fe_degree, typename number>
void soluteOperator<dim, fe_degree, number>::compute_diagonal()
{
this->inverse_diagonal_entries.reset(
new DiagonalMatrix<LinearAlgebra::distributed::Vector<number>>());
LinearAlgebra::distributed::Vector<number> &inverse_diagonal =
this->inverse_diagonal_entries->get_vector();
this->data->initialize_dof_vector(inverse_diagonal);
MatrixFreeTools::compute_diagonal(*this->data,
inverse_diagonal,
&soluteOperator::local_compute_diagonal,
this);
for (auto &diagonal_element : inverse_diagonal)
{
diagonal_element = (std::abs(diagonal_element) > 1.0e-10) ?
(1.0 / diagonal_element) :
1.0;
}
}
}
#endif
