/*
make
mpirun --use-hwthread-cpus -np 6 main
mpirun --use-hwthread-cpus -np 6 -- main -NN 6 -dt 0
.1
*/
//MFEM
#include "mfem
.hpp"
#include <fstream>
#include <
iostream>
#include <algorithm>
#include <iomanip>
#include "Integ
.hpp"
using namespace std;
using namespace mfem;
#define mfemPrintf(
...) { \
if (Mpi::Root()) { \
printf(__VA_ARGS__); \
} \
}//令mfemPrintf函数可用
const real_t PARAM_NU = 1
.0; // viscosity
const real_t PARAM_KAPPA = 1
.0; // thermal diffusivity
const real_t PARAM_BETA = 1
.0; // thermal expansion coefficient
//exampleC
void uFun_ex(const Vector & x,real_t t, Vector & u)
{
real_t xi(x(0));
real_t yi(x(1));
real_t zi(x(2));
u(0) = yi*sin(t)*sin(t);
u(1) = zi*sin(t)*sin(t);
u(2) = xi*sin(t)*sin(t);
}
void gradu_exact(const Vector &xx,real_t t, Vector &u)
{
real_t x(xx(0));
real_t y(xx(1));
real_t z(xx(2));
// grad u1
u[0] = 0
.0;
u[1] = sin(t)*sin(t);
u[2] = 0
.0;
// grad u2
u[3] = 0
.0;
u[4] = 0
.0;
u[5] = sin(t)*sin(t);
// grad u3
u[6] = sin(t)*sin(t);
u[7] = 0
.0;
u[8] = 0
.0;
}
real_t pFun_ex(const Vector & x,real_t t)
{
real_t xi(x(0));
real_t yi(x(1));
real_t zi(x(2));
return (xi+yi+zi-1
.5)*sin(t)*sin(t);
}
void fFun(const Vector & x,real_t t, Vector & f)
{
real_t xi(x(0));
real_t yi(x(1));
real_t zi(x(2));
f(0) = zi*pow(sin(t),4)+pow(sin(t),2)+2*yi*cos(t)*sin(t);
f(1) = xi*pow(sin(t),4)+pow(sin(t),2)+2*zi*cos(t)*sin(t);
f(2) = yi*pow(sin(t),4)+pow(sin(t),2)+2*xi*cos(t)*sin(t);
}
// 温度场的解析解和源项
real_t thetaFun_ex(const Vector & x, real_t t)
{
real_t xi(x(0));
real_t yi(x(1));
real_t zi(x(2));
return (xi + yi + zi - 1
.5) * sin(t) * sin(t);
}
real_t psiFun(const Vector & x, real_t t)
{
real_t xi(x(0));
real_t yi(x(1));
real_t zi(x(2));
// 根据方程推导的源项
return (xi + yi + zi - 1
.5) * (2 * sin(t) * cos(t)) +
PARAM_KAPPA * 3 * sin(t) * sin(t) +
(yi * sin(t) * sin(t) * sin(t) * sin(t) +
zi * sin(t) * sin(t) * sin(t) * sin(t) +
xi * sin(t) * sin(t) * sin(t) * sin(t));
}
// 浮力项系数
class BuoyancyCoefficient : public VectorCoefficient
{
private:
Coefficient &theta_coeff;
Coefficient &beta_coeff;
Vector gravity;
public:
BuoyancyCoefficient(int dim, Coefficient &theta_, Coefficient &beta_)
: VectorCoefficient(dim), theta_coeff(theta_), beta_coeff(beta_)
{
gravity
.SetSize(dim);
gravity = 0
.0;
gravity(2) = -1
.0; // 重力方向为负z轴
}
virtual void Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip)
{
real_t theta_val = theta_coeff
.Eval(T, ip);
real_t beta_val = beta_coeff
.Eval(T, ip);
V
.Set(beta_val * theta_val, gravity);
}
};
//向量函数的H1误差
real_t ComputeVectorGradError(ParGridFunction &sol,
VectorCoefficient &exgrad,
const IntegrationRule *irs[])
{
ParFiniteElementSpace *fes = sol
.ParFESpace();
real_t error = 0
.0;
const FiniteElement *fe;
ElementTransformation *Tr;
DenseMatrix grad;
int intorder;
int sdim = fes->GetMesh()->SpaceDimension();
int vdim = fes->GetVDim();
Vector vec(vdim * sdim);
for (int i = 0; i < fes->GetNE(); i++)
{
fe = fes->GetFE(i);
Tr = fes->GetElementTransformation(i);
intorder = 2 * fe->GetOrder() + 3; // <--------
const IntegrationRule *ir;
if (irs)
{
ir = irs[fe->GetGeomType()];
}
else
{
ir = &(IntRules
.Get(fe->GetGeomType(), intorder));
}
for (int j = 0; j < ir->GetNPoints(); j++)
{
const IntegrationPoint &ip = ir->IntPoint(j);
Tr->SetIntPoint(&ip);
sol
.GetVectorGradient(*Tr, grad);
exgrad
.Eval(vec, *Tr, ip);
for (int d = 0; d < vdim; d++)
{
for (int k = 0; k < sdim; k++)
{
vec(d * sdim + k) -= grad(d, k);
}
}
error += ip
.weight * Tr->Weight() * (vec * vec);
}
}
real_t glb_error = 0
.0;
MPI_Allreduce(&error, &glb_error, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
fes->GetComm());
return sqrt(glb_error);
}
// Function to compute the integral of a ParGridFunction and subtract its mean value
void Zeromean(ParGridFunction &p)
{
ParFiniteElementSpace *p_space = p
.ParFESpace();
ParLinearForm int_form(p_space);
ConstantCoefficient one_coeff(1
.0);
int_form
.AddDomainIntegrator(new DomainLFIntegrator(one_coeff));
int_form
.Assemble();
real_t int_p = 0
.0;
real_t int_p_loc = InnerProduct(int_form, p);
//使用 MPI 的 MPI_Allreduce 函数将所有进程中的局部积分结果 int_p_loc 进行全局归约（求和），并将结果存储在 int_p 中
MPI_Allreduce(&int_p_loc, &int_p, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
//mfemPrintf("Integral of p: %lg, ", int_p);
ParGridFunction one_gf(p_space);
one_gf
.ProjectCoefficient(one_coeff);
real_t vol_loc = InnerProduct(int_form, one_gf);
real_t vol = 0
.0;
MPI_Allreduce(&vol_loc, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
real_t mean_p = int_p / vol;
p
.Add(-mean_p, one_gf);
int_p_loc = InnerProduct(int_form, p);
MPI_Allreduce(&int_p_loc, &int_p, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
//mfemPrintf("after modification: %lg\n", int_p);
}
// 改进的预条件器类，处理速度-压力-温度三场耦合
class BoussinesqPC : public mfem::Solver
{
private:
ParFiniteElementSpace *V_space;
ParFiniteElementSpace *P_space;
ParFiniteElementSpace *T_space;
BlockOperator &block_op;
Array<int> &block_trueOffsets;
// 速度块预条件器组件
HypreParMatrix *pc_Fmat; // 修改名称，添加pc_前缀
HypreBoomerAMG *
prec_F;
CGSolver *solver_F;
// 压力Schur补预条件器组件
HypreParMatrix *pc_Kpmat; // 修改名称，添加pc_前缀
HypreBoomerAMG *
prec_Kp;
CGSolver *solver_Kp;
HypreParMatrix *pc_Mpmat; // 修改名称，添加pc_前缀
HypreDiagScale *
prec_Mp;
CGSolver *solver_Mp;
// 温度Schur补预条件器组件
HypreParMatrix *pc_Ktmat; // 修改名称，添加pc_前缀
HypreBoomerAMG *
prec_Kt;
CGSolver *solver_Kt;
HypreParMatrix *pc_Mtmat; // 修改名称，添加pc_前缀
HypreDiagScale *
prec_Mt;
CGSolver *solver_Mt;
// 耦合项矩阵
HypreParMatrix *pc_Gmat; // 修改名称，添加pc_前缀
real_t dt_, nu_, kappa_, beta_;
public:
BoussinesqPC(ParFiniteElementSpace *V_space_,
ParFiniteElementSpace *P_space_,
ParFiniteElementSpace *T_space_,
BlockOperator &block_op_,
Array<int> &offsets_,
real_t dt_, real_t nu_, real_t kappa_, real_t beta_)
: Solver(offsets_[3]),
V_space(V_space_), P_space(P_space_), T_space(T_space_),
block_op(block_op_), block_trueOffsets(offsets_),
dt_(dt_), nu_(nu_), kappa_(kappa_), beta_(beta_),
pc_Fmat(nullptr), pc_Kpmat(nullptr), pc_Mpmat(nullptr),
pc_Ktmat(nullptr), pc_Mtmat(nullptr), pc_Gmat(nullptr)
{
// 初始化速度块预条件器
ConstantCoefficient oneondt(1
.0 / dt_);
ConstantCoefficient nu_coeff(nu_);
ParBilinearForm F_varf(V_space);
F_varf
.AddDomainIntegrator(new VectorMassIntegrator(oneondt));
F_varf
.AddDomainIntegrator(new VectorDiffusionIntegrator(nu_coeff));
F_varf
.Assemble();
F_varf
.Finalize();
pc_Fmat = F_varf
.ParallelAssemble();
if (pc_Fmat == nullptr) {
if (Mpi::Root()) {
std::cerr << "错误: pc_Fmat 矩阵装配失败!" << std::endl;
}
MPI_Abort(MPI_COMM_WORLD, 1);
}
prec_F = new HypreBoomerAMG(*pc_Fmat);
prec_F->SetPrintLevel(0);
solver_F = new CGSolver(MPI_COMM_WORLD);
solver_F->Set
Preconditioner(*
prec_F);
solver_F->SetOperator(*pc_Fmat);
solver_F->SetRelTol(1e-5);
solver_F->SetMaxIter(100);
solver_F->SetPrintLevel(0);
// 初始化压力Schur补预条件器
ParBilinearForm Mp_varf(P_space);
Mp_varf
.AddDomainIntegrator(new MassIntegrator);
Mp_varf
.Assemble();
Mp_varf
.Finalize();
pc_Mpmat = Mp_varf
.ParallelAssemble();
if (pc_Mpmat == nullptr) {
if (Mpi::Root()) {
std::cerr << "错误: pc_Mpmat 矩阵装配失败!" << std::endl;
}
MPI_Abort(MPI_COMM_WORLD, 1);
}
prec_Mp = new HypreDiagScale(*pc_Mpmat);
solver_Mp = new CGSolver(MPI_COMM_WORLD);
solver_Mp->Set
Preconditioner(*
prec_Mp);
solver_Mp->SetOperator(*pc_Mpmat);
solver_Mp->SetRelTol(1e-5);
solver_Mp->SetMaxIter(100);
solver_Mp->SetPrintLevel(0);
ParBilinearForm Kp_varf(P_space);
Kp_varf
.AddDomainIntegrator(new DiffusionIntegrator);
Kp_varf
.Assemble();
Array<int> ess_vdofs_p;
if(Mpi::Root())
{
ess_vdofs_p
.Append(0);
}
Kp_varf
.EliminateVDofs(ess_vdofs_p);
Kp_varf
.Finalize();
pc_Kpmat = Kp_varf
.ParallelAssemble();
if (pc_Kpmat == nullptr) {
if (Mpi::Root()) {
std::cerr << "错误: pc_Kpmat 矩阵装配失败!" << std::endl;
}
MPI_Abort(MPI_COMM_WORLD, 1);
}
prec_Kp = new HypreBoomerAMG(*pc_Kpmat);
prec_Kp->SetPrintLevel(0);
solver_Kp = new CGSolver(MPI_COMM_WORLD);
solver_Kp->Set
Preconditioner(*
prec_Kp);
solver_Kp->SetOperator(*pc_Kpmat);
solver_Kp->SetRelTol(1e-5);
solver_Kp->SetMaxIter(100);
solver_Kp->SetPrintLevel(0);
// 初始化温度Schur补预条件器
ParBilinearForm Mt_varf(T_space);
Mt_varf
.AddDomainIntegrator(new MassIntegrator);
Mt_varf
.Assemble();
Mt_varf
.Finalize();
pc_Mtmat = Mt_varf
.ParallelAssemble();
if (pc_Mtmat == nullptr) {
if (Mpi::Root()) {
std::cerr << "错误: pc_Mtmat 矩阵装配失败!" << std::endl;
}
MPI_Abort(MPI_COMM_WORLD, 1);
}
prec_Mt = new HypreDiagScale(*pc_Mtmat);
solver_Mt = new CGSolver(MPI_COMM_WORLD);
solver_Mt->Set
Preconditioner(*
prec_Mt);
solver_Mt->SetOperator(*pc_Mtmat);
solver_Mt->SetRelTol(1e-5);
solver_Mt->SetMaxIter(100);
solver_Mt->SetPrintLevel(0);
ParBilinearForm Kt_varf(T_space);
Kt_varf
.AddDomainIntegrator(new DiffusionIntegrator);
Kt_varf
.Assemble();
Array<int> ess_vdofs_t;
if(Mpi::Root())
{
ess_vdofs_t
.Append(0);
}
Kt_varf
.EliminateVDofs(ess_vdofs_t);
Kt_varf
.Finalize();
pc_Ktmat = Kt_varf
.ParallelAssemble();
if (pc_Ktmat == nullptr) {
if (Mpi::Root()) {
std::cerr << "错误: pc_Ktmat 矩阵装配失败!" << std::endl;
}
MPI_Abort(MPI_COMM_WORLD, 1);
}
prec_Kt = new HypreBoomerAMG(*pc_Ktmat);
prec_Kt->SetPrintLevel(0);
solver_Kt = new CGSolver(MPI_COMM_WORLD);
solver_Kt->Set
Preconditioner(*
prec_Kt);
solver_Kt->SetOperator(*pc_Ktmat);
solver_Kt->SetRelTol(1e-5);
solver_Kt->SetMaxIter(100);
solver_Kt->SetPrintLevel(0);
// 初始化浮力耦合项矩阵
Vector gravity(3);
gravity = 0
.0;
gravity(2) = -1
.0; // 重力方向为负z轴
VectorConstantCoefficient g_coeff(gravity);
ParMixedBilinearForm G_varf(T_space, V_space);
G_varf
.AddDomainIntegrator(new VectorMassIntegrator(g_coeff));
G_varf
.Assemble();
G_varf
.Finalize();
pc_Gmat = G_varf
.ParallelAssemble();
if (pc_Gmat == nullptr) {
if (Mpi::Root()) {
std::cerr << "错误: pc_Gmat 矩阵装配失败!" << std::endl;
}
MPI_Abort(MPI_COMM_WORLD, 1);
}
};
~BoussinesqPC() {
delete pc_Fmat;
delete pc_Kpmat;
delete pc_Mpmat;
delete pc_Ktmat;
delete pc_Mtmat;
delete pc_Gmat;
delete
prec_F;
delete solver_F;
delete
prec_Kp;
delete solver_Kp;
delete
prec_Mp;
delete solver_Mp;
delete
prec_Kt;
delete solver_Kt;
delete
prec_Mt;
delete solver_Mt;
}
void Mult(const mfem::Vector &x, mfem::Vector &y) const
{
StopWatch timer;
Vector v_in(x
.GetData() + block_trueOffsets[0],
block_trueOffsets[1] - block_trueOffsets[0]);
Vector p_in(x
.GetData() + block_trueOffsets[1],
block_trueOffsets[2] - block_trueOffsets[1]);
Vector t_in(x
.GetData() + block_trueOffsets[2],
block_trueOffsets[3] - block_trueOffsets[2]);
Vector v_out(y
.GetData() + block_trueOffsets[0],
block_trueOffsets[1] - block_trueOffsets[0]);
Vector p_out(y
.GetData() + block_trueOffsets[1],
block_trueOffsets[2] - block_trueOffsets[1]);
Vector t_out(y
.GetData() + block_trueOffsets[2],
block_trueOffsets[3] - block_trueOffsets[2]);
// 第一步：求解温度Schur补系统
Vector t_temp(t_in
.Size());
timer
.Clear();
timer
.Start();
solver_Kt->Mult(t_in, t_temp);
pc_Mtmat->Mult(t_temp, t_temp);
solver_Mt->Mult(t_temp, t_out);
timer
.Stop();
// 第二步：求解压力Schur补系统
Vector p_temp(p_in
.Size());
timer
.Clear();
timer
.Start();
solver_Kp->Mult(p_in, p_temp);
pc_Mpmat->Mult(p_temp, p_temp);
solver_Mp->Mult(p_temp, p_out);
timer
.Stop();
// 第三步：求解速度系统，考虑浮力耦合
Vector v_temp(v_in
.Size());
// 添加浮力项贡献
Vector Gt_result(v_temp
.Size());
if (pc_Gmat != nullptr) {
pc_Gmat->Mult(t_out, Gt_result);
v_temp
.Add(-beta_, Gt_result);
} else {
if (Mpi::Root()) {
std::cerr << "警告: pc_Gmat 为空，跳过浮力项" << std::endl;
}
}
// 添加压力项贡献 - 修复这里
Operator& Bt_op = block_op
.GetBlock(0, 1);
Vector Bp_result(v_temp
.Size());
Bt_op
.Mult(p_out, Bp_result); // 直接使用运算符，不进行dynamic_cast
v_temp
.Add(-1
.0, Bp_result);
timer
.Clear();
timer
.Start();
solver_F->Mult(v_temp, v_out);
timer
.Stop();
};
void SetOperator(const Operator &op) {};
};
int main(int argc, char *argv[])
{
StopWatch chrono;
// 1
. Initialize MPI and HYPRE
.
Mpi::Init(argc, argv);
int myid = Mpi::WorldRank();
Hypre::Init();
bool verbose = (myid == 0);
// 2
. Parse command-line options
.
const char *mesh_file = "
../data/star
.mesh";
int ref_levels = -1;
int order = 1;
real_t dt = 0
.05;
real_t t_final = 0
.05;
int NN = 4;//空间每个坐标轴划分为NN等分
/*
make NS_1
mpirun --use-hwthread-cpus -np 6 NS_1
*/
bool par_format = false;
bool pa = false;
const char *device_config = "cpu";
bool visualization = 1;
bool ad
ios2 = false;
OptionsParser args(argc, argv);
args
.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use
.");
args
.AddOption(&ref_levels, "-r", "--refine",
"Number of times to refine the mesh uniformly
.");
args
.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree)
.");
args
.AddOption(&par_format, "-pf", "--parallel-format", "-sf",
"--serial-format",
"Format to use when saving the results for VisIt
.");
args
.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
"--no-partial-assembly", "Enable Partial Assembly
.");
args
.AddOption(&device_config, "-d", "--device",
"Device configuration string, see Device::Configure()
.");
args
.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization
.");
args
.AddOption(&ad
ios2, "-ad
ios2", "--ad
ios2-streams", "-no-ad
ios2",
"--no-ad
ios2-streams",
"Save data using ad
ios2 streams
.");
args
.AddOption(&dt, "-dt", "--dt",
"dt
.");
args
.AddOption(&NN, "-NN", "--NN",
"NN
.");
args
.AddOption(&t_final, "-t_final", "--t_final",
"t_final
.");
args
.Parse();
if (!args
.Good())
{
if (verbose)
{
args
.PrintUsage(cout);
}
return 1;
}
if (verbose)
{
args
.PrintOptions(cout);
}
int nsteps = ceil(t_final / dt);
dt = t_final / nsteps;
// 3
. Enable hardware devices such as GPUs, and programming models such as
// CUDA, OCCA, RAJA and OpenMP based on command line options
.
Device device(device_config);
if (myid == 0) { device
.Print(); }
//Mesh
int Nx = NN, Ny = NN, Nz = NN;
real_t Sx = 1
.0, Sy = 1
.0, Sz = 1
.0;
Mesh mesha = Mesh::MakeCartesian3D(Nx, Ny, Nz, Element::TETRAHEDRON, Sx, Sy, Sz);
Mesh *mesh=&mesha;
int dim = mesh->Dimension();
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
pmesh->PrintInfo(std::cout);
//delete mesh;
//FiniteElementSpace
FiniteElementCollection *P2(new H1_FECollection(order+1, dim));
FiniteElementCollection *P1(new H1_FECollection(order));
ParFiniteElementSpace *V_space = new ParFiniteElementSpace(pmesh, P2, dim);
ParFiniteElementSpace *P_space = new ParFiniteElementSpace(pmesh, P1);
ParFiniteElementSpace *T_space = new ParFiniteElementSpace(pmesh, P1); // 温度空间
HYPRE_BigInt dimR = V_space->GlobalTrueVSize();
HYPRE_BigInt dimW = P_space->GlobalTrueVSize();//查看全局的真实向量自由度（True Vector DoFs）总数
HYPRE_BigInt dimT = T_space->GlobalTrueVSize();
if (verbose)//若当前进程是主进程，查看维数
{
std::cout << "***********************************************************\n";
std::cout << "dim(R) = " << dimR << "\n";
std::cout << "dim(W) = " << dimW << "\n";
std::cout << "dim(T) = " << dimT << "\n";
std::cout << "dim(R+W+T) = " << dimR + dimW + dimT << "\n";
std::cout << "***********************************************************\n";
}
//Define the coefficients, analytical solution, and rhs of the PDE
.
ConstantCoefficient k(1
.0);
ConstantCoefficient nu(PARAM_NU);
ConstantCoefficient kappa(PARAM_KAPPA);
ConstantCoefficient beta_coeff(PARAM_BETA);
ConstantCoefficient m_one(-1
.0);
ConstantCoefficient oneondt1(1
.0 / dt);
VectorFunctionCoefficient fcoeff(dim, fFun);
VectorFunctionCoefficient ucoeff(dim, uFun_ex);
FunctionCoefficient pcoeff(pFun_ex);
FunctionCoefficient thetacoeff(thetaFun_ex);
FunctionCoefficient psicoeff(psiFun);
Array<int> block_Offsets(4);//分块方法
block_Offsets[0] = 0;
block_Offsets[1] = V_space->GetTrueVSize();
block_Offsets[2] = P_space->GetTrueVSize();
block_Offsets[3] = T_space->GetTrueVSize();
block_Offsets
.PartialSum();
//存储解向量和右端向量，每个时间步都更新
Vector rhs(block_Offsets[3]);
Vector sol(block_Offsets[3]);
Vector rhs_u(rhs
.GetData(), block_Offsets[1] - block_Offsets[0]);//使用rhs的内存，不分配新的内存
Vector rhs_p(rhs
.GetData() + block_Offsets[1], block_Offsets[2] - block_Offsets[1]);
Vector rhs_t(rhs
.GetData() + block_Offsets[2], block_Offsets[3] - block_Offsets[2]);
Vector sol_u(sol
.GetData(), block_Offsets[1] - block_Offsets[0]);
Vector sol_p(sol
.GetData() + block_Offsets[1], block_Offsets[2] - block_Offsets[1]);
Vector sol_t(sol
.GetData() + block_Offsets[2], block_Offsets[3] - block_Offsets[2]);
//使用低阶格式得到前面几步的值
ParGridFunction uh_0(V_space),ph(P_space),th(T_space);
ParGridFunction beta_kn(V_space),gamma_kn(V_space);
ucoeff
.SetTime(0
.0);
pcoeff
.SetTime(0
.0);
thetacoeff
.SetTime(0
.0);
real_t t = 0
.0;
uh_0
.ProjectCoefficient(ucoeff);//设置初值
th
.ProjectCoefficient(thetacoeff);
chrono
.Clear();
chrono
.Start();
real_t errp=0
.0;
real_t errt=0
.0;
Array<int> ess_bdr(mesh->bdr_attributes
.Max());
ess_bdr = 1; // Dirichlet boundary on all boundary faces
ParBilinearForm F_old_varf(V_space);
F_old_varf
.AddDomainIntegrator(new VectorMassIntegrator(oneondt1));
F_old_varf
.Assemble();
// 在时间循环外声明矩阵指针，使用不同的名称
HypreParMatrix *main_Fmat = nullptr;
HypreParMatrix *main_Bmat = nullptr;
HypreParMatrix *main_Btmat = nullptr;
HypreParMatrix *main_Cmat = nullptr;
HypreParMatrix *main_Gmat = nullptr;
HypreParMatrix *main_Gtmat = nullptr;
ParGridFunction uh_n(V_space),uh(V_space),theta_n(T_space),theta(T_space);
uh_n=uh_0;
theta_n=th;
if (verbose) {
std::cout << "开始时间循环
..." << std::endl;
}
for (int n = 0; n < nsteps; n++)
{
if (verbose && myid == 0) {
std::cout << "时间步 " << n + 1 << "/" << nsteps << std::endl;
}
// 删除已存在的矩阵
if (main_Fmat) { delete main_Fmat; main_Fmat = nullptr; }
if (main_Bmat) { delete main_Bmat; main_Bmat = nullptr; }
if (main_Btmat) { delete main_Btmat; main_Btmat = nullptr; }
if (main_Cmat) { delete main_Cmat; main_Cmat = nullptr; }
if (main_Gmat) { delete main_Gmat; main_Gmat = nullptr; }
if (main_Gtmat) { delete main_Gtmat; main_Gtmat = nullptr; }
if(n>0){
uh_n=uh;
theta_n=theta;
}
ucoeff
.SetTime(t+dt);
pcoeff
.SetTime(t+dt);
thetacoeff
.SetTime(t+dt);
uh=0
.0;
theta=0
.0;
uh
.ProjectBdrCoefficient(ucoeff, ess_bdr);
theta
.ProjectBdrCoefficient(thetacoeff, ess_bdr);
beta_kn=uh_n;
gamma_kn=uh_n;
fcoeff
.SetTime(t+dt);
ParLinearForm f_form(V_space);
f_form
.AddDomainIntegrator(new VectorDomainLFIntegrator(fcoeff));
f_form
.AddDomainIntegrator(new VectorConvectionLFIntegrator(gamma_kn,-1
.0));
f_form
.Assemble();
F_old_varf
.AddMult(beta_kn, f_form, 1
.0);//将F_old*beta_kn(矩阵乘向量)加到f_form上
// 添加浮力项
BuoyancyCoefficient buoyancy_coeff(dim, thetacoeff, beta_coeff);
ParLinearForm buoyancy_form(V_space);
buoyancy_form
.AddDomainIntegrator(new VectorDomainLFIntegrator(buoyancy_coeff));
buoyancy_form
.Assemble();
f_form += buoyancy_form;
ParLinearForm g_form(P_space);
g_form
.Assemble();
psicoeff
.SetTime(t+dt);
ParLinearForm psi_form(T_space);
psi_form
.AddDomainIntegrator(new DomainLFIntegrator(psicoeff));
psi_form
.Assemble();
//Bilinear forms
ParBilinearForm F_varf(V_space);
F_varf
.AddDomainIntegrator(new VectorDiffusionIntegrator(nu));
F_varf
.AddDomainIntegrator(new VectorMassIntegrator(oneondt1));
F_varf
.Assemble();
ParMixedBilinearForm B_varf(V_space,P_space);
B_varf
.AddDomainIntegrator(new VectorDivergenceIntegrator(m_one));
B_varf
.Assemble();
ParBilinearForm C_varf(T_space);
C_varf
.AddDomainIntegrator(new DiffusionIntegrator(kappa));
C_varf
.AddDomainIntegrator(new MassIntegrator(oneondt1));
C_varf
.Assemble();
Vector gravity(3);
gravity = 0
.0;
gravity(2) = -1
.0;
VectorConstantCoefficient g_coeff(gravity);
ParMixedBilinearForm G_varf(T_space,V_space);
G_varf
.AddDomainIntegrator(new VectorMassIntegrator(g_coeff));
G_varf
.Assemble();
// Dirichlet boundary
F_varf
.EliminateEssentialBC(ess_bdr, uh, f_form);
B_varf
.EliminateTrialEssentialBC(ess_bdr, uh, g_form);
C_varf
.EliminateEssentialBC(ess_bdr, theta, psi_form);
F_varf
.Finalize();
B_varf
.Finalize();
C_varf
.Finalize();
G_varf
.Finalize();
// 创建新矩阵
main_Fmat = F_varf
.ParallelAssemble();
main_Bmat = B_varf
.ParallelAssemble();
main_Btmat = main_Bmat->Transpose();
main_Cmat = C_varf
.ParallelAssemble();
main_Gmat = G_varf
.ParallelAssemble();
main_Gtmat = main_Gmat->Transpose();
// 检查矩阵是否成功创建
if (main_Fmat == nullptr || main_Bmat == nullptr || main_Cmat == nullptr || main_Gmat == nullptr ||
main_Btmat == nullptr || main_Gtmat == nullptr) {
if (myid == 0) {
std::cerr << "错误: 矩阵装配失败!" << std::endl;
if (main_Fmat == nullptr) std::cerr << "main_Fmat is null" << std::endl;
if (main_Bmat == nullptr) std::cerr << "main_Bmat is null" << std::endl;
if (main_Btmat == nullptr) std::cerr << "main_Btmat is null" << std::endl;
if (main_Cmat == nullptr) std::cerr << "main_Cmat is null" << std::endl;
if (main_Gmat == nullptr) std::cerr << "main_Gmat is null" << std::endl;
if (main_Gtmat == nullptr) std::cerr << "main_Gtmat is null" << std::endl;
}
MPI_Abort(MPI_COMM_WORLD, 1);
}
{
BlockOperator block_op(block_Offsets);
block_op
.SetBlock(0, 0, main_Fmat);
block_op
.SetBlock(0, 1, main_Btmat);
block_op
.SetBlock(0, 2, main_Gtmat);
block_op
.SetBlock(1, 0, main_Bmat);
block_op
.SetBlock(2, 0, main_Gmat);
block_op
.SetBlock(2, 2, main_Cmat);
sol = 0
.0;rhs=0
.0;
f_form
.ParallelAssemble(rhs_u);
g_form
.ParallelAssemble(rhs_p);
psi_form
.ParallelAssemble(rhs_t);//组装右端向量
// compatibility condition
Vector one_vec(rhs_p
.Size());
one_vec = 1
.0;
real_t rhs_p_sum_loc = InnerProduct(rhs_p, one_vec);
real_t rhs_p_sum;
MPI_Allreduce(&rhs_p_sum_loc, &rhs_p_sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
rhs_p -= (rhs_p_sum) / P_space->GlobalTrueVSize();
BoussinesqPC
prec(V_space, P_space, T_space, block_op, block_Offsets, dt, PARAM_NU, PARAM_KAPPA, PARAM_BETA);
FGMRESSolver solver(MPI_COMM_WORLD);
solver
.SetRelTol(1e-12);
solver
.SetMaxIter(500);
solver
.SetOperator(block_op);
solver
.Set
Preconditioner(
prec);
solver
.SetPrintLevel(1);
solver
.Mult(rhs, sol);
uh
.Distribute(sol_u);
ph
.Distribute(sol_p);
theta
.Distribute(sol_t);//提取解
if(myid == 0){
cout <<"t = "<< t+dt << endl;
}
t += dt;
}
}
// 在程序最后统一清理矩阵
delete main_Fmat;
delete main_Bmat;
delete main_Btmat;
delete main_Cmat;
delete main_Gmat;
delete main_Gtmat;
Zeromean(ph);
errp += pow(ph
.ComputeL2Error(pcoeff),2);
errt += pow(theta
.ComputeL2Error(thetacoeff),2);
//计算误差
//H1误差
int gradu_dim = dim*dim;
VectorFunctionCoefficient gradu_exact_coeff(gradu_dim, gradu_exact);
gradu_exact_coeff
.SetTime(t);
real_t GradError = ComputeVectorGradError(uh, gradu_exact_coeff,nullptr);
//计算L2误差
real_t err_u;
err_u = uh
.ComputeL2Error(ucoeff);
errp=sqrt(errp*dt);
real_t err_p=ph
.ComputeL2Error(pcoeff);
errt=sqrt(errt*dt);
real_t err_theta=theta
.ComputeL2Error(thetacoeff);
int order_quad = max(2, 2*order+1);
const IntegrationRule *irs[Geometry::NumGeom];
for (int i=0; i < Geometry::NumGeom; ++i)
{
irs[i] = &(IntRules
.Get(i, order_quad));
}
real_t norm_u = ComputeGlobalLpNorm(2, ucoeff, *pmesh, irs);
chrono
.Stop();
if (myid == 0)
{
cout <<"t = "<< t << " || u_h - u_ex || = " << err_u << endl;
cout <<"t = "<< t << " err_p_L2 = " << err_p << endl;
cout <<"t = "<< t << " err_theta_L2 = " << err_theta << endl;
cout <<"t= "<< t << " H1_error = "<< GradError <<endl;
cout << "t = "<< t <<" norm_u = "<< norm_u << "\n";
std::cout << "It took " << chrono
.RealTime() << "s
. \n";
std::ofstream out_file("parameters1
.txt", std::
ios::app);
if (out_file
.is_open()) {
out_file << "\n##################### " << "\n";
out_file << std::scientific << std::set
precision(10);
out_file << "k = 1" << "\n";
out_file << "T = " << t_final << "\n";
out_file << "dt = " << dt << "\n";
out_file << "NN = " << NN << "\n";
out_file << "err_u = " << err_u << "\n";
out_file << "err_p_L2 = " << err_p << "\n";
out_file << "err_theta_L2 = " << err_theta << "\n";
out_file << "H1_error = "<< GradError << "\n";
out_file << "It took " << chrono
.RealTime() << "\n";
out_file
.close();
std::cout << "参数已保存到 parameters1
.txt" << std::endl;
} else {
std::cerr << "错误：无法打开文件！" << std::endl;
}
}
// Free the used memory
.
delete V_space;
delete P_space;
delete T_space;
delete P2;
delete P1;
delete pmesh;
return 0;
}该代码可以编译但不能运行，运行后出现下面的报错Options used:
--mesh
../data/star
.mesh
--refine -1
--order 1
--serial-format
--no-partial-assembly
--device cpu
--visualization
--no-ad
ios2-streams
--dt 0
.05
--NN 4
--t_final 0
.05
Device configuration: cpu
Memory configuration: host-std
Parallel Mesh Stats:
minimum average maximum total
vertices 48 62 77 125
edges 268 302 336 604
faces 411 432 453 864
elements 191 192 193 384
neighbors 1 1 1
minimum maximum
h 0
.280616 0
.280616
kappa 2
.41421 2
.41421
***********************************************************
dim(R) = 2187
dim(W) = 125
dim(T) = 125
dim(R+W+T) = 2437
***********************************************************
开始时间循环
...
时间步 1/1
[LAPTOP-QL5QOURA:00495] *** Process received signal ***
[LAPTOP-QL5QOURA:00495] Signal: Segmentation fault (11)
[LAPTOP-QL5QOURA:00495] Signal code: Address not mapped (1)
[LAPTOP-QL5QOURA:00495] Failing at address: (nil)
[LAPTOP-QL5QOURA:00496] *** Process received signal ***
[LAPTOP-QL5QOURA:00496] Signal: Segmentation fault (11)
[LAPTOP-QL5QOURA:00496] Signal code: Address not mapped (1)
[LAPTOP-QL5QOURA:00496] Failing at address: (nil)
[LAPTOP-QL5QOURA:00495] [ 0] /lib/x86_64-linux-gnu/libc
.so
.6(+0x45330)[0x7e1f49645330]
[LAPTOP-QL5QOURA:00495] [ 1]
./NS_1(+0x44ca78)[0x5675352dba78]
[LAPTOP-QL5QOURA:00495] [ 2]
./NS_1(+0x44179f)[0x5675352d079f]
[LAPTOP-QL5QOURA:00495] [ 3]
./NS_1(+0x4986ed)[0x5675353276ed]
[LAPTOP-QL5QOURA:00495] [ 4]
./NS_1(+0xa06ef)[0x567534f2f6ef]
[LAPTOP-QL5QOURA:00495] [ 5] [LAPTOP-QL5QOURA:00496] [ 0] /lib/x86_64-linux-gnu/libc
.so
.6(+0x45330)[0x7e6862445330]
[LAPTOP-QL5QOURA:00496] [ 1]
./NS_1(+0x44ca78)[0x555f7df7ea78]
[LAPTOP-QL5QOURA:00496] [ 2]
./NS_1(+0x44179f)[0x555f7df7379f]
[LAPTOP-QL5QOURA:00496] [ 3]
./NS_1(+0x4986ed)[0x555f7dfca6ed]
[LAPTOP-QL5QOURA:00496] [ 4]
./NS_1(+0xa06ef)[0x555f7dbd26ef]
[LAPTOP-QL5QOURA:00496] [ 5] /lib/x86_64-linux-gnu/libc
.so
.6(+0x2a1ca)[0x7e686242a1ca]
[LAPTOP-QL5QOURA:00496] [ 6] /lib/x86_64-linux-gnu/libc
.so
.6(__libc_start_main+0x8b)[0x7e686242a28b]
[LAPTOP-QL5QOURA:00496] [ 7]
./NS_1(+0xa7e15)[0x555f7dbd9e15]
[LAPTOP-QL5QOURA:00496] *** End of error message ***
/lib/x86_64-linux-gnu/libc
.so
.6(+0x2a1ca)[0x7e1f4962a1ca]
[LAPTOP-QL5QOURA:00495] [ 6] /lib/x86_64-linux-gnu/libc
.so
.6(__libc_start_main+0x8b)[0x7e1f4962a28b]
[LAPTOP-QL5QOURA:00495] [ 7]
./NS_1(+0xa7e15)[0x567534f36e15]
[LAPTOP-QL5QOURA:00495] *** End of error message ***
--------------------------------------------------------------------------
Primary job terminated normally, but 1 process returned
a non-zero exit code
. Per user-direction, the job has been aborted
.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 0 on node LAPTOP-QL5QOURA exited on signal 11 (Segmentation fault)
.
--------------------------------------------------------------------------给我提供详细的修改方案使其可以运行
由于博客内容为空，暂无法提供包含关键信息的摘要。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
2866
