Compare commits

...
Author SHA1 Message Date
Stowell, Mark L. ef5f729245 Merge remote-tracking branch 'origin/master' into aniso-diffusion-dev
# Conflicts:
#	fem/bilininteg.hpp
#	fem/coefficient.cpp
#	fem/coefficient.hpp
#	linalg/hypre.cpp
#	linalg/hypre.hpp
#	linalg/solvers.cpp
#	makefile
#	miniapps/common/pfem_extras.hpp
#	miniapps/electromagnetics/tesla_solver.hpp
2025-03-13 17:43:24 -07:00
Stowell, Mark L 6fa5a0b096 Merge remote-tracking branch 'origin/master' into aniso-diffusion-dev
# Conflicts:
#	fem/coefficient.cpp
#	fem/coefficient.hpp
2019-04-01 11:27:41 -07:00
Stowell, Mark L 3bd8349909 Merge remote-tracking branch 'origin/master' into aniso-diffusion-dev 2018-11-08 18:00:03 -08:00
Stowell, Mark L af82ee8560 Merge remote-tracking branch 'origin/master' into aniso-diffusion-dev 2018-10-20 12:18:07 -07:00
Stowell, Mark L 2ac542c720 Attempting to support 2D curl cleaning 2018-10-20 12:17:00 -07:00
Mark L. Stowell e6f828a5fe Attempting to add curl free projection... 2018-10-18 13:07:44 -07:00
Mark L. Stowell 4d756edd80 Adding DivergenceFree/Irrotational projectors for RT spaces 2018-10-18 10:29:50 -07:00
Stowell, Mark L cf5bd1f5cc make style 2018-10-18 00:17:32 -07:00
Stowell, Mark L 7173dd2002 Adding H1 diffusion solver 2018-10-18 00:16:56 -07:00
Stowell, Mark L 50182bf440 Adding perturbed elliptic case 2018-10-17 19:51:12 -07:00
Stowell, Mark L 24bfcc5165 Initializing a solution vector before solve 2018-10-17 10:20:38 -07:00
Stowell, Mark L 302f22f297 Switching to analytic evaluation of b vector field 2018-10-16 15:34:29 -07:00
Stowell, Mark L 1637fcd933 Adding argument to control lower bound of mesh size 2018-10-16 13:13:56 -07:00
Stowell, Mark L 2563506174 make style 2018-10-14 16:14:31 -07:00
Stowell, Mark L 60640c3f7e Adding computation of full thermal flux 2018-10-14 10:39:39 -07:00
Stowell, Mark L f221521203 make style 2018-10-14 10:12:37 -07:00
Stowell, Mark L 1d9e736af6 Adding a miniapp which solve for thermal flux in HDiv 2018-10-14 10:10:55 -07:00
Stowell, Mark L d80dbfd99a Adding flux computation 2018-10-10 16:46:17 -07:00
Stowell, Mark L 3f44043e60 Adding steady state anisotropic diffusion solver 2018-10-10 12:48:07 -07:00
Stowell, Mark L b218959bca Inserting the thermal flux solver 2018-10-03 10:33:45 -07:00
Stowell, Mark L aee7bc9d43 Adding first draft of hybrid diffusion solver 2018-10-01 16:14:53 -07:00
Stowell, Mark L 31cac320d4 Bugfix in activation of nonlinear solver 2018-10-01 15:08:15 -07:00
Stowell, Mark L e7e0fb0a88 Adding a specialized miniapp to duplicate results from the van Es papper 2018-09-30 21:47:54 -07:00
Stowell, Mark L ba71d13980 Merge remote-tracking branch 'origin/master' into aniso-diffusion-dev 2018-09-25 13:04:15 -07:00
Stowell, Mark L 5c326a5535 Fixing a typo in a comment 2018-09-25 13:02:39 -07:00
Stowell, Mark L 9457f7e5b6 Switching to nonlinear solver 2018-09-24 15:45:49 -07:00
Stowell, Mark L ff030ee970 Adding another time dependent test case 2018-09-24 12:45:36 -07:00
Stowell, Mark L 2be9e1f36c Adding a steady state solver to the thermal miniapps 2018-09-24 12:45:03 -07:00
Stowell, Mark L 7671cd9f36 Merge remote-tracking branch 'origin/master' into aniso-diffusion-dev 2018-09-16 13:17:34 -07:00
Stowell, Mark L f89a633fda make style 2018-09-14 14:41:13 -07:00
Stowell, Mark L 4ce1cef6b8 Adding SetOperator methods to HyprePCG, HypreGMRES, HypreDiagScale, and HypreParaSails 2018-09-14 14:34:13 -07:00
Stowell, Mark L c667bf3025 Fixing HypreGMRES::SetOperator method in the presence of a preconditioner 2018-09-14 13:47:29 -07:00
Stowell, Mark L e1678afe40 Using new HypreGMRES with SetOperator method 2018-09-10 16:46:25 -07:00
Stowell, Mark L c0291398ed Implementing HypreGMRES::SetOperator method 2018-09-10 16:45:59 -07:00
Stowell, Mark L 0b4f10d79d Debugging gradient check 2018-09-09 16:29:15 -07:00
Stowell, Mark L 8dfd0e1547 Adding NewtonSolver method to validate gradient 2018-09-09 16:28:18 -07:00
Stowell, Mark L caf239c99a Updating with time dependent source and exact solution 2018-09-09 00:39:58 -07:00
Stowell, Mark L 44c33aece0 Merge remote-tracking branch 'origin/master' into aniso-diffusion-dev 2018-09-08 23:59:16 -07:00
Stowell, Mark L ed49856390 Merge branch 'aniso-diffusion-dev' of github.com:mfem/mfem into aniso-diffusion-dev 2018-09-06 14:13:18 -07:00
Stowell, Mark L 4fef6ca298 make style 2018-09-06 14:12:28 -07:00
Stowell, Mark L d9d809e81c Adding "thermal" to miniapps subdirectories 2018-09-06 14:12:16 -07:00
Stowell, Mark L 51e85ccd84 Fixing nonlinear solve and applying 'make style' 2018-09-06 14:11:54 -07:00
Mark L. Stowell 4304159303 Merge branch 'aniso-diffusion-dev' of github.com:mfem/mfem into aniso-diffusion-dev 2018-09-05 16:56:05 -07:00
Stowell, Mark L 6276268e52 Retain zeros to maintain sparsity pattern 2018-09-05 16:54:59 -07:00
Mark L. Stowell 580ae34842 Retaining zeros to maintain sparsity pattern 2018-09-05 16:51:10 -07:00
Stowell, Mark L 97eaf8efbc Parallelizing the linear solves 2018-09-05 15:48:28 -07:00
Stowell, Mark L e6621c9b0c Parallel bug 2018-09-05 15:25:18 -07:00
Stowell, Mark L 461246f80e Adding a missing overload 2018-09-05 13:29:10 -07:00
Stowell, Mark L a5941ee72f Bugfix: reinitializing matrices before reassembling 2018-09-05 11:10:09 -07:00
Stowell, Mark L c4c2ceab59 Linear case now working 2018-09-05 10:15:17 -07:00
Stowell, Mark L 88b99a1719 Fixing vector dimension in vector grid functions 2018-09-05 10:14:52 -07:00
Stowell, Mark L 7aa7b4ee53 Adjusting initialization order so that vector size is known earlier 2018-09-03 11:19:27 -07:00
Stowell, Mark L bcf87fee29 Modifying VectorGridFunctionCoefs to accept NULL pointers 2018-09-03 11:06:32 -07:00
Stowell, Mark L 2949dc5a46 Adding makefile for miniapps/thermal 2018-09-03 10:48:11 -07:00
Stowell, Mark L b0dbadd007 Adding first draft of non-linear thermal diffusion solver 2018-09-03 10:19:34 -07:00
Stowell, Mark L 6ca1f95979 Adding scalar multiplication by a constant 2018-08-31 22:54:36 -07:00
Stowell, Mark L 39794585c4 Adding an Identity Matrix Coefficient 2018-08-31 16:50:02 -07:00
Stowell, Mark L 65a71259f1 Merge remote-tracking branch 'origin/elementwise-error-dev' into aniso-diffusion-dev 2018-08-31 16:49:37 -07:00
Stowell, Mark L 3f4e8324d4 Adding a coefficient which computes a unit vector field from a vector field 2018-08-29 14:26:29 -07:00
Stowell, Mark L 9fca398741 Adding ability to alter derived coefficients 2018-08-29 13:59:38 -07:00
Stowell, Mark L a2b8f7a129 Merge remote-tracking branch 'origin/master' into aniso-diffusion-dev 2018-08-29 09:28:15 -07:00
Stowell, Mark L a367631ce5 Adding various coefficients which are sums or products of other coefficients 2018-08-28 16:34:13 -07:00
Stowell, Mark L d7718f5c57 make style 2018-08-28 14:51:48 -07:00
Stowell, Mark L fe88c4685d Adding coefficients to compute div, grad, or curl of grid functions. 2018-08-28 14:24:15 -07:00
24 changed files with 11066 additions and 43 deletions
File diff suppressed because it is too large Load Diff
+33
View File
@@ -2003,6 +2003,39 @@ void NewtonSolver::Mult(const Vector &b, Vector &x) const
Monitor(final_iter, final_norm, r, x, true);
}
double NewtonSolver::CheckGradient(const Vector &x, const Vector &h) const
{
Vector x1(x.Size());
Vector b0(x.Size());
// Evaluate operator and its gradient at x
oper->Mult(x, b0);
oper->GetGradient(x).Mult(h, c);
// Evaluate operator at x+h
add(x, 1.0, h, x1);
oper->Mult(x1, r);
// Compute error in F(x) + G * h
r.Add(-1.0, b0);
r.Add(-1.0, c);
double norm1 = Norm(r);
// Evaluate operator at x+h/2
add(x, 0.5, h, x1);
oper->Mult(x1, r);
// Compute error in F(x) + G * h / 2
r.Add(-1.0, b0);
r.Add(-0.5, c);
double norm2 = Norm(r);
if (norm1 == 0.0 ) { return -1.0; }
return 2.0 * norm2 / norm1;
}
void NewtonSolver::SetAdaptiveLinRtol(const int type,
const real_t rtol0,
const real_t rtol_max,
+10
View File
@@ -737,6 +737,16 @@ public:
/** If `b.Size() != Height()`, then @a b is assumed to be zero. */
void Mult(const Vector &b, Vector &x) const override;
/// Verify that the operator returns a valid gradient
/** The gradient should satisfy the definition of a Frechet Derivative
i.e. lim_{h->0} ||F(x+H)-F(x)-G(x)*h||/||h|| = 0. This method
returns 2 * ||F(x+h/2)-F(x)-G(x)*h/2|| / ||F(x+h)-F(x)-G(x)*h||
which should be less than or equal to 1 for any valid gradient
provided h is sufficiently small. This method returns -1 if the
operator appears to be linear in which case the ratio would be 0/0.
*/
virtual double CheckGradient(const Vector &x, const Vector &h) const;
/** @brief This method can be overloaded in derived classes to implement line
search algorithms. */
/** The base class implementation (NewtonSolver) simply returns 1. A return
+1 -1
View File
@@ -125,7 +125,7 @@ EXAMPLE_TEST_DIRS := examples
MINIAPP_SUBDIRS = common electromagnetics meshing navier performance tools \
toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff hooke \
multidomain dpg hdiv-linear-solver spde
multidomain dpg hdiv-linear-solver spde thermal
MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
+205 -24
View File
@@ -94,13 +94,13 @@ ParDiscreteDivOperator::ParDiscreteDivOperator(ParFiniteElementSpace *dfes,
this->AddDomainInterpolator(new DivergenceInterpolator);
}
IrrotationalProjector
::IrrotationalProjector(ParFiniteElementSpace & H1FESpace,
ParFiniteElementSpace & HCurlFESpace,
const int & irOrder,
ParBilinearForm * s0,
ParMixedBilinearForm * weakDiv,
ParDiscreteGradOperator * grad)
IrrotationalNDProjector
::IrrotationalNDProjector(ParFiniteElementSpace & H1FESpace,
ParFiniteElementSpace & HCurlFESpace,
const int & irOrder,
ParBilinearForm * s0,
ParMixedBilinearForm * weakDiv,
ParDiscreteGradOperator * grad)
: H1FESpace_(&H1FESpace),
HCurlFESpace_(&HCurlFESpace),
s0_(s0),
@@ -152,7 +152,7 @@ IrrotationalProjector
xDiv_ = new ParGridFunction(H1FESpace_);
}
IrrotationalProjector::~IrrotationalProjector()
IrrotationalNDProjector::~IrrotationalNDProjector()
{
delete psi_;
delete xDiv_;
@@ -167,7 +167,7 @@ IrrotationalProjector::~IrrotationalProjector()
}
void
IrrotationalProjector::InitSolver() const
IrrotationalNDProjector::InitSolver() const
{
delete pcg_;
delete amg_;
@@ -182,7 +182,7 @@ IrrotationalProjector::InitSolver() const
}
void
IrrotationalProjector::Mult(const Vector &x, Vector &y) const
IrrotationalNDProjector::Mult(const Vector &x, Vector &y) const
{
// Compute the divergence of x
weakDiv_->Mult(x,*xDiv_); *xDiv_ *= -1.0;
@@ -203,7 +203,7 @@ IrrotationalProjector::Mult(const Vector &x, Vector &y) const
}
void
IrrotationalProjector::Update()
IrrotationalNDProjector::Update()
{
delete pcg_; pcg_ = NULL;
delete amg_; amg_ = NULL;
@@ -234,31 +234,212 @@ IrrotationalProjector::Update()
H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
}
DivergenceFreeProjector
::DivergenceFreeProjector(ParFiniteElementSpace & H1FESpace,
ParFiniteElementSpace & HCurlFESpace,
const int & irOrder,
ParBilinearForm * s0,
ParMixedBilinearForm * weakDiv,
ParDiscreteGradOperator * grad)
: IrrotationalProjector(H1FESpace,HCurlFESpace, irOrder, s0, weakDiv, grad)
DivergenceFreeNDProjector
::DivergenceFreeNDProjector(ParFiniteElementSpace & H1FESpace,
ParFiniteElementSpace & HCurlFESpace,
const int & irOrder,
ParBilinearForm * s0,
ParMixedBilinearForm * weakDiv,
ParDiscreteGradOperator * grad)
: IrrotationalNDProjector(H1FESpace,HCurlFESpace, irOrder, s0, weakDiv, grad)
{}
DivergenceFreeProjector::~DivergenceFreeProjector()
DivergenceFreeNDProjector::~DivergenceFreeNDProjector()
{}
void
DivergenceFreeProjector::Mult(const Vector &x, Vector &y) const
DivergenceFreeNDProjector::Mult(const Vector &x, Vector &y) const
{
this->IrrotationalProjector::Mult(x, y);
this->IrrotationalNDProjector::Mult(x, y);
y -= x;
y *= -1.0;
}
void
DivergenceFreeProjector::Update()
DivergenceFreeNDProjector::Update()
{
this->IrrotationalProjector::Update();
this->IrrotationalNDProjector::Update();
}
DivergenceFreeRTProjector
::DivergenceFreeRTProjector(ParFiniteElementSpace & HCurlFESpace,
ParFiniteElementSpace & HDivFESpace,
const int & irOrder,
ParBilinearForm * s1,
ParMixedBilinearForm * weakCurl,
ParDiscreteCurlOperator * curl)
: HCurlFESpace_(&HCurlFESpace),
HDivFESpace_(&HDivFESpace),
s1_(s1),
weakCurl_(weakCurl),
curl_(curl),
psi_(NULL),
xCurl_(NULL),
S1_(NULL),
pc_(NULL),
pcg_(NULL),
dim_(HCurlFESpace_->GetFE(0)->GetDim()),
ownsS1_(s1 == NULL),
ownsWeakCurl_(weakCurl == NULL),
ownsCurl_(curl == NULL)
{
ess_bdr_.SetSize(HCurlFESpace_->GetParMesh()->bdr_attributes.Max());
ess_bdr_ = 1;
HCurlFESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
int geom = HCurlFESpace_->GetFE(0)->GetGeomType();
const IntegrationRule * ir = &IntRules.Get(geom, irOrder);
if ( s1 == NULL )
{
s1_ = new ParBilinearForm(HCurlFESpace_);
BilinearFormIntegrator * ccInteg = (dim_==2) ?
dynamic_cast<BilinearFormIntegrator*>(new DiffusionIntegrator) :
dynamic_cast<BilinearFormIntegrator*>(new CurlCurlIntegrator);
ccInteg->SetIntRule(ir);
s1_->AddDomainIntegrator(ccInteg);
s1_->Assemble();
s1_->Finalize();
S1_ = new HypreParMatrix;
}
if ( weakCurl_ == NULL )
{
weakCurl_ = new ParMixedBilinearForm(HDivFESpace_, HCurlFESpace_);
BilinearFormIntegrator * wcurlInteg = new MixedVectorWeakCurlIntegrator;
wcurlInteg->SetIntRule(ir);
weakCurl_->AddDomainIntegrator(wcurlInteg);
weakCurl_->Assemble();
weakCurl_->Finalize();
}
if ( curl_ == NULL )
{
curl_ = new ParDiscreteCurlOperator(HCurlFESpace_, HDivFESpace_);
curl_->Assemble();
curl_->Finalize();
}
psi_ = new ParGridFunction(HCurlFESpace_);
xCurl_ = new ParGridFunction(HCurlFESpace_);
}
DivergenceFreeRTProjector::~DivergenceFreeRTProjector()
{
delete psi_;
delete xCurl_;
delete pc_;
delete pcg_;
delete S1_;
delete s1_;
delete weakCurl_;
}
void
DivergenceFreeRTProjector::InitSolver() const
{
delete pcg_;
delete pc_;
if (dim_ == 2)
{
HypreBoomerAMG * amg = new HypreBoomerAMG(*S1_);
amg->SetPrintLevel(0);
pc_ = amg;
}
else
{
HypreAMS * ams = new HypreAMS(*S1_, HCurlFESpace_);
ams->SetPrintLevel(0);
pc_ = ams;
}
pcg_ = new HyprePCG(*S1_);
pcg_->SetTol(1e-14);
pcg_->SetMaxIter(200);
pcg_->SetPrintLevel(0);
pcg_->SetPreconditioner(*pc_);
}
void
DivergenceFreeRTProjector::Mult(const Vector &x, Vector &y) const
{
// Compute the curl of x
weakCurl_->Mult(x,*xCurl_);
// Apply essential BC and form linear system
*psi_ = 0.0;
s1_->FormLinearSystem(ess_bdr_tdofs_, *psi_, *xCurl_, *S1_, Psi_, RHS_);
// Solve the linear system for Psi
if ( pcg_ == NULL ) { this->InitSolver(); }
pcg_->Mult(RHS_, Psi_);
// Compute the parallel grid function correspoinding to Psi
s1_->RecoverFEMSolution(Psi_, *xCurl_, *psi_);
// Compute the divergence free portion of x
curl_->Mult(*psi_, y);
}
void
DivergenceFreeRTProjector::Update()
{
delete pcg_; pcg_ = NULL;
delete pc_; pc_ = NULL;
delete S1_; S1_ = new HypreParMatrix;
psi_->Update();
xCurl_->Update();
if ( ownsS1_ )
{
s1_->Update();
s1_->Assemble();
s1_->Finalize();
}
if ( ownsWeakCurl_ )
{
weakCurl_->Update();
weakCurl_->Assemble();
weakCurl_->Finalize();
}
if ( ownsCurl_ )
{
curl_->Update();
curl_->Assemble();
curl_->Finalize();
}
HCurlFESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
}
IrrotationalRTProjector
::IrrotationalRTProjector(ParFiniteElementSpace & HCurlFESpace,
ParFiniteElementSpace & HDivFESpace,
const int & irOrder,
ParBilinearForm * s1,
ParMixedBilinearForm * weakCurl,
ParDiscreteCurlOperator * curl)
: DivergenceFreeRTProjector(HCurlFESpace, HDivFESpace, irOrder,
s1, weakCurl, curl)
{}
IrrotationalRTProjector::~IrrotationalRTProjector()
{}
void
IrrotationalRTProjector::Mult(const Vector &x, Vector &y) const
{
this->DivergenceFreeRTProjector::Mult(x, y);
y -= x;
y *= -1.0;
}
void
IrrotationalRTProjector::Update()
{
this->DivergenceFreeRTProjector::Update();
}
void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
+89 -16
View File
@@ -115,16 +115,16 @@ public:
/// This class computes the irrotational portion of a vector field.
/// This vector field must be discretized using Nedelec basis
/// functions.
class IrrotationalProjector : public Operator
class IrrotationalNDProjector : public Operator
{
public:
IrrotationalProjector(ParFiniteElementSpace & H1FESpace,
ParFiniteElementSpace & HCurlFESpace,
const int & irOrder,
ParBilinearForm * s0 = NULL,
ParMixedBilinearForm * weakDiv = NULL,
ParDiscreteGradOperator * grad = NULL);
virtual ~IrrotationalProjector();
IrrotationalNDProjector(ParFiniteElementSpace & H1FESpace,
ParFiniteElementSpace & HCurlFESpace,
const int & irOrder,
ParBilinearForm * s0 = NULL,
ParMixedBilinearForm * weakDiv = NULL,
ParDiscreteGradOperator * grad = NULL);
virtual ~IrrotationalNDProjector();
// Given a GridFunction 'x' of Nedelec DoFs for an arbitrary vector field,
// compute the Nedelec DoFs of the irrotational portion, 'y', of
@@ -164,16 +164,16 @@ private:
/// This class computes the divergence free portion of a vector field.
/// This vector field must be discretized using Nedelec basis
/// functions.
class DivergenceFreeProjector : public IrrotationalProjector
class DivergenceFreeNDProjector : public IrrotationalNDProjector
{
public:
DivergenceFreeProjector(ParFiniteElementSpace & H1FESpace,
ParFiniteElementSpace & HCurlFESpace,
const int & irOrder,
ParBilinearForm * s0 = NULL,
ParMixedBilinearForm * weakDiv = NULL,
ParDiscreteGradOperator * grad = NULL);
virtual ~DivergenceFreeProjector();
DivergenceFreeNDProjector(ParFiniteElementSpace & H1FESpace,
ParFiniteElementSpace & HCurlFESpace,
const int & irOrder,
ParBilinearForm * s0 = NULL,
ParMixedBilinearForm * weakDiv = NULL,
ParDiscreteGradOperator * grad = NULL);
virtual ~DivergenceFreeNDProjector();
// Given a vector 'x' of Nedelec DoFs for an arbitrary vector field,
// compute the Nedelec DoFs of the divergence free portion, 'y', of
@@ -185,6 +185,79 @@ public:
};
/// This class computes the divergence free portion of a vector field.
/// This vector field must be discretized using Raviart-Thomas basis
/// functions.
class DivergenceFreeRTProjector : public Operator
{
public:
DivergenceFreeRTProjector(ParFiniteElementSpace & HCurlFESpace,
ParFiniteElementSpace & HDivFESpace,
const int & irOrder,
ParBilinearForm * s1 = NULL,
ParMixedBilinearForm * weakCurl = NULL,
ParDiscreteCurlOperator * curl = NULL);
virtual ~DivergenceFreeRTProjector();
// Given a GridFunction 'x' of Raviart-Thomas DoFs for an arbitrary vector
// field, compute the Raviart-Thomas DoFs of the divergence free portion,
// 'y', of this vector field. The resulting GridFunction will satisfy
// Div y = 0 to machine precision.
virtual void Mult(const Vector &x, Vector &y) const;
void Update();
private:
void InitSolver() const;
ParFiniteElementSpace * HCurlFESpace_;
ParFiniteElementSpace * HDivFESpace_;
ParBilinearForm * s1_;
ParMixedBilinearForm * weakCurl_;
ParDiscreteCurlOperator * curl_;
ParGridFunction * psi_;
ParGridFunction * xCurl_;
HypreParMatrix * S1_;
mutable Vector Psi_;
mutable Vector RHS_;
mutable HypreSolver * pc_;
mutable HyprePCG * pcg_;
Array<int> ess_bdr_, ess_bdr_tdofs_;
int dim_;
bool ownsS1_;
bool ownsWeakCurl_;
bool ownsCurl_;
};
/// This class computes the irrotational portion of a vector field.
/// This vector field must be discretized using Nedelec basis
/// functions.
class IrrotationalRTProjector : public DivergenceFreeRTProjector
{
public:
IrrotationalRTProjector(ParFiniteElementSpace & HCurlFESpace,
ParFiniteElementSpace & HDivFESpace,
const int & irOrder,
ParBilinearForm * s1 = NULL,
ParMixedBilinearForm * weakCurl = NULL,
ParDiscreteCurlOperator * curl = NULL);
virtual ~IrrotationalRTProjector();
// Given a GridFunction 'x' of Raviart-Thomas DoFs for an arbitrary vector
// field, compute the Raviart-Thomas DoFs of the irrotational portion,
// 'y', of this vector field. The resulting GridFunction will satisfy
// Curl y = 0 to machine precision.
virtual void Mult(const Vector &x, Vector &y) const;
void Update();
};
/// Visualize the given parallel mesh object, using a GLVis server on the
/// specified host and port. Set the visualization window title, and optionally,
/// its geometry.
+2 -2
View File
@@ -152,8 +152,8 @@ TeslaSolver::TeslaSolver(ParMesh & pmesh, int order,
{
jr_ = new ParGridFunction(HCurlFESpace_);
j_ = new ParGridFunction(HCurlFESpace_);
DivFreeProj_ = new DivergenceFreeProjector(*H1FESpace_, *HCurlFESpace_,
irOrder, NULL, NULL, grad_);
DivFreeProj_ = new DivergenceFreeNDProjector(*H1FESpace_, *HCurlFESpace_,
irOrder, NULL, NULL, grad_);
}
if ( kbcs.Size() > 0 )
+691
View File
@@ -0,0 +1,691 @@
// MFEM Example 1 - Parallel Version
//
// Compile with: make ex1p
//
// Sample runs: mpirun -np 4 ex1p -m ../data/square-disc.mesh
// mpirun -np 4 ex1p -m ../data/star.mesh
// mpirun -np 4 ex1p -m ../data/escher.mesh
// mpirun -np 4 ex1p -m ../data/fichera.mesh
// mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
// mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
// mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
// mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
// mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
// mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
// mpirun -np 4 ex1p -m ../data/star-surf.mesh
// mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
// mpirun -np 4 ex1p -m ../data/inline-segment.mesh
// mpirun -np 4 ex1p -m ../data/amr-quad.mesh
// mpirun -np 4 ex1p -m ../data/amr-hex.mesh
// mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
// mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
//
// Description: This example code demonstrates the use of MFEM to define a
// simple finite element discretization of the Laplace problem
// -Delta u = 1 with homogeneous Dirichlet boundary conditions.
// Specifically, we discretize using a FE space of the specified
// order, or if order < 1 using an isoparametric/isogeometric
// space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
// NURBS mesh, etc.)
//
// The example highlights the use of mesh refinement, finite
// element grid functions, as well as linear and bilinear forms
// corresponding to the left-hand side and right-hand side of the
// discrete linear system. We also cover the explicit elimination
// of essential boundary conditions, static condensation, and the
// optional connection to the GLVis tool for visualization.
#include "mfem.hpp"
#include "../common/pfem_extras.hpp"
#include <fstream>
#include <iostream>
using namespace std;
using namespace mfem;
static double nl_exp_ = 2.5;
static double theta_ = 0.0;
static double chi_perp_ = 1.0;
static double chi_para_min_ = 100.0;
static double chi_para_max_ = 1000.0;
double uFunc(const Vector &x)
{
return sin(M_PI * x[0]) * sin(M_PI * x[1]);
}
double QFunc(const Vector &x)
{
double chi_ratio = (nl_exp_ > 0.0) ?
pow(chi_para_min_ / chi_para_max_, 1.0 / nl_exp_) : 1.0;
double u = uFunc(x);
double T = chi_ratio + (1.0 - chi_ratio) * u;
double cx = cos(M_PI * x[0]);
double sx = sin(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sy = sin(M_PI * x[1]);
double ct = cos(theta_);
double st = sin(theta_);
double s2t = sin(2.0 * theta_);
return M_PI * M_PI * (chi_perp_ * (u + cx * cy * s2t) +
chi_para_max_ * (u - cx * cy * s2t) * pow(T, nl_exp_) +
chi_para_max_ * nl_exp_ * (1.0 - chi_ratio) *
(u * u - sx * sx * st * st - sy * sy * ct * ct -
u * cx * cy * s2t) * pow(T, nl_exp_ - 1.0) );
}
void unitVectorField(const Vector &, Vector &u)
{
u.SetSize(2);
u[0] = cos(theta_);
u[1] = sin(theta_);
}
class ChiParaCoef : public MatrixCoefficient
{
private:
MatrixCoefficient * bbT_;
GridFunctionCoefficient * T_;
double nl_exp_;
double chi_min_;
double chi_max_;
double gamma_;
public:
ChiParaCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double nl_exp, double chi_min, double chi_max)
: MatrixCoefficient(2), bbT_(&bbT), T_(&T), nl_exp_(nl_exp),
chi_min_(chi_min), chi_max_(chi_max),
gamma_(pow(chi_min/chi_max, 1.0 / nl_exp_))
{
// cout << "(chi_min/chi_max)^nl_exp = " << gamma_ << endl;
}
void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
if ( nl_exp_ == 0.0)
{
K *= chi_max_;
}
else
{
double Tval = T_->Eval(T, ip);
// cout << "Tval = " << Tval << endl;
// cout << "Multiplier: " << pow(gamma_ + (1.0 - gamma_) * Tval, nl_exp_) << endl;
double u = gamma_ + (1.0 - gamma_) * Tval;
u = max(gamma_, min(u, 1.0));
K *= chi_max_ * pow(u, nl_exp_);
}
}
};
class ChiCoef : public MatrixSumCoefficient
{
private:
ChiParaCoef * chiParaCoef_;
public:
ChiCoef(MatrixCoefficient & chiPerp, ChiParaCoef & chiPara)
: MatrixSumCoefficient(chiPerp, chiPara), chiParaCoef_(&chiPara) {}
void SetTemp(GridFunction & T) { chiParaCoef_->SetTemp(T); }
};
class dChiCoef : public MatrixCoefficient
{
private:
MatrixCoefficient * bbT_;
GridFunctionCoefficient * T_;
double nl_exp_;
double chi_min_;
double chi_max_;
double gamma_;
public:
dChiCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double nl_exp, double chi_min, double chi_max)
: MatrixCoefficient(2), bbT_(&bbT), T_(&T), nl_exp_(nl_exp),
chi_min_(chi_min), chi_max_(chi_max),
gamma_(pow(chi_min/chi_max, 1.0 / nl_exp_))
{}
void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
double Tval = T_->Eval(T, ip);
double u = gamma_ + (1.0 - gamma_) * Tval;
u = max(gamma_, min(u, 1.0));
K *= nl_exp_ * chi_max_ * (1.0 - gamma_) * pow(u, nl_exp_ - 1.0);
}
};
class ImplicitDiffOp : public Operator
{
public:
ImplicitDiffOp(ParFiniteElementSpace & H1_FESpace,
Coefficient & TBdr,
Array<int> & bdr_attr,
ChiCoef & chi,
dChiCoef & dchi,
Coefficient & heatSource);
~ImplicitDiffOp();
// void SetState(ParGridFunction & T);
void Mult(const Vector &x, Vector &y) const;
Operator & GetGradient(const Vector &x) const;
Solver & GetGradientSolver() const;
const Vector & GetRHS() const { return RHS_; }
private:
bool first_;
// bool nonLinear_;
Array<int> & ess_bdr_attr_;
Array<int> ess_bdr_tdofs_;
Coefficient * bdrCoef_;
ChiCoef * chiCoef_;
dChiCoef * dChiCoef_;
Coefficient * QCoef_;
// ScalarMatrixProductCoefficient dtChiCoef_;
mutable ParGridFunction T_;
// mutable ParGridFunction T1_;
// mutable ParGridFunction dT_;
mutable GradientGridFunctionCoefficient gradTCoef_;
// ScalarVectorProductCoefficient dtGradTCoef_;
// MatVecCoefficient dtdChiGradTCoef_;
MatVecCoefficient dChiGradTCoef_;
mutable ParBilinearForm s0chi_;
mutable ParBilinearForm a0_;
mutable HypreParMatrix A_;
// mutable ParGridFunction dTdt_;
mutable ParLinearForm Q_;
mutable ParLinearForm Qs_;
mutable ParLinearForm rhs_;
mutable Vector SOL_;
mutable Vector RHS_;
// Vector RHS0_; // Dummy RHS vector which hase length zero
mutable Solver * AInv_;
mutable HypreBoomerAMG * APrecond_;
};
int main(int argc, char *argv[])
{
// 1. Initialize MPI.
int num_procs, myid;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
// 2. Parse command-line options.
int n = 1;
int el_type = Element::QUADRILATERAL;
int order = 1;
int max_iter = 100;
int ser_ref_levels = 0;
int par_ref_levels = 0;
bool static_cond = false;
bool visualization = 1;
OptionsParser args(argc, argv);
args.AddOption(&n, "-n", "--num-elems-1d",
"Number of elements in x and y directions. "
"Total number of elements is n^2.");
args.AddOption(&el_type, "-e", "--element-type",
"Element type: 2-Triangle, 3-Quadrilateral.");
args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
"Number of times to refine the mesh uniformly in serial.");
args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
"Number of times to refine the mesh uniformly in parallel.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree) or -1 for"
" isoparametric space.");
args.AddOption(&max_iter, "-mit", "--max-iter",
"Maximum number of Newton iterations.");
args.AddOption(&chi_perp_, "-chi-perp", "--chi-perpendicular",
"Chi_perp.");
args.AddOption(&chi_para_max_, "-chi-max", "--chi-para-max",
"Maximum value of chi along field lines.");
args.AddOption(&chi_para_min_, "-chi-min", "--chi-para-min",
"Minimum value of chi along field lines.");
args.AddOption(&theta_, "-t", "--theta",
"Angle of strong diffusion in degrees.");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
theta_ *= M_PI / 180.0;
// 3. Read the (serial) mesh from the given mesh file on all processors. We
// can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
// and volume meshes with the same code.
Mesh *mesh = new Mesh(n, n, (Element::Type)el_type, 1);
int dim = mesh->Dimension();
// 4. Refine the serial mesh on all processors to increase the resolution. In
// this example we do 'ref_levels' of uniform refinement. We choose
// 'ref_levels' to be the largest number that gives a final mesh with no
// more than 10,000 elements.
for (int lev = 0; lev < ser_ref_levels; lev++)
{
mesh->UniformRefinement();
}
// 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
// this mesh further in parallel to increase the resolution. Once the
// parallel mesh is defined, the serial mesh can be deleted.
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
delete mesh;
for (int lev = 0; lev < par_ref_levels; lev++)
{
pmesh->UniformRefinement();
}
// 6. Define a parallel finite element space on the parallel mesh. Here we
// use continuous Lagrange finite elements of the specified order. If
// order < 1, we instead use an isoparametric/isogeometric space.
FiniteElementCollection *fec;
if (order > 0)
{
fec = new H1_FECollection(order, dim);
}
else if (pmesh->GetNodes())
{
fec = pmesh->GetNodes()->OwnFEC();
if (myid == 0)
{
cout << "Using isoparametric FEs: " << fec->Name() << endl;
}
}
else
{
fec = new H1_FECollection(order = 1, dim);
}
ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
HYPRE_Int size = fespace->GlobalTrueVSize();
if (myid == 0)
{
cout << "Number of finite element unknowns: " << size << endl;
}
L2_FECollection L2FEC0(0, dim);
ParFiniteElementSpace L2FESpace0(pmesh, &L2FEC0);
// 7. Determine the list of true (i.e. parallel conforming) essential
// boundary dofs. In this example, the boundary conditions are defined
// by marking all the boundary attributes from the mesh as essential
// (Dirichlet) and converting them to a list of true dofs.
Array<int> ess_tdof_list;
Array<int> ess_bdr;
if (pmesh->bdr_attributes.Size())
{
ess_bdr.SetSize(pmesh->bdr_attributes.Max());
ess_bdr = 1;
fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}
// 8. Set up the parallel linear form b(.) which corresponds to the
// right-hand side of the FEM linear system, which in this case is
// (1,phi_i) where phi_i are the basis functions in fespace.
ConstantCoefficient zeroCoef(0.0);
ConstantCoefficient oneCoef(1.0);
FunctionCoefficient uCoef(uFunc);
FunctionCoefficient QCoef(QFunc);
ParLinearForm *Q = new ParLinearForm(fespace);
Q->AddDomainIntegrator(new DomainLFIntegrator(QCoef));
Q->Assemble();
// 9. Define the solution vector x as a parallel finite element grid function
// corresponding to fespace. Initialize x with initial guess of zero,
// which satisfies the boundary conditions.
ParGridFunction u(fespace);
ParGridFunction u_error(&L2FESpace0);
ParGridFunction Q_gf(fespace);
//u = 0.0;
u.ProjectCoefficient(uCoef);
Q_gf.ProjectCoefficient(QCoef);
// 10. Set up the parallel bilinear form a(.,.) on the finite element space
// corresponding to the Laplacian operator -Delta, by adding the Diffusion
// domain integrator.
VectorFunctionCoefficient vCoef(2, unitVectorField);
OuterProductCoefficient vvTCoef(vCoef, vCoef);
IdentityMatrixCoefficient ICoef(2);
GridFunctionCoefficient uGFCoef(&u);
ChiParaCoef chiPara(vvTCoef, uGFCoef, nl_exp_, chi_para_min_, chi_para_max_);
MatrixSumCoefficient chiPerp(ICoef, vvTCoef, chi_perp_, -chi_perp_);
ChiCoef chiCoef(chiPerp, chiPara);
dChiCoef dchiCoef(vvTCoef, uGFCoef, nl_exp_, chi_para_min_, chi_para_max_);
ImplicitDiffOp ido(*fespace, zeroCoef, ess_bdr,
chiCoef, dchiCoef, QCoef);
/*
ParBilinearForm *a = new ParBilinearForm(fespace);
a->AddDomainIntegrator(new DiffusionIntegrator(chiCoef));
// 11. Assemble the parallel bilinear form and the corresponding linear
// system, applying any necessary transformations such as: parallel
// assembly, eliminating boundary conditions, applying conforming
// constraints for non-conforming AMR, static condensation, etc.
if (static_cond) { a->EnableStaticCondensation(); }
a->Assemble();
HypreParMatrix A;
Vector Q_dof, u_dof;
a->FormLinearSystem(ess_tdof_list, u, *Q, A, u_dof, Q_dof);
if (myid == 0)
{
cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
}
// 12. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
// preconditioner from hypre.
HypreSolver *amg = new HypreBoomerAMG(A);
HyprePCG *pcg = new HyprePCG(A);
pcg->SetTol(1e-12);
pcg->SetMaxIter(200);
pcg->SetPrintLevel(2);
pcg->SetPreconditioner(*amg);
pcg->Mult(Q_dof, u_dof);
// 13. Recover the parallel grid function corresponding to X. This is the
// local finite element solution on each processor.
a->RecoverFEMSolution(u_dof, *Q, u);
*/
// ido.SetState(u);
Solver & solver = ido.GetGradientSolver();
NewtonSolver newton(MPI_COMM_WORLD);
newton.SetPrintLevel(2);
// newton.SetRelTol(1e-10);
newton.SetAbsTol(1e-10);
// newton.SetMaxIter(max_iter);
newton.SetOperator(ido);
newton.SetSolver(solver);
Vector uVec(fespace->GetTrueVSize());
Vector duVec(fespace->GetTrueVSize());
uVec = 1.0;
duVec = 0.001;
cout << "Gradient verification: " << newton.CheckGradient(uVec, duVec)
<< endl;
uVec = 0.0;
socketstream vis_T, vis_Q, vis_errT;
for (int it = 0; it<max_iter; it++)
{
newton.SetMaxIter(1);
newton.Mult(ido.GetRHS(), uVec);
bool conv = newton.GetConverged();
u.Distribute(uVec);
u.GridFunction::ComputeElementL2Errors(uCoef, u_error);
double err = u.ComputeL2Error(uCoef);
if (myid == 0)
{
cout << "Range of solution vector: "
<< uVec.Min() << " -> " << uVec.Max() << endl;
cout << "L2 Error of Solution: " << err << endl;
}
// 14. Save the refined mesh and the solution in parallel. This output can
// be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
{
ostringstream mesh_name, sol_name;
mesh_name << "mesh." << setfill('0') << setw(6) << myid;
sol_name << "sol." << setfill('0') << setw(6) << myid;
ofstream mesh_ofs(mesh_name.str().c_str());
mesh_ofs.precision(8);
pmesh->Print(mesh_ofs);
ofstream sol_ofs(sol_name.str().c_str());
sol_ofs.precision(8);
u.Save(sol_ofs);
}
// 15. Send the solution by socket to a GLVis server.
if (visualization)
{
char vishost[] = "localhost";
int visport = 19916;
// Make sure all ranks have sent their 'v' solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
vis_T.precision(8);
vis_Q.precision(8);
vis_errT.precision(8);
int Wx = 0, Wy = 0; // window position
int Ww = 350, Wh = 350; // window size
int offx = Ww+10;//, offy = Wh+45; // window offsets
miniapps::VisualizeField(vis_Q, vishost, visport,
Q_gf, "Heat Soruce", Wx, Wy, Ww, Wh);
Wx += offx;
miniapps::VisualizeField(vis_T, vishost, visport,
u, "Temperature", Wx, Wy, Ww, Wh);
Wx += offx;
miniapps::VisualizeField(vis_errT, vishost, visport,
u_error, "Error in T", Wx, Wy, Ww, Wh);
}
if (conv)
{
cout << "Number of Newton Iterations: " << it+1 << endl;
break;
}
}
// 16. Free the used memory.
// delete pcg;
// delete amg;
// delete a;
delete Q;
delete fespace;
if (order > 0) { delete fec; }
delete pmesh;
MPI_Finalize();
return 0;
}
ImplicitDiffOp::ImplicitDiffOp(ParFiniteElementSpace & H1_FESpace,
Coefficient & TBdr,
Array<int> & bdr_attr,
ChiCoef & chi,
dChiCoef & dchi,
Coefficient & heatSource)
: Operator(H1_FESpace.GetTrueVSize()),
first_(true),
ess_bdr_attr_(bdr_attr),
bdrCoef_(&TBdr),
chiCoef_(&chi),
dChiCoef_(&dchi),
QCoef_(&heatSource),
// dtChiCoef_(1.0, *chiCoef_),
T_(&H1_FESpace),
gradTCoef_(&T_),
// dtGradTCoef_(-1.0, gradTCoef_),
dChiGradTCoef_(*dChiCoef_, gradTCoef_),
s0chi_(&H1_FESpace),
a0_(&H1_FESpace),
// dTdt_(&H1_FESpace),
Q_(&H1_FESpace),
Qs_(&H1_FESpace),
rhs_(&H1_FESpace),
RHS_(H1_FESpace.GetTrueVSize()),
// RHS0_(0),
AInv_(NULL),
APrecond_(NULL)
{
H1_FESpace.GetEssentialTrueDofs(ess_bdr_attr_, ess_bdr_tdofs_);
s0chi_.AddDomainIntegrator(new DiffusionIntegrator(*chiCoef_));
a0_.AddDomainIntegrator(new DiffusionIntegrator(*chiCoef_));
//a0_.AddDomainIntegrator(new MixedScalarWeakDivergenceIntegrator(
// dChiGradTCoef_));
Qs_.AddDomainIntegrator(new DomainLFIntegrator(*QCoef_));
Qs_.Assemble();
Qs_.ParallelAssemble(RHS_);
}
ImplicitDiffOp::~ImplicitDiffOp()
{
delete AInv_;
delete APrecond_;
}
/*
void ImplicitDiffOp::SetState(ParGridFunction & T)
{
T_ = T;
if (first_)
{
s0chi_.Assemble();
s0chi_.Finalize();
ofstream ofsS0("s0_const_initial.mat");
s0chi_.SpMat().Print(ofsS0);
a0_.Assemble();
a0_.Finalize();
cout << "Assembling Q" << endl;
Qs_.Assemble();
Qs_.ParallelAssemble(RHS_);
cout << "Norm of Q: " << Qs_.Norml2() << endl;
}
first_ = false;
}
*/
void ImplicitDiffOp::Mult(const Vector &T, Vector &Q) const
{
T_.Distribute(T);
// add(T0_, dt_, dT_, T1_);
chiCoef_->SetTemp(T_);
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
s0chi_.Mult(T_, Q_);
Q_.ParallelAssemble(Q);
Q.SetSubVector(ess_bdr_tdofs_, 0.0);
}
Operator & ImplicitDiffOp::GetGradient(const Vector &T) const
{
T_.Distribute(T);
chiCoef_->SetTemp(T_);
dChiCoef_->SetTemp(T_);
gradTCoef_.SetGridFunction(&T_);
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
a0_.Update();
a0_.Assemble(0);
a0_.Finalize(0);
rhs_ = Qs_;
T_.ProjectBdrCoefficient(*bdrCoef_, ess_bdr_attr_);
a0_.FormLinearSystem(ess_bdr_tdofs_, T_, rhs_, A_, SOL_, RHS_);
return A_;
}
Solver & ImplicitDiffOp::GetGradientSolver() const
{
if (AInv_ == NULL)
{
/*
HypreSmoother *J_hypreSmoother = new HypreSmoother;
J_hypreSmoother->SetType(HypreSmoother::l1Jacobi);
J_hypreSmoother->SetPositiveDiagonal(true);
JPrecond_ = J_hypreSmoother;
GMRESSolver * AInv_gmres = NULL;
cout << "Building GMRES" << endl;
AInv_gmres = new GMRESSolver(T0_.ParFESpace()->GetComm());
AInv_gmres->SetRelTol(1e-12);
AInv_gmres->SetAbsTol(0.0);
AInv_gmres->SetMaxIter(20000);
AInv_gmres->SetPrintLevel(2);
AInv_gmres->SetPreconditioner(*JPrecond_);
AInv_ = AInv_gmres;
*/
HypreGMRES * AInv_gmres = NULL;
cout << "Building HypreGMRES" << endl;
AInv_gmres = new HypreGMRES(T_.ParFESpace()->GetComm());
AInv_gmres->SetTol(1e-12);
AInv_gmres->SetMaxIter(200);
AInv_gmres->SetPrintLevel(2);
if ( APrecond_ == NULL )
{
cout << "Building AMG" << endl;
APrecond_ = new HypreBoomerAMG();
APrecond_->SetPrintLevel(0);
AInv_gmres->SetPreconditioner(*APrecond_);
}
AInv_ = AInv_gmres;
}
return *AInv_;
}
+958
View File
@@ -0,0 +1,958 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
//
// -----------------------------------------------------
// Fourier Miniapp: Thermal Diffusion
// -----------------------------------------------------
//
// This miniapp solves a time dependent heat equation.
//
#include "fourier_solver.hpp"
#include <cassert>
#include <memory>
#include <iostream>
#include <fstream>
using namespace std;
using namespace mfem;
using namespace mfem::thermal;
void display_banner(ostream & os);
static int prob_ = 1;
static int gamma_ = 10;
static double alpha_ = NAN;
static double chi_max_ratio_ = 1.0;
static double chi_min_ratio_ = 1.0;
double QFunc(const Vector &x, double t)
{
switch (prob_)
{
case 1:
{
return 2.0 * M_PI * M_PI * sin(M_PI * x[0]) * sin(M_PI * x[1]);
}
case 2:
case 4:
{
double a = 0.4;
double b = 0.8;
double r = pow(x[0] / a, 2) + pow(x[1] / b, 2);
double e = exp(-0.25 * t * M_PI * M_PI / (a * b) );
if ( r == 0.0 )
return 0.25 * M_PI * M_PI *
( (1.0 - e) * ( pow(a, -2) + pow(b, -2) ) + e / (a * b));
return ( M_PI / r ) *
( 0.25 * M_PI * pow(a * b, -4) *
( pow(b * b * x[0],2) + pow(a * a * x[1], 2) +
(a - b) * (b * pow(b * x[0], 2) - a * pow(a*x[1],2)) * e) *
cos(0.5 * M_PI * sqrt(r)) +
0.5 * pow(a * b, -2) * (x * x) * (1.0 - e) *
sin(0.5 * M_PI * sqrt(r)) / sqrt(r)
);
}
case 3:
{
double cx = cos(M_PI * (x[0]-0.5));
double cy = cos(M_PI * (x[1]-0.5));
double c2x = cos(2.0 * M_PI * (x[0]-0.5));
double s2x = sin(2.0 * M_PI * (x[0]-0.5));
double c2y = cos(2.0 * M_PI * (x[1]-0.5));
double s2y = sin(2.0 * M_PI * (x[1]-0.5));
double c2a = cos(2.0 * alpha_);
double s2a = sin(2.0 * alpha_);
double ccg = 0.5 * M_PI * M_PI * gamma_ * pow(cx * cy, gamma_ - 2);
double perp = 1.0 * gamma_ * (c2x * c2y - 1.0) + c2x + c2y + 2.0;
double para = 0.5 * (gamma_ * (c2x * c2y - s2a * s2x * s2y - 1.0) +
(gamma_ - 1.0) * c2a * (c2x - c2y) +
c2x + c2y + 2.0);
return ccg * (1.0 * perp + (chi_max_ratio_ - 1.0) * para);
}
}
}
//static double chi_ratio_ = 1.0;
double TFunc(const Vector &x, double t)
{
switch (prob_)
{
case 1:
{
double e = exp(-2.0 * M_PI * M_PI * t);
return sin(M_PI * x[0]) * sin(M_PI * x[1]) * (1.0 - e);
}
case 2:
{
double a = 0.4;
double b = 0.8;
double r = pow(x[0] / a, 2) + pow(x[1] / b, 2);
double e = exp(-0.25 * t * M_PI * M_PI / (a * b) );
return cos(0.5 * M_PI * sqrt(r)) * (1.0 - e);
}
case 3:
return pow(sin(M_PI * x[0]) * sin(M_PI * x[1]), gamma_);
case 4:
{
double a = 0.4;
double b = 0.8;
double r = pow(x[0] / a, 2) + pow(x[1] / b, 2);
double rs = pow(x[0] - 0.5 * a, 2) + pow(x[1] - 0.5 * b, 2);
return cos(0.5 * M_PI * sqrt(r)) + 0.5 * exp(-400.0 * rs);
}
}
}
void dTFunc(const Vector &x, double t, Vector &dT)
{
dT.SetSize(x.Size());
dT = 0.0;
switch (prob_)
{
case 1:
{
double e = exp(-2.0 * M_PI * M_PI * t);
dT[0] = M_PI * cos(M_PI * x[0]) * sin(M_PI * x[1]);
dT[1] = M_PI * sin(M_PI * x[0]) * cos(M_PI * x[1]);
dT *= (1.0 - e);
}
break;
case 2:
{
double a = 0.4;
double b = 0.8;
double r = pow(x[0] / a, 2) + pow(x[1] / b, 2);
double r_2 = sqrt(r);
double sr = sin(0.5 * M_PI * r_2);
double e = exp(-0.25 * t * M_PI * M_PI / (a * b) );
dT[0] = -0.5 * M_PI * x[0] * sr / ( a * a * r_2 );
dT[1] = -0.5 * M_PI * x[1] * sr / ( b * b * r_2 );
dT *= (1.0 - e);
}
break;
case 3:
{
double cx = cos(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sx = sin(M_PI * x[0]);
double sy = sin(M_PI * x[1]);
// T = pow(sin(M_PI * x[0]) * sin(M_PI * x[1]), gamma_);
dT[0] = cx * sy;
dT[1] = sx * cy;
dT *= M_PI * gamma_ * pow(sx * sy, gamma_ - 1);
}
break;
case 4:
{
double a = 0.4;
double b = 0.8;
double r = pow(x[0] / a, 2) + pow(x[1] / b, 2);
double rs = pow(x[0] - 0.5 * a, 2) + pow(x[1] - 0.5 * b, 2);
double ers = exp(-400.0 * rs);
double r_2 = sqrt(r);
double sr = sin(0.5 * M_PI * r_2);
// T = cos(0.5 * M_PI * sqrt(r)) + 0.5 * exp(-400.0 * rs);
dT[0] = -0.5 * M_PI * x[0] * sr / ( a * a * r_2 );
dT[1] = -0.5 * M_PI * x[1] * sr / ( b * b * r_2 );
dT[0] -= 400.0 * (x[0] - 0.5 * a) * ers;
dT[1] -= 400.0 * (x[1] - 0.5 * b) * ers;
}
break;
}
}
void ChiFunc(const Vector &x, DenseMatrix &M)
{
M.SetSize(2);
switch (prob_)
{
case 1:
{
double cx = cos(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sx = sin(M_PI * x[0]);
double sy = sin(M_PI * x[1]);
double den = cx * cx * sy * sy + sx * sx * cy * cy;
M(0,0) = chi_max_ratio_ * sx * sx * cy * cy + sy * sy * cx * cx;
M(1,1) = chi_max_ratio_ * sy * sy * cx * cx + sx * sx * cy * cy;
M(0,1) = (1.0 - chi_max_ratio_) * cx * cy * sx * sy;
M(1,0) = M(0,1);
M *= 1.0 / den;
}
break;
case 2:
case 4:
{
double a = 0.4;
double b = 0.8;
double den = pow(b * b * x[0], 2) + pow(a * a * x[1], 2);
M(0,0) = chi_max_ratio_ * pow(a * a * x[1], 2) + pow(b * b * x[0], 2);
M(1,1) = chi_max_ratio_ * pow(b * b * x[0], 2) + pow(a * a * x[1], 2);
M(0,1) = (1.0 - chi_max_ratio_) * pow(a * b, 2) * x[0] * x[1];
M(1,0) = M(0,1);
M *= 1.0 / den;
}
break;
case 3:
{
double ca = cos(alpha_);
double sa = sin(alpha_);
M(0,0) = 1.0 + (chi_max_ratio_ - 1.0) * ca * ca;
M(1,1) = 1.0 + (chi_max_ratio_ - 1.0) * sa * sa;
M(0,1) = (chi_max_ratio_ - 1.0) * ca * sa;
M(1,0) = (chi_max_ratio_ - 1.0) * ca * sa;
}
break;
}
}
void bbTFunc(const Vector &x, DenseMatrix &M)
{
M.SetSize(2);
switch (prob_)
{
case 1:
{
double cx = cos(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sx = sin(M_PI * x[0]);
double sy = sin(M_PI * x[1]);
double den = cx * cx * sy * sy + sx * sx * cy * cy;
M(0,0) = sx * sx * cy * cy;
M(1,1) = sy * sy * cx * cx;
M(0,1) = -1.0 * cx * cy * sx * sy;
M(1,0) = M(0,1);
M *= 1.0 / den;
}
break;
case 2:
case 4:
{
double a = 0.4;
double b = 0.8;
double den = pow(b * b * x[0], 2) + pow(a * a * x[1], 2);
M(0,0) = pow(a * a * x[1], 2);
M(1,1) = pow(b * b * x[0], 2);
M(0,1) = -1.0 * pow(a * b, 2) * x[0] * x[1];
M(1,0) = M(0,1);
M *= 1.0 / den;
}
break;
case 3:
{
double ca = cos(alpha_);
double sa = sin(alpha_);
M(0,0) = ca * ca;
M(1,1) = sa * sa;
M(0,1) = ca * sa;
M(1,0) = ca * sa;
}
break;
}
}
class ChiGridFuncCoef : public MatrixCoefficient
{
private:
GridFunction * T_;
public:
ChiGridFuncCoef(GridFunction & T) : MatrixCoefficient(2), T_(&T) {}
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
void qFunc(const Vector &x, double t, Vector &q)
{
DenseMatrix Chi(x.Size());
Vector dT(x.Size());
dTFunc(x, t, dT);
ChiFunc(x, Chi);
Chi.Mult(dT, q);
q *= -1.0;
}
long int factorial(unsigned int n)
{
long int fact = 1;
for (unsigned int i=2; i<=n; i++)
{
fact *= i;
}
return fact;
}
// Returns the Gamma(n) function for a positive integer n
long int gamma(unsigned int n)
{
assert(n > 0);
return factorial(n-1);
}
// Returns Gamma(n+1/2) for a positive integer n
double gamma1_2(unsigned int n)
{
return sqrt(M_PI) * factorial(2*n) / (pow(4, n) * factorial(n));
}
double TNorm()
{
switch (prob_)
{
case 1:
return 0.5;
case 2:
return (gamma1_2((unsigned int)gamma_) /
gamma((unsigned int)gamma_+1)) / sqrt(M_PI);
}
}
double qPerpNorm()
{
switch (prob_)
{
case 1:
return M_PI * M_SQRT1_2 * chi_max_ratio_;
case 3:
return sqrt(M_PI * gamma_) * M_SQRT1_2 *
sqrt(gamma1_2((unsigned int)gamma_-1) *
gamma1_2((unsigned int)gamma_)) /
sqrt(gamma((unsigned int)gamma_) * gamma((unsigned int)gamma_+1));
}
}
double qParaNorm()
{
switch (prob_)
{
case 1:
return 0.0;
case 3:
return chi_max_ratio_ * qPerpNorm();
}
}
int main(int argc, char *argv[])
{
// 1. Initialize MPI.
MPI_Session mpi(argc, argv);
int myid = mpi.WorldRank();
// print the cool banner
if (mpi.Root()) { display_banner(cout); }
// 2. Parse command-line options.
int n = -1;
int order = 1;
int irOrder = -1;
int el_type = Element::QUADRILATERAL;
int ode_solver_type = 1;
int vis_steps = 1;
double dt = 0.5;
double t_final = 5.0;
double tol = 1e-4;
const char *basename = "Fourier";
const char *mesh_file = "";
bool zero_start = true;
bool static_cond = false;
bool gfprint = true;
bool visit = true;
bool visualization = true;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&n, "-n", "--num-elems-1d",
"Number of elements in x and y directions. "
"Total number of elements is n^2.");
args.AddOption(&prob_, "-p", "--problem",
"Specify problem type: 1 - Square, 2 - Ellipse, 3 - van Es.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree) or -1 for"
" isoparametric space.");
args.AddOption(&irOrder, "-iro", "--int-rule-order",
"Integration Rule Order.");
args.AddOption(&alpha_, "-alpha", "--constant-angle",
"Angle for constant B field (in degrees)");
args.AddOption(&gamma_, "-gamma", "--exponent",
"Exponent used in problem 2");
args.AddOption(&chi_max_ratio_, "-chi-max", "--chi-max-ratio",
"Ratio of chi_max_parallel/chi_perp.");
args.AddOption(&chi_min_ratio_, "-chi-min", "--chi-min-ratio",
"Ratio of chi_min_parallel/chi_perp.");
args.AddOption(&dt, "-dt", "--time-step",
"Time step.");
args.AddOption(&t_final, "-tf", "--final-time",
"Final Time.");
args.AddOption(&tol, "-tol", "--tolerance",
"Tolerance used to determine convergence to steady state.");
args.AddOption(&el_type, "-e", "--element-type",
"Element type: 2-Triangle, 3-Quadrilateral.");
args.AddOption(&ode_solver_type, "-s", "--ode-solver",
"ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3\n\t."
"\t 22 - Mid-Point, 23 - SDIRK23, 34 - SDIRK34.");
args.AddOption(&zero_start, "-z", "--zero-start", "-no-z",
"--no-zero-start",
"Initial guess of zero or exact solution.");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&gfprint, "-print", "--print","-no-print","--no-print",
"Print results (grid functions) to disk.");
args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
"Enable or disable VisIt visualization.");
args.AddOption(&vis_steps, "-vs", "--visualization-steps",
"Visualize every n-th timestep.");
args.AddOption(&basename, "-k", "--outputfilename",
"Name of the visit dump files");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
if (irOrder < 0)
{
irOrder = std::max(4, 2 * order - 2);
}
if (isnan(alpha_))
{
alpha_ = 0.0;
}
else
{
alpha_ *= M_PI / 180.0;
}
// 3. Construct a (serial) mesh of the given size on all processors. We
// can handle triangular and quadrilateral surface meshes with the
// same code.
Mesh *mesh = (n > 0) ?
new Mesh(n, n, (Element::Type)el_type, 1) :
new Mesh(mesh_file, 1, 1);
int dim = mesh->Dimension();
// 4. This step is no longer needed
// 5. Define a parallel mesh by a partitioning of the serial mesh. Once the
// parallel mesh is defined, the serial mesh can be deleted.
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
delete mesh;
// 7. Determine the list of true (i.e. parallel conforming) essential
// boundary dofs. In this example, the boundary conditions are defined
// by marking all the boundary attributes from the mesh as essential
// (Dirichlet) and converting them to a list of true dofs.
Array<int> ess_tdof_list;
Array<int> ess_bdr(0);
if (pmesh->bdr_attributes.Size())
{
ess_bdr.SetSize(pmesh->bdr_attributes.Max());
ess_bdr = 1;
}
// The following is required for mesh refinement
// mesh->EnsureNCMesh();
// 6. Define the ODE solver used for time integration. Several implicit
// methods are available, including singly diagonal implicit Runge-Kutta
// (SDIRK).
ODESolver *ode_solver;
switch (ode_solver_type)
{
// Implicit L-stable methods
case 1: ode_solver = new BackwardEulerSolver; break;
case 2: ode_solver = new SDIRK23Solver(2); break;
case 3: ode_solver = new SDIRK33Solver; break;
// Implicit A-stable methods (not L-stable)
case 22: ode_solver = new ImplicitMidpointSolver; break;
case 23: ode_solver = new SDIRK23Solver; break;
case 34: ode_solver = new SDIRK34Solver; break;
default:
if (mpi.Root())
{
cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
}
delete mesh;
return 3;
}
// 12. Define the parallel finite element spaces. We use:
//
// H(curl) for electric field,
// H(div) for magnetic flux,
// H(div) for thermal flux,
// H(grad)/H1 for electrostatic potential,
// L2 for temperature
// L2 contains discontinuous "cell-center" finite elements, type 2 is
// "positive"
L2_FECollection L2FEC0(0, dim);
L2_FECollection L2FEC(order-1, dim);
// RT contains Raviart-Thomas "face-centered" vector finite elements with
// continuous normal component.
RT_FECollection HDivFEC(order-1, dim);
// ND contains Nedelec "edge-centered" vector finite elements with
// continuous tangential component.
ND_FECollection HCurlFEC(order, dim);
// H1 contains continuous "node-centered" Lagrange finite elements.
H1_FECollection HGradFEC(order, dim);
ParFiniteElementSpace L2FESpace0(pmesh, &L2FEC0);
ParFiniteElementSpace L2FESpace(pmesh, &L2FEC);
ParFiniteElementSpace HDivFESpace(pmesh, &HDivFEC);
ParFiniteElementSpace HCurlFESpace(pmesh, &HCurlFEC);
ParFiniteElementSpace HGradFESpace(pmesh, &HGradFEC);
// The terminology is TrueVSize is the unique (non-redundant) number of dofs
// HYPRE_Int glob_size_l2 = L2FESpace.GlobalTrueVSize();
// HYPRE_Int glob_size_rt = HDivFESpace.GlobalTrueVSize();
HYPRE_Int glob_size_h1 = HGradFESpace.GlobalTrueVSize();
if (mpi.Root())
{
cout << "Number of Temperature unknowns: " << glob_size_h1 << endl;
}
// int Vsize_l2 = L2FESpace.GetVSize();
// int Vsize_rt = HDivFESpace.GetVSize();
// int Vsize_h1 = HGradFESpace.GetVSize();
// grid functions E, B, T, F, P, and w which is the Joule heating
ParGridFunction q(&HCurlFESpace);
ParGridFunction qPara(&HCurlFESpace);
ParGridFunction qPerp(&HCurlFESpace);
ParGridFunction Q(&L2FESpace);
ParGridFunction T1(&HGradFESpace);
ParGridFunction T0(&HGradFESpace);
ParGridFunction dT(&HGradFESpace);
ParGridFunction errorq(&L2FESpace0);
ParGridFunction errorqPara(&L2FESpace0);
ParGridFunction errorqPerp(&L2FESpace0);
ParGridFunction errorT(&L2FESpace0);
T0 = 0.0;
T1 = 0.0;
dT = 1.0;
// 13. Get the boundary conditions, set up the exact solution grid functions
// These VectorCoefficients have an Eval function. Note that e_exact and
// b_exact in this case are exact analytical solutions, taking a 3-vector
// point as input and returning a 3-vector field
FunctionCoefficient TCoef(TFunc);
VectorFunctionCoefficient qCoef(2, qFunc);
Vector zeroVec(dim); zeroVec = 0.0;
ConstantCoefficient zeroCoef(0.0);
VectorConstantCoefficient zeroVecCoef(zeroVec);
IdentityMatrixCoefficient ICoef(2);
MatrixFunctionCoefficient bbTCoef(2, bbTFunc);
MatrixSumCoefficient ImbbTCoef(bbTCoef, ICoef, -1.0);
MatVecCoefficient qParaCoef(bbTCoef, qCoef);
MatVecCoefficient qPerpCoef(ImbbTCoef, qCoef);
ConstantCoefficient SpecificHeatCoef(1.0);
MatrixFunctionCoefficient ConductionCoef(2, ChiFunc);
FunctionCoefficient HeatSourceCoef(QFunc);
Q.ProjectCoefficient(HeatSourceCoef);
if (!zero_start)
{
T1.ProjectCoefficient(TCoef);
q.ProjectCoefficient(qCoef);
}
T1.GridFunction::ComputeElementL2Errors(TCoef, errorT);
q.GridFunction::ComputeElementL2Errors(qCoef, errorq);
qPara.GridFunction::ComputeElementL2Errors(qParaCoef, errorqPara);
qPerp.GridFunction::ComputeElementL2Errors(qPerpCoef, errorqPerp);
ParBilinearForm m1(&HCurlFESpace);
m1.AddDomainIntegrator(new VectorFEMassIntegrator);
m1.Assemble();
ParMixedBilinearForm gPara(&HGradFESpace, &HCurlFESpace);
gPara.AddDomainIntegrator(new MixedVectorGradientIntegrator(bbTCoef));
gPara.Assemble();
ParMixedBilinearForm gPerp(&HGradFESpace, &HCurlFESpace);
gPerp.AddDomainIntegrator(new MixedVectorGradientIntegrator(ImbbTCoef));
gPerp.Assemble();
HypreParMatrix M1C;
Vector RHS1(HCurlFESpace.GetTrueVSize()), X1(HCurlFESpace.GetTrueVSize());
Array<int> ess_tdof_list_q(0);
// Array<int> ess_bdr_q;
// HCurlFESpace.GetEssentialTrueDofs(ess_bdr_q, ess_tdof_list_q);
m1.FormSystemMatrix(ess_tdof_list_q, M1C);
HypreDiagScale Precond(M1C);
HyprePCG M1Inv(M1C);
M1Inv.SetTol(1e-12);
M1Inv.SetMaxIter(200);
M1Inv.SetPrintLevel(0);
M1Inv.SetPreconditioner(Precond);
// 14. Initialize the Diffusion operator, the GLVis visualization and print
// the initial energies.
ThermalDiffusionOperator oper(HGradFESpace,
zeroCoef, ess_bdr,
SpecificHeatCoef, false,
ConductionCoef, false,
HeatSourceCoef, false);
// This function initializes all the fields to zero or some provided IC
// oper.Init(F);
socketstream vis_Q;
socketstream vis_q, vis_errq;
socketstream vis_qPara, vis_errqPara;
socketstream vis_qPerp, vis_errqPerp;
socketstream vis_T, vis_errT;
char vishost[] = "localhost";
int visport = 19916;
if (visualization)
{
// Make sure all ranks have sent their 'v' solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
vis_Q.precision(8);
vis_T.precision(8);
vis_errT.precision(8);
vis_q.precision(8);
vis_errq.precision(8);
vis_qPara.precision(8);
vis_errqPara.precision(8);
vis_qPerp.precision(8);
vis_errqPerp.precision(8);
int Wx = 0, Wy = 0; // window position
int Ww = 280, Wh = 280; // window size
int offx = Ww+10, offy = Wh+45; // window offsets
miniapps::VisualizeField(vis_Q, vishost, visport,
Q, "Heat Source", Wx, Wy, Ww, Wh);
Wy += offy;
// miniapps::VisualizeField(vis_U, vishost, visport,
// U1, "Energy", Wx, Wy, Ww, Wh);
Wx += offx;
Wy -= offy;
miniapps::VisualizeField(vis_T, vishost, visport,
T1, "Temperature", Wx, Wy, Ww, Wh);
Wy += offy;
miniapps::VisualizeField(vis_errT, vishost, visport,
errorT, "Error in T", Wx, Wy, Ww, Wh);
Wx += offx;
Wy -= offy;
miniapps::VisualizeField(vis_q, vishost, visport,
q, "Heat Flux", Wx, Wy, Ww, Wh);
Wy += offy;
miniapps::VisualizeField(vis_errq, vishost, visport,
errorq, "Error in q", Wx, Wy, Ww, Wh);
Wx += offx;
Wy -= offy;
miniapps::VisualizeField(vis_qPara, vishost, visport,
qPara, "Parallel Heat Flux", Wx, Wy, Ww, Wh);
Wy += offy;
miniapps::VisualizeField(vis_errqPara, vishost, visport,
errorqPara, "Error in q para", Wx, Wy, Ww, Wh);
Wx += offx;
Wy -= offy;
miniapps::VisualizeField(vis_qPerp, vishost, visport,
qPerp, "Perpendicular Heat Flux", Wx, Wy, Ww, Wh);
Wy += offy;
miniapps::VisualizeField(vis_errqPerp, vishost, visport,
errorqPerp, "Error in q perp", Wx, Wy, Ww, Wh);
}
// VisIt visualization
VisItDataCollection visit_dc(basename, pmesh);
if ( visit )
{
visit_dc.RegisterField("Q", &Q);
visit_dc.RegisterField("q", &q);
visit_dc.RegisterField("qPara", &qPara);
visit_dc.RegisterField("qPerp", &qPerp);
visit_dc.RegisterField("T", &T1);
visit_dc.RegisterField("L2 Error T", &errorT);
visit_dc.RegisterField("L2 Error q", &errorq);
visit_dc.RegisterField("L2 Error q para", &errorqPara);
visit_dc.RegisterField("L2 Error q perp", &errorqPerp);
visit_dc.SetCycle(0);
visit_dc.SetTime(0.0);
visit_dc.Save();
}
// 15. Perform time-integration (looping over the time iterations, ti, with a
// time-step dt). The object oper is the MagneticDiffusionOperator which
// has a Mult() method and an ImplicitSolve() method which are used by
// the time integrators.
ode_solver->Init(oper);
double t = 0.0;
bool last_step = false;
for (int ti = 1; !last_step; ti++)
{
if (t + dt >= t_final - dt/2)
{
if (myid == 0)
{
cout << "Final Time Reached" << endl;
}
last_step = true;
}
// F is the vector of dofs, t is the current time, and dt is the time step
// to advance.
T0 = T1;
ode_solver->Step(T1, t, dt);
add(1.0, T1, -1.0, T0, dT);
double maxT = T1.ComputeMaxError(zeroCoef);
double maxDiff = dT.ComputeMaxError(zeroCoef);
if ( !last_step )
{
if ( maxT == 0.0 )
{
last_step = (maxDiff < tol) ? true:false;
}
else if ( maxDiff/maxT < tol )
{
last_step = true;
}
if (last_step && myid == 0)
{
cout << "Converged to Steady State" << endl;
}
}
/*
if (debug == 1)
{
oper.Debug(basename,t);
}
*/
gPara.Mult(T1, qPara);
gPerp.Mult(T1, qPerp);
qPara.ParallelAssemble(RHS1);
X1 = 0.0;
M1Inv.Mult(RHS1, X1);
qPara.Distribute(X1);
qPara *= -chi_max_ratio_;
qPerp.ParallelAssemble(RHS1);
X1 = 0.0;
M1Inv.Mult(RHS1, X1);
qPerp.Distribute(X1);
qPerp *= -1.0;
q = qPara;
q += qPerp;
if (gfprint)
{
ostringstream q_name, T_name, mesh_name;
q_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "q." << setfill('0') << setw(6) << myid;
T_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "T." << setfill('0') << setw(6) << myid;
mesh_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "mesh." << setfill('0') << setw(6) << myid;
ofstream mesh_ofs(mesh_name.str().c_str());
mesh_ofs.precision(8);
pmesh->Print(mesh_ofs);
mesh_ofs.close();
ofstream q_ofs(q_name.str().c_str());
q_ofs.precision(8);
q.Save(q_ofs);
q_ofs.close();
ofstream T_ofs(T_name.str().c_str());
T_ofs.precision(8);
T1.Save(T_ofs);
T_ofs.close();
}
if (last_step || (ti % vis_steps) == 0)
{
// Make sure all ranks have sent their 'v' solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
if (visualization)
{
int Wx = 0, Wy = 0; // window position
int Ww = 350, Wh = 350; // window size
int offx = Ww+10, offy = Wh+45; // window offsets
miniapps::VisualizeField(vis_q, vishost, visport,
q, "Heat Flux", Wx, Wy, Ww, Wh);
miniapps::VisualizeField(vis_qPara, vishost, visport,
qPara, "Parallel Heat Flux",
Wx, Wy, Ww, Wh);
miniapps::VisualizeField(vis_qPerp, vishost, visport,
qPerp, "Perpendicular Heat Flux",
Wx, Wy, Ww, Wh);
// Wx += offx;
// miniapps::VisualizeField(vis_U, vishost, visport,
// U1, "Energy", Wx, Wy, Ww, Wh);
// Wx -= offx;
// Wy += offy;
miniapps::VisualizeField(vis_T, vishost, visport,
T1, "Temperature", Wx, Wy, Ww, Wh);
// Wx += offx;
miniapps::VisualizeField(vis_errT, vishost, visport,
errorT, "Error in T", Wx, Wy, Ww, Wh);
// Wx += offx;
miniapps::VisualizeField(vis_errq, vishost, visport,
errorq, "Error in q", Wx, Wy, Ww, Wh);
miniapps::VisualizeField(vis_errqPara, vishost, visport,
errorqPara, "Error in q para",
Wx, Wy, Ww, Wh);
miniapps::VisualizeField(vis_errqPerp, vishost, visport,
errorqPerp, "Error in q perp",
Wx, Wy, Ww, Wh);
}
if (visit)
{
visit_dc.SetCycle(ti);
visit_dc.SetTime(t);
visit_dc.Save();
}
}
}
if (visualization)
{
vis_Q.close();
vis_q.close();
vis_T.close();
vis_errT.close();
vis_errq.close();
vis_errqPara.close();
vis_errqPerp.close();
}
double loc_T_max = T1.Normlinf();
double T_max = -1.0;
MPI_Allreduce(&loc_T_max, &T_max, 1, MPI_DOUBLE, MPI_MAX,
MPI_COMM_WORLD);
double err1 = T1.ComputeL2Error(TCoef);
if (myid == 0)
{
cout << "L2 Error of Solution: " << err1 << endl;
cout << "Maximum Temperature: " << T_max << endl;
cout << "| chi_eff - 1 | = " << fabs(1.0/T_max - 1) << endl;
}
// 16. Free the used memory.
delete ode_solver;
delete pmesh;
return 0;
}
void display_banner(ostream & os)
{
os << "___________ .__ " << endl
<< "\\_ _____/___ __ _________|__| ___________ " << endl
<< " | __)/ _ \\| | \\_ __ \\ |/ __ \\_ __ \\" << endl
<< " | | ( <_> ) | /| | \\/ \\ ___/| | \\/" << endl
<< " \\__ | \\____/|____/ |__| |__|\\___ >__| " << endl
<< " \\/ \\/ " << endl
<< flush;
}
File diff suppressed because it is too large Load Diff
+441
View File
@@ -0,0 +1,441 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#include "fourier_flux_solver.hpp"
#ifdef MFEM_USE_MPI
using namespace std;
namespace mfem
{
using namespace miniapps;
namespace thermal
{
ThermalDiffusionFluxOperator::ThermalDiffusionFluxOperator(
ParMesh & pmesh,
ParFiniteElementSpace &HDiv_FES,
ParFiniteElementSpace &L2_FES,
VectorCoefficient & dqdtBdr,
Array<int> & bdr_attr,
Coefficient & c, bool td_c,
Coefficient & k, bool td_k,
Coefficient & Q, bool td_Q)
: TimeDependentOperator(HDiv_FES.GetVSize() + L2_FES.GetVSize(), 0.0),
init_(false), //initA_(false), initAInv_(false),
dim_(pmesh.Dimension()),
multCount_(0), solveCount_(0),
HDiv_FESpace_(&HDiv_FES),
L2_FESpace_(&L2_FES),
mK_(NULL), sC_(NULL), dC_(NULL), a_(NULL), Div_(NULL),
dqdt_gf_(NULL), Qs_(NULL),
MKInv_(NULL), MKDiag_(NULL),
AInv_(NULL), APrecond_(NULL),
// rhs_(NULL),
bdr_attr_(&bdr_attr), ess_bdr_tdofs_(0), dqdtBdrCoef_(&dqdtBdr),
tdQ_(td_Q), tdC_(td_c), tdK_(td_k),
QCoef_(&Q), CCoef_(&c), kCoef_(&k), KCoef_(NULL),
// CInvCoef_(NULL), kInvCoef_(NULL), KInvCoef_(NULL)
CInvCoef_(new InverseCoefficient(c)),
kInvCoef_(new InverseCoefficient(k)), KInvCoef_(NULL),
dtCInvCoef_(NULL)
{
this->init();
}
ThermalDiffusionFluxOperator::ThermalDiffusionFluxOperator(
ParMesh & pmesh,
ParFiniteElementSpace &HDiv_FES,
ParFiniteElementSpace &L2_FES,
VectorCoefficient & dqdtBdr,
Array<int> & bdr_attr,
Coefficient & c, bool td_c,
MatrixCoefficient & K, bool td_k,
Coefficient & Q, bool td_Q)
: TimeDependentOperator(HDiv_FES.GetVSize() + L2_FES.GetVSize(), 0.0),
init_(false),
dim_(pmesh.Dimension()),
multCount_(0), solveCount_(0),
HDiv_FESpace_(&HDiv_FES),
L2_FESpace_(&L2_FES),
mK_(NULL), sC_(NULL), dC_(NULL), a_(NULL), Div_(NULL),
dqdt_gf_(NULL), Qs_(NULL),
MKInv_(NULL), MKDiag_(NULL),
AInv_(NULL), APrecond_(NULL),
// rhs_(NULL),
bdr_attr_(&bdr_attr), ess_bdr_tdofs_(0), dqdtBdrCoef_(&dqdtBdr),
tdQ_(td_Q), tdC_(td_c), tdK_(td_k),
QCoef_(&Q), CCoef_(&c), kCoef_(NULL), KCoef_(&K),
CInvCoef_(new InverseCoefficient(c)),
kInvCoef_(NULL),
KInvCoef_(new MatrixInverseCoefficient(K)),
dtCInvCoef_(NULL)
{
this->init();
}
ThermalDiffusionFluxOperator::~ThermalDiffusionFluxOperator()
{
delete CInvCoef_;
delete kInvCoef_;
delete KInvCoef_;
delete dtCInvCoef_;
delete Div_;
delete dC_;
delete a_;
delete mK_;
delete sC_;
delete dqdt_gf_;
delete Qs_;
delete MKInv_;
delete MKDiag_;
delete AInv_;
delete APrecond_;
}
void
ThermalDiffusionFluxOperator::init()
{
if ( init_ ) { return; }
if ( mK_ == NULL )
{
mK_ = new ParBilinearForm(HDiv_FESpace_);
if ( kCoef_ != NULL )
{
mK_->AddDomainIntegrator(new VectorFEMassIntegrator(*kInvCoef_));
}
else
{
mK_->AddDomainIntegrator(new VectorFEMassIntegrator(*KInvCoef_));
}
mK_->Assemble();
}
if ( sC_ == NULL )
{
sC_ = new ParBilinearForm(HDiv_FESpace_);
sC_->AddDomainIntegrator(new DivDivIntegrator(*CInvCoef_));
sC_->Assemble();
}
if ( dC_ == NULL )
{
dC_ = new ParMixedBilinearForm(L2_FESpace_, HDiv_FESpace_);
dC_->AddDomainIntegrator(
new MixedScalarWeakGradientIntegrator(*CInvCoef_));
dC_->Assemble();
}
if ( dqdt_gf_ == NULL )
{
dqdt_gf_ = new ParGridFunction(HDiv_FESpace_);
}
if ( Qs_ == NULL && QCoef_ != NULL )
{
Qs_ = new ParGridFunction(L2_FESpace_);
Qs_->ProjectCoefficient(*QCoef_);
}
Div_ = new ParDiscreteDivOperator(HDiv_FESpace_, L2_FESpace_);
Div_->Assemble();
Div_->Finalize();
rhs_.SetSize(HDiv_FESpace_->GetVSize());
dQs_.SetSize(HDiv_FESpace_->GetVSize());
tmp_.SetSize(L2_FESpace_->GetVSize());
HDiv_FESpace_->GetEssentialTrueDofs(*bdr_attr_, ess_bdr_tdofs_);
init_ = true;
}
void
ThermalDiffusionFluxOperator::SetTime(const double time)
{
this->TimeDependentOperator::SetTime(time);
dqdtBdrCoef_->SetTime(t);
if ( tdQ_ )
{
QCoef_->SetTime(t);
Qs_->ProjectCoefficient(*QCoef_);
}
if ( tdC_ )
{
// CCoef_->SetTime(t);
// CInvCoef_->SetTime(t);
dtCInvCoef_->SetTime(t);
sC_->Assemble();
}
if ( tdK_ )
{
if ( kCoef_ != NULL ) { kCoef_->SetTime(t); kInvCoef_->SetTime(t); }
if ( KCoef_ != NULL ) { KCoef_->SetTime(t); KInvCoef_->SetTime(t); }
mK_->Assemble();
}
if ( ( tdC_ || tdK_ ) && a_ != NULL )
{
a_->Assemble();
}
newTime_ = true;
}
/*
void
ThermalDiffusionFluxOperator::SetHeatSource(Coefficient & Q, bool time_dep)
{
if ( ownsQ_ )
{
delete QCoef_;
}
tdQ_ = time_dep;
QCoef_ = &Q;
}
void
ThermalDiffusionFluxOperator::SetConductivityCoefficient(Coefficient & k,
bool time_dep)
{
if ( ownsK_ )
{
delete kCoef_;
delete KCoef_;
}
tdK_ = time_dep;
kCoef_ = &k;
KCoef_ = NULL;
}
void
ThermalDiffusionFluxOperator::SetConductivityCoefficient(MatrixCoefficient & K,
bool time_dep)
{
if ( ownsK_ )
{
delete kCoef_;
delete KCoef_;
}
tdK_ = time_dep;
kCoef_ = NULL;
KCoef_ = &K;
}
void
ThermalDiffusionFluxOperator::SetSpecificHeatCoefficient(Coefficient & c,
bool time_dep)
{
if ( ownsC_ )
{
delete CCoef_;
}
tdC_ = time_dep;
CCoef_ = &c;
}
*/
void
ThermalDiffusionFluxOperator::initMult() const
{
if ( tdC_ || MKInv_ == NULL || MKDiag_ == NULL )
{
if ( MKInv_ == NULL )
{
MKInv_ = new HyprePCG(MK_);
MKInv_->SetTol(1e-12);
MKInv_->SetMaxIter(200);
MKInv_->SetPrintLevel(0);
}
else
{
MKInv_->SetOperator(MK_);
}
if ( MKDiag_ == NULL )
{
MKDiag_ = new HypreDiagScale(MK_);
MKInv_->SetPreconditioner(*MKDiag_);
}
else
{
MKDiag_->SetOperator(MK_);
}
}
}
void
ThermalDiffusionFluxOperator::Mult(const Vector &y, Vector &dy_dt) const
{
cout << "Entering Mult" << endl;
dy_dt = 0.0;
q_.MakeRef(const_cast<ParFiniteElementSpace*>(HDiv_FESpace_),
const_cast<Vector&>(y), 0);
u_.MakeRef(const_cast<ParFiniteElementSpace*>(L2_FESpace_),
const_cast<Vector&>(y), HDiv_FESpace_->GetVSize());
dqdt_.MakeRef(HDiv_FESpace_, dy_dt, 0);
dudt_.MakeRef(L2_FESpace_, dy_dt, HDiv_FESpace_->GetVSize());
sC_->Mult(q_, rhs_);
dC_->Mult(*Qs_, dQs_);
rhs_ += dQs_;
rhs_.Neg();
dqdt_gf_->ProjectBdrCoefficientNormal(*dqdtBdrCoef_, *bdr_attr_);
mK_->FormLinearSystem(ess_bdr_tdofs_, *dqdt_gf_, rhs_, MK_, X_, RHS_);
this->initMult();
MKInv_->Mult(RHS_, X_);
mK_->RecoverFEMSolution(X_, rhs_, dqdt_);
Div_->Mult(q_, dudt_);
dudt_ *= -1.0;
dudt_ += *Qs_;
multCount_++;
cout << "Leaving Mult" << endl;
}
void
ThermalDiffusionFluxOperator::initA(double dt)
{
if ( CInvCoef_ != NULL )
{
dtCInvCoef_ = new ScaledCoefficient(dt, *CInvCoef_);
}
if ( a_ == NULL)
{
a_ = new ParBilinearForm(HDiv_FESpace_);
if ( kInvCoef_ != NULL)
{
a_->AddDomainIntegrator(new VectorFEMassIntegrator(*kInvCoef_));
}
else
{
a_->AddDomainIntegrator(new VectorFEMassIntegrator(*KInvCoef_));
}
a_->AddDomainIntegrator(new DivDivIntegrator(*dtCInvCoef_));
a_->Assemble();
}
else if ( tdK_ )
{
a_->Update();
a_->Assemble();
}
}
void
ThermalDiffusionFluxOperator::initImplicitSolve()
{
if ( tdC_ || tdK_ || AInv_ == NULL || APrecond_ == NULL )
{
delete AInv_;
AInv_ = new HyprePCG(A_);
AInv_->SetTol(1e-12);
AInv_->SetMaxIter(200);
AInv_->SetPrintLevel(0);
delete APrecond_;
APrecond_ = (dim_==2) ?
(HypreSolver*)(new HypreAMS(A_, HDiv_FESpace_)):
(HypreSolver*)(new HypreADS(A_, HDiv_FESpace_));
if ( dim_ == 2 )
{
dynamic_cast<HypreAMS*>(APrecond_)->SetPrintLevel(0);
}
else
{
dynamic_cast<HypreADS*>(APrecond_)->SetPrintLevel(0);
}
AInv_->SetPreconditioner(*APrecond_);
}
}
void
ThermalDiffusionFluxOperator::ImplicitSolve(const double dt,
const Vector &y, Vector &dy_dt)
{
dy_dt = 0.0;
q_.MakeRef(const_cast<ParFiniteElementSpace*>(HDiv_FESpace_),
const_cast<Vector&>(y), 0);
u_.MakeRef(const_cast<ParFiniteElementSpace*>(L2_FESpace_),
const_cast<Vector&>(y), HDiv_FESpace_->GetVSize());
dqdt_.MakeRef(HDiv_FESpace_, dy_dt, 0);
dudt_.MakeRef(L2_FESpace_, dy_dt, HDiv_FESpace_->GetVSize());
// cout << "sC size: " << sC_->Width() << ", q_ size: " << q_.Size() << ", rhs_ size: " << rhs_.Size() << endl;
sC_->Mult(q_, rhs_);
dC_->Mult(*Qs_, dQs_);
rhs_ += dQs_;
rhs_ *= -1.0;
// dqdt_gf_->ProjectBdrCoefficientNormal(*dqdtBdrCoef_, *bdr_attr_);
dqdt_.ProjectBdrCoefficientNormal(*dqdtBdrCoef_, *bdr_attr_);
this->initA(dt);
// a_->FormLinearSystem(ess_bdr_tdofs_, *dqdt_gf_, rhs_, A_, X_, RHS_);
a_->FormLinearSystem(ess_bdr_tdofs_, dqdt_, rhs_, A_, X_, RHS_);
this->initImplicitSolve();
AInv_->Mult(RHS_, X_);
a_->RecoverFEMSolution(X_, rhs_, dqdt_);
Div_->Mult(q_, dudt_);
Div_->Mult(dqdt_, tmp_);
tmp_ *= dt;
dudt_ += tmp_;
dudt_ *= -1.0;
dudt_ += *Qs_;
solveCount_++;
}
} // namespace thermal
void
MatrixInverseCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K.Invert();
}
void
ScaledMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K *= a_;
}
} // namespace mfem
#endif // MFEM_USE_MPI
+269
View File
@@ -0,0 +1,269 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#ifndef MFEM_FOURIER_FLUX_SOLVER
#define MFEM_FOURIER_FLUX_SOLVER
#include "../common/pfem_extras.hpp"
#ifdef MFEM_USE_MPI
#include <memory>
#include <iostream>
#include <fstream>
namespace mfem
{
namespace thermal
{
/**
The thermal diffusion equation can be written:
dcT/dt = Div (chi Grad T) + Q_s
We would like to rewrite this using the flux formulation which solves for
the heat flux vector q. The primary equations are:
q = chi Grad T
u = c T
du/dt + Div q = Q_s
Which lead to:
dq/dt = chi Grad (c^{-1} Div q) - Grad(c^{-1} Q_s)
where
T is the temperature.
q is the heat flux
u is the thermal energy density
Div is the divergence operator,
Grad is the gradient operator,
chi is the thermal conductivity,
c is the heat capacity,
Q_s is the heat source
Class ThermalDiffusionFluxOperator represents the right-hand side of
the above system of ODEs.
f(t, T) = -M_0(c)^{-1}(S_0(chi)T - M_0 Q_s)
where
M_0(c) is an H_1 mass matrix
S_0(sigma) is the diffusion operator
The implicit solve method will solve
(M_0(c)+dt S_0(sigma))k = -S_0(sigma)T + M_0 Q_s
*/
class ThermalDiffusionFluxOperator : public TimeDependentOperator
{
public:
ThermalDiffusionFluxOperator(ParMesh & pmesh,
ParFiniteElementSpace &HDiv_FES,
ParFiniteElementSpace &L2_FES,
VectorCoefficient & dqdtBdr,
Array<int> & bdr_attr,
Coefficient & c, bool td_c,
Coefficient & k, bool td_k,
Coefficient & Q, bool td_Q);
ThermalDiffusionFluxOperator(ParMesh & pmesh,
ParFiniteElementSpace &HDiv_FES,
ParFiniteElementSpace &L2_FES,
VectorCoefficient & dqdtBdr,
Array<int> & bdr_attr,
Coefficient & c, bool td_c,
MatrixCoefficient & K, bool td_k,
Coefficient & Q, bool td_Q);
void SetTime(const double time);
/*
void SetHeatSource(Coefficient & Q, bool time_dep = false);
void SetConductivityCoefficient(Coefficient & k,
bool time_dep = false);
void SetConductivityCoefficient(MatrixCoefficient & K,
bool time_dep = false);
void SetSpecificHeatCoefficient(
bool time_dep = false);
*/
/** @brief Perform the action of the operator: @a q = f(@a y, t), where
q solves the algebraic equation F(@a y, q, t) = G(@a y, t) and t is the
current time. */
virtual void Mult(const Vector &y, Vector &q) const;
/** @brief Solve the equation: @a q = f(@a y + @a dt @a q, t), for the
unknown @a q at the current time t.
For general F and G, the equation for @a q becomes:
F(@a y + @a dt @a q, @a q, t) = G(@a y + @a dt @a q, t).
The input vector @a y corresponds to time index (or cycle) n, while the
currently set time, #t, and the result vector @a q correspond to time
index n+1. The time step @a dt corresponds to the time interval between
cycles n and n+1.
This method allows for the abstract implementation of some time
integration methods, including diagonal implicit Runge-Kutta (DIRK)
methods and the backward Euler method in particular.
If not re-implemented, this method simply generates an error. */
virtual void ImplicitSolve(const double dt, const Vector &y, Vector &q);
virtual ~ThermalDiffusionFluxOperator();
private:
void init();
void initMult() const;
void initA(double dt);
void initImplicitSolve();
bool init_;
// bool initA_;
// bool initAInv_;
bool newTime_;
int dim_;
mutable int multCount_;
int solveCount_;
ParFiniteElementSpace * HDiv_FESpace_;
ParFiniteElementSpace * L2_FESpace_;
ParBilinearForm * mK_;
ParBilinearForm * sC_;
ParMixedBilinearForm * dC_;
ParBilinearForm * a_;
ParDiscreteLinearOperator * Div_;
ParGridFunction * dqdt_gf_;
ParGridFunction * Qs_;
mutable HypreParMatrix MK_;
mutable HyprePCG * MKInv_;
mutable HypreDiagScale * MKDiag_;
HypreParMatrix A_;
HyprePCG * AInv_;
HypreSolver * APrecond_;
// HypreParVector * T_;
mutable ParGridFunction q_;
mutable ParGridFunction u_;
mutable ParGridFunction dqdt_;
mutable ParGridFunction dudt_;
mutable Vector X_;
mutable Vector RHS_;
mutable Vector rhs_;
mutable Vector dQs_;
mutable Vector tmp_;
Array<int> * bdr_attr_;
Array<int> ess_bdr_tdofs_;
VectorCoefficient * dqdtBdrCoef_;
bool tdQ_;
bool tdC_;
bool tdK_;
/*
bool ownsQ_;
bool ownsC_;
bool ownsK_;
*/
Coefficient * QCoef_;
Coefficient * CCoef_;
Coefficient * kCoef_;
MatrixCoefficient * KCoef_;
Coefficient * CInvCoef_;
Coefficient * kInvCoef_;
MatrixCoefficient * KInvCoef_;
Coefficient * dtCInvCoef_;
// MatrixCoefficient * dtKCoef_;
};
} // namespace thermal
class InverseCoefficient : public Coefficient
{
public:
InverseCoefficient(Coefficient & c) : c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return 1.0 / c_->Eval(T, ip); }
private:
Coefficient * c_;
};
class MatrixInverseCoefficient :public MatrixCoefficient
{
public:
MatrixInverseCoefficient(MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
MatrixCoefficient * M_;
};
class ScaledCoefficient : public Coefficient
{
public:
ScaledCoefficient(double a, Coefficient & c) : a_(a), c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return a_ * c_->Eval(T, ip); }
private:
double a_;
Coefficient * c_;
};
class ScaledMatrixCoefficient :public MatrixCoefficient
{
public:
ScaledMatrixCoefficient(double a, MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), a_(a), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
double a_;
MatrixCoefficient * M_;
};
} // namespace mfem
#endif // MFEM_USE_MPI
#endif // MFEM_FOURIER_FLUX_SOLVER
+818
View File
@@ -0,0 +1,818 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
//
// -----------------------------------------------------
// Fourier Miniapp: Thermal Diffusion
// -----------------------------------------------------
//
// This miniapp solves a time dependent heat equation.
//
#include "fourier_hybrid_solver.hpp"
#include <memory>
#include <iostream>
#include <fstream>
using namespace std;
using namespace mfem;
using namespace mfem::thermal;
void display_banner(ostream & os);
static int prob_ = 1;
static int unit_vec_type_ = 1;
static bool non_linear_ = false;
static double alpha_ = NAN;
static double theta_ = NAN;
static double gamma_ = 10.0;
static double chi_perp_ = 1.0;
static double chi_para_ = 1.0;
static double a_ = 0.15;
static double b_ = 0.85;
static double xc_ = 0.0;
static double yc_ = 0.0;
double TFunc(const Vector &x, double t)
{
switch (prob_)
{
case 1:
return x[0] * x[1] * pow(sin(M_PI * x[0]) * sin(M_PI * x[1]), gamma_);
case 2:
return 1.0 - pow(pow(x[0] - xc_, 2) + pow(x[1] - yc_, 2), 1.5);
case 3:
return 1.0 + (a_ * x[0] + b_ * x[1]) * pow(x[0] * x[0] + x[1] * x[1], 1.5);
case 4:
return 1.0 - pow(a_ * pow(x[0] * cos(theta_) + x[1] * sin(theta_), 2) +
b_ * pow(x[0] * sin(theta_) - x[1] * cos(theta_), 2), 1.5);
default:
return 0.0;
}
}
void qFunc(const Vector &x, Vector &q)
{
q.SetSize(2);
switch (prob_)
{
case 1:
{
double ssg = pow(sin(M_PI * x[0]) * sin(M_PI * x[1]), gamma_ - 1.0);
double ca = cos(alpha_);
double sa = sin(alpha_);
double cx = cos(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sx = sin(M_PI * x[0]);
double sy = sin(M_PI * x[1]);
double xcx = sx + M_PI * gamma_ * x[0] * cx;
double ycy = sy + M_PI * gamma_ * x[1] * cy;
double cd = chi_para_ - chi_perp_;
double cdca = cd * ca * ca + chi_perp_;
double cdsa = cd * sa * sa + chi_perp_;
q[0] = - x[0] * cd * ca * sa * sx * ycy - x[1] * cdca * sy * xcx;
q[1] = - x[1] * cd * ca * sa * sy * xcx - x[0] * cdsa * sx * ycy;
q *= ssg;
}
break;
default:
q = 0.0;
}
}
void qParaFunc(const Vector &x, Vector &q)
{
q.SetSize(2);
switch (prob_)
{
case 1:
{
double ssg = pow(sin(M_PI * x[0]) * sin(M_PI * x[1]), gamma_ - 1.0);
double ca = cos(alpha_);
double sa = sin(alpha_);
double cx = cos(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sx = sin(M_PI * x[0]);
double sy = sin(M_PI * x[1]);
double xcx = sx + M_PI * gamma_ * x[0] * cx;
double ycy = sy + M_PI * gamma_ * x[1] * cy;
double cd = chi_para_;
double cdca = cd * ca * ca;
double cdsa = cd * sa * sa;
q[0] = - x[0] * cd * ca * sa * sx * ycy - x[1] * cdca * sy * xcx;
q[1] = - x[1] * cd * ca * sa * sy * xcx - x[0] * cdsa * sx * ycy;
q *= ssg;
}
break;
default:
q = 0.0;
}
}
void qPerpFunc(const Vector &x, Vector &q)
{
q.SetSize(2);
switch (prob_)
{
case 1:
{
double ssg = pow(sin(M_PI * x[0]) * sin(M_PI * x[1]), gamma_ - 1.0);
double ca = cos(alpha_);
double sa = sin(alpha_);
double cx = cos(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sx = sin(M_PI * x[0]);
double sy = sin(M_PI * x[1]);
double xcx = sx + M_PI * gamma_ * x[0] * cx;
double ycy = sy + M_PI * gamma_ * x[1] * cy;
double cd = - chi_perp_;
double cdca = cd * ca * ca + chi_perp_;
double cdsa = cd * sa * sa + chi_perp_;
q[0] = - x[0] * cd * ca * sa * sx * ycy - x[1] * cdca * sy * xcx;
q[1] = - x[1] * cd * ca * sa * sy * xcx - x[0] * cdsa * sx * ycy;
q *= ssg;
}
break;
default:
q = 0.0;
}
}
void UnitBFunc(const Vector &x, Vector &b)
{
switch (unit_vec_type_)
{
case 2:
{
b[0] = -x[1] + yc_;
b[1] = x[0] - xc_;
}
break;
case 3:
{
b[0] = -3.0 * a_ * x[0] * x[1] -
b_ * (x[0] * x[0] + 4.0 * x[1] * x[1]);
b[1] = a_ * (4.0 * x[0] * x[0] + x[1] * x[1]) + 3.0 * b_ * x[0] * x[1];
}
break;
case 4:
{
double ct = cos(theta_);
double st = sin(theta_);
double ctst = 0.5 * sin(2.0 * theta_);
b[0] = x[1] * (a_ * st * st + b_ * ct * ct) + (a_ - b_) * x[0] * ctst;
b[1] = -x[0] * (a_ * ct * ct + b_ * st * st) - (a_ - b_) * x[1] * ctst;
}
break;
default:
b[0] = cos(alpha_);
b[1] = sin(alpha_);
}
double nrm = b.Norml2();
if ( nrm > 0.0 ) { b /= nrm; }
}
double QFunc(const Vector &x, double t)
{
switch (prob_)
{
case 1:
{
double cx = cos(M_PI * x[0]);
double sx = sin(M_PI * x[0]);
double s2x = sin(2.0 * M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sy = sin(M_PI * x[1]);
double s2y = sin(2.0 * M_PI * x[1]);
double ca = cos(alpha_);
double sa = sin(alpha_);
double s2a = sin(2.0 * alpha_);
double chi_sc = chi_perp_ * sa * sa + chi_para_ * ca * ca;
double chi_cs = chi_perp_ * ca * ca + chi_para_ * sa * sa;
double chi_s2 = (chi_para_ - chi_perp_) * s2a;
double s2gcx = s2x + M_PI * x[0] * (gamma_ * cx * cx - 1.0);
double s2gcy = s2y + M_PI * x[1] * (gamma_ * cy * cy - 1.0);
double sgcx = sx + M_PI * x[0] * gamma_ * cx;
double sgcy = sy + M_PI * x[1] * gamma_ * cy;
return -1.0 * (M_PI * gamma_ * x[0] * chi_cs * s2gcy * sx * sx +
M_PI * gamma_ * x[1] * chi_sc * s2gcx * sy * sy +
chi_s2 * sgcx * sgcy * sx * sy) *
pow(sx * sy, gamma_ - 2.0);
}
case 2:
{
return 9.0 * chi_perp_ * sqrt(pow(x[0] - xc_, 2) + pow(x[1] - yc_, 2));
}
default:
return 0.0;
}
}
void shiftUnitSquare(const Vector &x, Vector &p)
{
p[0] = x[0] - 0.5;
p[1] = x[1] - 0.5;
}
int main(int argc, char *argv[])
{
// 1. Initialize MPI.
MPI_Session mpi(argc, argv);
int myid = mpi.WorldRank();
// print the cool banner
if (mpi.Root()) { display_banner(cout); }
// 2. Parse command-line options.
int n = -1;
int order = 1;
int irOrder = -1;
int el_type = Element::QUADRILATERAL;
int ode_solver_type = 1;
int coef_type = 0;
int vis_steps = 1;
double dt = 0.5;
double t_final = 5.0;
double tol = 1e-4;
const char *basename = "FourierHybrid";
const char *mesh_file = "";
bool zero_start = true;
bool static_cond = false;
bool gfprint = true;
bool visit = true;
bool visualization = true;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&n, "-n", "--num-elems-1d",
"Number of elements in x and y directions. "
"Total number of elements is n^2.");
args.AddOption(&prob_, "-p", "--problem",
"Specify problem type: 1 - Square, 2 - Ellipse.");
args.AddOption(&coef_type, "-c", "--coef",
"Specify diffusion coefficient type: "
"0 - Constant, 1 - Linearized, 2 - Non-Linear.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree) or -1 for"
" isoparametric space.");
args.AddOption(&irOrder, "-iro", "--int-rule-order",
"Integration Rule Order.");
args.AddOption(&alpha_, "-alpha", "--constant-angle",
"Angle for constant B field (in degrees)");
args.AddOption(&theta_, "-theta", "--tilt-angle",
"Angle for orientation of ellipse (in degrees)");
args.AddOption(&a_, "-a", "--ellipse-a",
"First size parameter for ellipse");
args.AddOption(&b_, "-b", "--ellipse-b",
"Second size parameter for ellipse");
args.AddOption(&xc_, "-xc", "--x-center",
"x coordinate of field center");
args.AddOption(&yc_, "-yc", "--y-center",
"y coordinate of field center");
args.AddOption(&chi_perp_, "-chi-perp", "--chi-perpendicular",
"Chi perpendicular to field lines.");
args.AddOption(&chi_para_, "-chi-para", "--chi-parallel",
"Chi along field lines.");
args.AddOption(&dt, "-dt", "--time-step",
"Time step.");
args.AddOption(&t_final, "-tf", "--final-time",
"Final Time.");
args.AddOption(&tol, "-tol", "--tolerance",
"Tolerance used to determine convergence to steady state.");
args.AddOption(&el_type, "-e", "--element-type",
"Element type: 2-Triangle, 3-Quadrilateral.");
args.AddOption(&ode_solver_type, "-s", "--ode-solver",
"ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3\n\t."
"\t 22 - Mid-Point, 23 - SDIRK23, 34 - SDIRK34.");
args.AddOption(&zero_start, "-z", "--zero-start", "-no-z",
"--no-zero-start",
"Initial guess of zero or exact solution.");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&gfprint, "-print", "--print","-no-print","--no-print",
"Print results (grid functions) to disk.");
args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
"Enable or disable VisIt visualization.");
args.AddOption(&vis_steps, "-vs", "--visualization-steps",
"Visualize every n-th timestep.");
args.AddOption(&basename, "-k", "--outputfilename",
"Name of the visit dump files");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
if (irOrder < 0)
{
irOrder = std::max(4, 2 * order - 2);
}
if (isnan(alpha_))
{
alpha_ = 0.0;
}
else
{
alpha_ *= M_PI / 180.0;
}
unit_vec_type_ = prob_;
non_linear_ = coef_type > 0;
// 3. Construct a (serial) mesh of the given size on all processors. We
// can handle triangular and quadrilateral surface meshes with the
// same code.
Mesh *mesh = (n > 0) ?
new Mesh(n, n, (Element::Type)el_type, 1) :
new Mesh(mesh_file, 1, 1);
int dim = mesh->Dimension();
if (prob_ > 1) { mesh->Transform(shiftUnitSquare); }
// 4. This step is no longer needed
// 5. Define a parallel mesh by a partitioning of the serial mesh. Once the
// parallel mesh is defined, the serial mesh can be deleted.
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
delete mesh;
// 7. Determine the list of true (i.e. parallel conforming) essential
// boundary dofs. In this example, the boundary conditions are defined
// by marking all the boundary attributes from the mesh as essential
// (Dirichlet) and converting them to a list of true dofs.
Array<int> ess_tdof_list;
Array<int> ess_bdr(0);
if (pmesh->bdr_attributes.Size())
{
ess_bdr.SetSize(pmesh->bdr_attributes.Max());
ess_bdr = 1;
}
// The following is required for mesh refinement
// mesh->EnsureNCMesh();
// 6. Define the ODE solver used for time integration. Several implicit
// methods are available, including singly diagonal implicit Runge-Kutta
// (SDIRK).
ODESolver *ode_solver;
switch (ode_solver_type)
{
// Implicit L-stable methods
case 1: ode_solver = new BackwardEulerSolver; break;
case 2: ode_solver = new SDIRK23Solver(2); break;
case 3: ode_solver = new SDIRK33Solver; break;
// Implicit A-stable methods (not L-stable)
case 22: ode_solver = new ImplicitMidpointSolver; break;
case 23: ode_solver = new SDIRK23Solver; break;
case 34: ode_solver = new SDIRK34Solver; break;
default:
if (mpi.Root())
{
cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
}
delete mesh;
return 3;
}
// 12. Define the parallel finite element spaces. We use:
//
// H(curl) for electric field,
// H(div) for magnetic flux,
// H(div) for thermal flux,
// H(grad)/H1 for electrostatic potential,
// L2 for temperature
// L2 contains discontinuous "cell-center" finite elements, type 2 is
// "positive"
L2_FECollection L2FEC0(0, dim);
L2_FECollection L2FEC(order-1, dim);
// RT contains Raviart-Thomas "face-centered" vector finite elements with
// continuous normal component.
RT_FECollection HDivFEC(order-1, dim);
ND_FECollection HCurlFEC(order, dim);
// H1 contains continuous "node-centered" Lagrange finite elements.
H1_FECollection HGradFEC(order, dim);
ParFiniteElementSpace L2FESpace0(pmesh, &L2FEC0);
ParFiniteElementSpace L2FESpace(pmesh, &L2FEC);
ParFiniteElementSpace HDivFESpace(pmesh, &HDivFEC);
ParFiniteElementSpace HCurlFESpace(pmesh, &HCurlFEC);
ParFiniteElementSpace HGradFESpace(pmesh, &HGradFEC);
// The terminology is TrueVSize is the unique (non-redundant) number of dofs
// HYPRE_Int glob_size_l2 = L2FESpace.GlobalTrueVSize();
// HYPRE_Int glob_size_rt = HDivFESpace.GlobalTrueVSize();
HYPRE_Int glob_size_h1 = HGradFESpace.GlobalTrueVSize();
HYPRE_Int glob_size_rt = HDivFESpace.GlobalTrueVSize();
HYPRE_Int glob_size_l2 = L2FESpace.GlobalTrueVSize();
if (mpi.Root())
{
cout << "Number of Temperature unknowns: " << glob_size_h1 << endl;
cout << "Number of Heat Flux unknowns: " << glob_size_rt << endl;
cout << "Number of Thermal Energy unknowns: " << glob_size_l2 << endl;
}
// int Vsize_l2 = L2FESpace.GetVSize();
// int Vsize_rt = HDivFESpace.GetVSize();
// int Vsize_h1 = HGradFESpace.GetVSize();
// grid functions E, B, T, F, P, and w which is the Joule heating
ParGridFunction T_gf(&HGradFESpace);
ParGridFunction q_gf(&HDivFESpace);
ParGridFunction qPerpT_gf(&HDivFESpace);
ParGridFunction qParaT_gf(&HDivFESpace);
ParGridFunction qPerp_gf(&HDivFESpace);
ParGridFunction qPara_gf(&HDivFESpace);
ParGridFunction b_gf(&HDivFESpace);
ParGridFunction dT_gf(&HGradFESpace);
ParGridFunction Qs_gf(&HGradFESpace);
ParGridFunction errorT(&L2FESpace0);
ParGridFunction errorq(&L2FESpace0);
ParGridFunction errorqPerp(&L2FESpace0);
ParGridFunction errorqPara(&L2FESpace0);
ParGridFunction errorqPerpT(&L2FESpace0);
ParGridFunction errorqParaT(&L2FESpace0);
T_gf = 0.0;
q_gf = 0.0;
dT_gf = 1.0;
// 13. Get the boundary conditions, set up the exact solution grid functions
// These VectorCoefficients have an Eval function. Note that e_exact and
// b_exact in this case are exact analytical solutions, taking a 3-vector
// point as input and returning a 3-vector field
FunctionCoefficient TCoef(TFunc);
VectorFunctionCoefficient qCoef(2, qFunc);
VectorFunctionCoefficient qParaCoef(2, qParaFunc);
VectorFunctionCoefficient qPerpCoef(2, qPerpFunc);
Vector zeroVec(2); zeroVec = 0.0;
ConstantCoefficient zeroCoef(0.0);
VectorConstantCoefficient zeroVecCoef(zeroVec);
ConstantCoefficient SpecificHeatCoef(1.0);
// MatrixFunctionCoefficient ConductionCoef(2, ChiFunc);
FunctionCoefficient HeatSourceCoef(QFunc);
VectorFunctionCoefficient UnitBCoef(2, UnitBFunc);
b_gf.ProjectCoefficient(UnitBCoef);
Qs_gf.ProjectCoefficient(HeatSourceCoef);
T_gf.ProjectCoefficient(TCoef);
q_gf.ProjectCoefficient(qCoef);
qPara_gf.ProjectCoefficient(qParaCoef);
qPerp_gf.ProjectCoefficient(qPerpCoef);
double T_nrm = T_gf.ComputeL2Error(zeroCoef);
double q_nrm = q_gf.ComputeL2Error(zeroVecCoef);
double qPara_nrm = qPara_gf.ComputeL2Error(zeroVecCoef);
double qPerp_nrm = qPerp_gf.ComputeL2Error(zeroVecCoef);
T_gf.ProjectBdrCoefficient(TCoef, ess_bdr);
q_gf.ProjectBdrCoefficientNormal(qCoef, ess_bdr);
T_gf.GridFunction::ComputeElementL2Errors(TCoef, errorT);
q_gf.GridFunction::ComputeElementL2Errors(qCoef, errorq);
// 14. Initialize the Diffusion operator, the GLVis visualization and print
// the initial energies.
cout << "Building TDO" << endl;
HybridThermalDiffusionTDO oper(HGradFESpace,
HCurlFESpace,
HDivFESpace,
L2FESpace,
zeroVecCoef,
zeroCoef, ess_bdr,
chi_perp_,
chi_para_,
prob_,
coef_type,
UnitBCoef,
SpecificHeatCoef, false,
// ConductionCoef, false,
HeatSourceCoef, false);
// This function initializes all the fields to zero or some provided IC
// oper.Init(F);
socketstream vis_T, vis_q, vis_b, vis_Q, vis_errT, vis_errq;
char vishost[] = "localhost";
int visport = 19916;
if (visualization)
{
// Make sure all ranks have sent their 'v' solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
vis_T.precision(8);
vis_Q.precision(8);
vis_q.precision(8);
vis_b.precision(8);
vis_errT.precision(8);
vis_errq.precision(8);
int Wx = 0, Wy = 0; // window position
int Ww = 350, Wh = 350; // window size
int offx = Ww+10, offy = Wh+45; // window offsets
miniapps::VisualizeField(vis_Q, vishost, visport,
Qs_gf, "Heat Source", Wx, Wy, Ww, Wh);
Wy += offy;
miniapps::VisualizeField(vis_b, vishost, visport,
b_gf, "Unit B Field", Wx, Wy, Ww, Wh, true);
Wx += offx; Wy -= offy;
miniapps::VisualizeField(vis_T, vishost, visport,
T_gf, "Temperature", Wx, Wy, Ww, Wh);
Wy += offy;
miniapps::VisualizeField(vis_errT, vishost, visport,
errorT, "Error in T", Wx, Wy, Ww, Wh);
Wx += offx; Wy -= offy;
miniapps::VisualizeField(vis_q, vishost, visport,
q_gf, "Heat Flux", Wx, Wy, Ww, Wh, true);
Wy += offy;
miniapps::VisualizeField(vis_errq, vishost, visport,
errorq, "Error in q", Wx, Wy, Ww, Wh);
}
// VisIt visualization
VisItDataCollection visit_dc(basename, pmesh);
if ( visit )
{
visit_dc.RegisterField("T", &T_gf);
visit_dc.RegisterField("Qs", &Qs_gf);
visit_dc.RegisterField("q", &q_gf);
visit_dc.RegisterField("qPerp", &qPerp_gf);
visit_dc.RegisterField("qPara", &qPara_gf);
visit_dc.RegisterField("qPerpT", &qPerpT_gf);
visit_dc.RegisterField("qParaT", &qParaT_gf);
visit_dc.RegisterField("b", &b_gf);
visit_dc.RegisterField("L2 Error T", &errorT);
visit_dc.RegisterField("L2 Error q", &errorq);
visit_dc.RegisterField("L2 Error qPerp", &errorqPerp);
visit_dc.RegisterField("L2 Error qPara", &errorqPara);
visit_dc.RegisterField("L2 Error qPerpT", &errorqPerpT);
visit_dc.RegisterField("L2 Error qParaT", &errorqParaT);
oper.SetVisItDC(visit_dc);
visit_dc.SetCycle(0);
visit_dc.SetTime(0.0);
visit_dc.Save();
}
ostringstream oss_errs;
oss_errs << "fourier_hybrid_errs"
<< "_p" << prob_ << "_c" << coef_type
<< "_e" << (int)floor(log10(chi_para_/chi_perp_));
if (n > 0) { oss_errs << "_n" << n; }
oss_errs << "_o" << order << ".dat";
ofstream ofs_errs;
if (myid == 0) { ofs_errs.open(oss_errs.str().c_str()); }
// 15. Perform time-integration (looping over the time iterations, ti, with a
// time-step dt). The object oper is the MagneticDiffusionOperator which
// has a Mult() method and an ImplicitSolve() method which are used by
// the time integrators.
ode_solver->Init(oper);
double t = 0.0;
int tsize = HGradFESpace.GetTrueVSize();
int qsize = HDivFESpace.GetTrueVSize();
Vector X0(tsize+qsize), X1(tsize+qsize), dX(tsize+qsize);
Vector T1(X1.GetData(), tsize);
Vector q1(&(X1.GetData())[tsize], qsize);
X0 = 0.0; X1 = 0.0; dX = 0.0;
T_gf.ParallelProject(T1);
bool last_step = false;
for (int ti = 1; !last_step; ti++)
{
if (t + dt >= t_final - dt/2)
{
if (myid == 0)
{
cout << "Final Time Reached" << endl;
}
last_step = true;
}
// F is the vector of dofs, t is the current time, and dt is the time step
// to advance.
X0 = X1;
ode_solver->Step(X1, t, dt);
T_gf.Distribute(T1);
q_gf.Distribute(q1);
TCoef.SetTime(t);
oper.GetParaFluxFromTemp(T_gf, qParaT_gf);
oper.GetPerpFluxFromTemp(T_gf, qPerpT_gf);
oper.GetParaFluxFromFlux(q_gf, qPara_gf);
oper.GetPerpFluxFromFlux(q_gf, qPerp_gf);
T_gf.GridFunction::ComputeElementL2Errors(TCoef, errorT);
q_gf.GridFunction::ComputeElementL2Errors(qCoef, errorq);
qPerp_gf.GridFunction::ComputeElementL2Errors(qPerpCoef, errorqPerp);
qPara_gf.GridFunction::ComputeElementL2Errors(qParaCoef, errorqPara);
qPerpT_gf.GridFunction::ComputeElementL2Errors(qPerpCoef, errorqPerpT);
qParaT_gf.GridFunction::ComputeElementL2Errors(qParaCoef, errorqParaT);
double l2_error_T = T_gf.ComputeL2Error(TCoef);
double l2_error_q = q_gf.ComputeL2Error(qCoef);
if ( myid == 0 )
{
ofs_errs << t << '\t' << l2_error_T << '\t' << l2_error_q << endl;
cout << t << '\t' << l2_error_T << '\t' << l2_error_q << endl;
}
add(1.0, X1, -1.0, X0, dX);
Vector dT(dX.GetData(), tsize);
dT_gf.Distribute(dT);
double maxT = T_gf.ComputeMaxError(zeroCoef);
double maxDiff = dT_gf.ComputeMaxError(zeroCoef);
if ( !last_step )
{
if ( maxT == 0.0 )
{
last_step = (maxDiff < tol) ? true:false;
}
else if ( maxDiff/maxT < tol )
{
last_step = true;
}
if (last_step && myid == 0)
{
cout << "Converged to Steady State" << endl;
}
}
/*
if (debug == 1)
{
oper.Debug(basename,t);
}
*/
if (gfprint)
{
ostringstream T_name, q_name, mesh_name;
T_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "T." << setfill('0') << setw(6) << myid;
q_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "q." << setfill('0') << setw(6) << myid;
mesh_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "mesh." << setfill('0') << setw(6) << myid;
ofstream mesh_ofs(mesh_name.str().c_str());
mesh_ofs.precision(8);
pmesh->Print(mesh_ofs);
mesh_ofs.close();
ofstream T_ofs(T_name.str().c_str());
T_ofs.precision(8);
T_gf.Save(T_ofs);
T_ofs.close();
ofstream q_ofs(q_name.str().c_str());
q_ofs.precision(8);
q_gf.Save(q_ofs);
q_ofs.close();
}
if (last_step || (ti % vis_steps) == 0)
{
// Make sure all ranks have sent their 'v' solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
if (visualization)
{
int Wx = 0, Wy = 0; // window position
int Ww = 350, Wh = 350; // window size
miniapps::VisualizeField(vis_T, vishost, visport,
T_gf, "Temperature", Wx, Wy, Ww, Wh);
miniapps::VisualizeField(vis_q, vishost, visport,
q_gf, "Heat Flux", Wx, Wy, Ww, Wh);
miniapps::VisualizeField(vis_errT, vishost, visport,
errorT, "Error in T", Wx, Wy, Ww, Wh);
miniapps::VisualizeField(vis_errq, vishost, visport,
errorq, "Error in q", Wx, Wy, Ww, Wh);
}
if (visit)
{
visit_dc.SetCycle(ti);
visit_dc.SetTime(t);
visit_dc.Save();
}
}
}
// oper.GetParaFluxFromTemp(T_gf, qParaT_gf);
// oper.GetPerpFluxFromTemp(T_gf, qPerpT_gf);
if (visualization)
{
vis_T.close();
vis_q.close();
vis_errT.close();
vis_errq.close();
}
if (myid == 0) { ofs_errs.close(); }
/*
double loc_T_max = T1.Normlinf();
double T_max = -1.0;
MPI_Allreduce(&loc_T_max, &T_max, 1, MPI_DOUBLE, MPI_MAX,
MPI_COMM_WORLD);
*/
double err1 = T_gf.ComputeL2Error(TCoef);
double errq = q_gf.ComputeL2Error(qCoef);
double errqParaT = qParaT_gf.ComputeL2Error(qParaCoef);
double errqPerpT = qPerpT_gf.ComputeL2Error(qPerpCoef);
double errqPara = qPara_gf.ComputeL2Error(qParaCoef);
double errqPerp = qPerp_gf.ComputeL2Error(qPerpCoef);
double T_max = T_gf.ComputeMaxError(zeroCoef);
double q_max = q_gf.ComputeMaxError(zeroVecCoef);
// double qParaT_max = qParaT_gf.ComputeMaxError(zeroVecCoef);
// double qPerpT_max = qPerpT_gf.ComputeMaxError(zeroVecCoef);
// double qPara_max = qPara_gf.ComputeMaxError(zeroVecCoef);
// double qPerp_max = qPerp_gf.ComputeMaxError(zeroVecCoef);
if (myid == 0)
{
cout << "Maximum Temperature: " << T_max << endl;
cout << "Maximum Flux Magnitude: " << q_max << endl;
cout << "L2 Error of Temperature: " << err1
<< ", (relative " << err1 / T_nrm << ")"
<< endl;
cout << "L2 Error of Flux: " << errq
<< ", (relative " << errq / q_nrm << ")"
<< endl;
cout << "L2 Error of Para Flux: " << errqPara
<< ", (relative " << errqPara / qPara_nrm << ")"
<< endl;
cout << "L2 Error of Perp Flux: " << errqPerp
<< ", (relative " << errqPerp / qPerp_nrm << ")"
<< endl;
cout << "L2 Error of Para Flux T: " << errqParaT
<< ", (relative " << errqParaT / qPara_nrm << ")"
<< endl;
cout << "L2 Error of Perp Flux T: " << errqPerpT
<< ", (relative " << errqPerpT / qPerp_nrm << ")"
<< endl;
cout << "| chi_eff - 1 | = " << fabs(1.0/T_max - 1) << endl;
}
// 16. Free the used memory.
delete ode_solver;
delete pmesh;
return 0;
}
void display_banner(ostream & os)
{
os << "___________ .__ " << endl
<< "\\_ _____/___ __ _________|__| ___________ " << endl
<< " | __)/ _ \\| | \\_ __ \\ |/ __ \\_ __ \\" << endl
<< " | | ( <_> ) | /| | \\/ \\ ___/| | \\/" << endl
<< " \\__ | \\____/|____/ |__| |__|\\___ >__| " << endl
<< " \\/ \\/ " << endl
<< flush;
}
+914
View File
@@ -0,0 +1,914 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#include "fourier_hybrid_solver.hpp"
#ifdef MFEM_USE_MPI
using namespace std;
namespace mfem
{
using namespace miniapps;
void ChiPerpCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
K *= -1.0;
K(0,0) += 1.0;
K(1,1) += 1.0;
if (nonlin_)
{
K *= 1.0 / sqrt(fabs(T_->Eval(T, ip)));
}
K *= chi_perp_;
}
void ChiParaCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
if (nonlin_)
{
K *= pow(fabs(T_->Eval(T, ip)), 2.5);
}
K *= chi_para_;
}
void dChiParaCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
double temp = T_->Eval(T, ip);
double para_factor = 2.5 * chi_para_ * pow(fabs(temp), 1.5);
bbT_->Eval(K, T, ip);
K *= para_factor;
}
void dChiCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
double temp = fabs(T_->Eval(T, ip));
double perp_factor = 0.5 * chi_perp_ * pow(temp, -1.5);
double para_factor = 2.5 * chi_para_ * pow(temp, 1.5);
bbT_->Eval(K, T, ip);
K *= perp_factor + para_factor;
K(0,0) -= perp_factor;
K(1,1) -= perp_factor;
}
void ChiInvPerpCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
K *= -1.0;
K(0,0) += 1.0;
K(1,1) += 1.0;
if (nonlin_)
{
K *= sqrt(fabs(T_->Eval(T, ip)));
}
K *= 1.0 / chi_perp_;
}
void ChiInvParaCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
if (nonlin_)
{
K *= pow(fabs(T_->Eval(T, ip)), -2.5);
}
K *= 1.0 / chi_para_;
}
namespace thermal
{
HybridThermalDiffusionTDO::HybridThermalDiffusionTDO(
ParFiniteElementSpace &H1_FESpace,
ParFiniteElementSpace &HCurl_FESpace,
ParFiniteElementSpace &HDiv_FESpace,
ParFiniteElementSpace &L2_FESpace,
VectorCoefficient & dqdtBdr,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
double chi_perp,
double chi_para,
int prob,
int coef_type,
VectorCoefficient & UnitB,
Coefficient & c, bool td_c,
Coefficient & Q, bool td_Q)
: TimeDependentOperator(H1_FESpace.GetTrueVSize() +
HDiv_FESpace.GetTrueVSize(), 0.0),
init_(false),
nonLinear_(coef_type == 2),
testGradient_(false),
dim_(H1_FESpace.GetParMesh()->Dimension()),
tsize_(H1_FESpace.GetTrueVSize()),
qsize_(HDiv_FESpace.GetTrueVSize()),
multCount_(0), solveCount_(0),
T_(&H1_FESpace),
dT_(&H1_FESpace),
q_(&HDiv_FESpace),
Q_perp_(&L2_FESpace),
TCoef_(&T_),
unitBCoef_(&UnitB),
bbTCoef_(*unitBCoef_, *unitBCoef_),
ICoef_(dim_),
PPerpCoef_(bbTCoef_, ICoef_, -1.0),
chiPerpCoef_(bbTCoef_, TCoef_, chi_perp, coef_type != 0),
chiParaCoef_(bbTCoef_, TCoef_, chi_para, coef_type != 0),
chiCoef_(chiPerpCoef_, chiParaCoef_),
dChiCoef_(bbTCoef_, TCoef_, chi_perp, chi_para),
dChiParaCoef_(bbTCoef_, TCoef_, chi_para),
chiInvPerpCoef_(bbTCoef_, TCoef_, chi_perp, coef_type != 0),
chiInvParaCoef_(bbTCoef_, TCoef_, chi_para, coef_type != 0),
chiInvCoef_(chiInvPerpCoef_, chiInvParaCoef_),
H1_FESpace_(&H1_FESpace),
HCurl_FESpace_(&HCurl_FESpace),
HDiv_FESpace_(&HDiv_FESpace),
L2_FESpace_(&L2_FESpace),
m2_(NULL), mPara_(NULL), mPerp_(NULL), sC_(NULL), dC_(NULL), a_(NULL),
gPerp_(NULL), gPara_(NULL),
Div_(NULL),
Grad_(NULL),
dqdt_gf_(NULL), Qs_(NULL),
M2Inv_(NULL), M2Diag_(NULL),
AInv_(NULL), APrecond_(NULL),
dqdt_(&HDiv_FESpace),
dqdt_perp_(&HDiv_FESpace),
dqdt_para_(&HDiv_FESpace),
dqdt_from_T_(&HCurl_FESpace),
dqdt_para_from_T_(&HDiv_FESpace),
q1_perp_(&HDiv_FESpace),
dqdt_perp_dual_(&HDiv_FESpace),
dqdt_para_dual_(&HDiv_FESpace),
// rhs_(NULL),
bdr_attr_(&bdr_attr), ess_bdr_tdofs_(0), dqdtBdrCoef_(&dqdtBdr),
tdQ_(td_Q), tdC_(td_c),
QCoef_(&Q), CCoef_(&c),
CInvCoef_(new InverseCoefficient(c)),
dtCInvCoef_(NULL),
impOp_(H1_FESpace,
dTdtBdr, false,
bdr_attr,
c, td_c,
chiCoef_, coef_type > 0,
dChiCoef_, coef_type > 0,
Q, td_Q || true,
coef_type == 2),
newton_(H1_FESpace.GetComm())
{
this->init();
}
HybridThermalDiffusionTDO::~HybridThermalDiffusionTDO()
{
delete CInvCoef_;
delete dtCInvCoef_;
delete Div_;
delete Grad_;
delete dC_;
delete a_;
delete gPara_;
delete gPerp_;
delete m2_;
delete mPara_;
delete mPerp_;
delete sC_;
delete dqdt_gf_;
delete Qs_;
delete M2Inv_;
delete M2Diag_;
delete AInv_;
delete APrecond_;
}
void
HybridThermalDiffusionTDO::SetVisItDC(VisItDataCollection & visit_dc)
{
visit_dc.RegisterField("Q_perp", &Q_perp_);
visit_dc.RegisterField("dqdt_para", &dqdt_para_);
visit_dc.RegisterField("dqdt_perp", &dqdt_perp_);
visit_dc.RegisterField("dqdt T", &dqdt_from_T_);
visit_dc.RegisterField("dqdt_para T", &dqdt_para_from_T_);
}
void
HybridThermalDiffusionTDO::init()
{
cout << "Entering TDO::Init" << endl;
if ( init_ ) { return; }
if ( m2_ == NULL )
{
m2_ = new ParBilinearForm(HDiv_FESpace_);
m2_->AddDomainIntegrator(new VectorFEMassIntegrator());
m2_->Assemble();
}
if ( mPerp_ == NULL )
{
mPerp_ = new ParBilinearForm(HDiv_FESpace_);
mPerp_->AddDomainIntegrator(new VectorFEMassIntegrator(PPerpCoef_));
mPerp_->Assemble();
}
if ( mPara_ == NULL )
{
mPara_ = new ParBilinearForm(HDiv_FESpace_);
mPara_->AddDomainIntegrator(new VectorFEMassIntegrator(bbTCoef_));
mPara_->Assemble();
}
if ( sC_ == NULL )
{
sC_ = new ParBilinearForm(HDiv_FESpace_);
sC_->AddDomainIntegrator(new DivDivIntegrator(*CInvCoef_));
sC_->Assemble();
}
if ( dC_ == NULL )
{
dC_ = new ParMixedBilinearForm(L2_FESpace_, HDiv_FESpace_);
dC_->AddDomainIntegrator(
new MixedScalarWeakGradientIntegrator(*CInvCoef_));
dC_->Assemble();
}
if ( gPara_ == NULL )
{
gPara_ = new ParMixedBilinearForm(H1_FESpace_, HDiv_FESpace_);
gPara_->AddDomainIntegrator(
new MixedVectorGradientIntegrator(chiParaCoef_));
gPara_->Assemble();
}
if ( gPerp_ == NULL )
{
gPerp_ = new ParMixedBilinearForm(H1_FESpace_, HDiv_FESpace_);
gPerp_->AddDomainIntegrator(
new MixedVectorGradientIntegrator(chiPerpCoef_));
gPerp_->Assemble();
}
if ( dqdt_gf_ == NULL )
{
dqdt_gf_ = new ParGridFunction(HDiv_FESpace_);
}
if ( Qs_ == NULL && QCoef_ != NULL )
{
Qs_ = new ParGridFunction(L2_FESpace_);
Qs_->ProjectCoefficient(*QCoef_);
}
Div_ = new ParDiscreteDivOperator(HDiv_FESpace_, L2_FESpace_);
Div_->Assemble();
Div_->Finalize();
Grad_ = new ParDiscreteGradOperator(H1_FESpace_, HCurl_FESpace_);
Grad_->Assemble();
Grad_->Finalize();
rhs_.SetSize(HDiv_FESpace_->GetVSize());
dQs_.SetSize(HDiv_FESpace_->GetVSize());
// tmp_.SetSize(L2_FESpace_->GetVSize());
HDiv_FESpace_->GetEssentialTrueDofs(*bdr_attr_, ess_bdr_tdofs_);
newton_.SetPrintLevel(2);
newton_.SetRelTol(1e-10);
newton_.SetAbsTol(0.0);
if ( nonLinear_ && testGradient_ )
{
Vector x(impOp_.Height());
Vector dx(impOp_.Height());
T_.Distribute(x);
Q_perp_ = 0.0;
cout << "GetTime " << this->GetTime() << endl;
impOp_.SetState(T_, Q_perp_, this->GetTime(), 0.1);
cout << "init 0" << endl;
newton_.SetOperator(impOp_);
cout << "init 1" << endl;
cout << "init 2" << endl;
x.Randomize(1);
x.Print(cout);
dx.Randomize(2);
dx *= 0.01;
dx.Print(cout);
cout << "init 3" << endl;
double ratio = newton_.CheckGradient(x, dx);
cout << "CheckGradient returns: " << ratio << endl;
}
init_ = true;
cout << "Leaving TDO::Init" << endl;
}
void
HybridThermalDiffusionTDO::SetTime(const double time)
{
this->TimeDependentOperator::SetTime(time);
dqdtBdrCoef_->SetTime(t);
if ( tdQ_ )
{
QCoef_->SetTime(t);
Qs_->ProjectCoefficient(*QCoef_);
}
if ( tdC_ )
{
// CCoef_->SetTime(t);
// CInvCoef_->SetTime(t);
dtCInvCoef_->SetTime(t);
sC_->Assemble();
}
chiInvCoef_.SetTime(t);
if ( tdC_ && a_ != NULL )
{
a_->Assemble();
}
newTime_ = true;
}
void
HybridThermalDiffusionTDO::Mult(const Vector &T, Vector &dT_dt) const
{
MFEM_ABORT("HybridThermalDiffusionTDO::Mult should not be called");
}
void
HybridThermalDiffusionTDO::initA(double dt)
{
cout << "Entering initA" << endl;
if ( CInvCoef_ != NULL )
{
dtCInvCoef_ = new ScaledCoefficient(dt, *CInvCoef_);
}
if ( a_ == NULL)
{
a_ = new ParBilinearForm(HDiv_FESpace_);
a_->AddDomainIntegrator(new VectorFEMassIntegrator(chiInvCoef_));
a_->AddDomainIntegrator(new DivDivIntegrator(*dtCInvCoef_));
a_->Assemble();
}
else
{
a_->Update();
a_->Assemble();
}
cout << "Leaving initA" << endl;
}
void
HybridThermalDiffusionTDO::initImplicitSolve()
{
cout << "Entering initImplicitSolve" << endl;
// if ( tdC_ || AInv_ == NULL || APrecond_ == NULL )
{
delete AInv_;
AInv_ = new HyprePCG(A_);
AInv_->SetTol(1e-12);
AInv_->SetMaxIter(200);
AInv_->SetPrintLevel(0);
delete APrecond_;
APrecond_ = (dim_==2) ?
(HypreSolver*)(new HypreAMS(A_, HDiv_FESpace_)):
(HypreSolver*)(new HypreADS(A_, HDiv_FESpace_));
if ( dim_ == 2 )
{
dynamic_cast<HypreAMS*>(APrecond_)->SetPrintLevel(0);
}
else
{
dynamic_cast<HypreADS*>(APrecond_)->SetPrintLevel(0);
}
AInv_->SetPreconditioner(*APrecond_);
}
/*
else
{
AInv_->SetOperator(A_);
}
*/
if ( M2Inv_ == NULL )
{
Array<int> ess_tdof(0);
m2_->FormSystemMatrix(ess_tdof, M2_);
M2Inv_ = new HyprePCG(M2_);
M2Inv_->SetTol(1e-12);
M2Inv_->SetMaxIter(200);
M2Inv_->SetPrintLevel(0);
M2Diag_ = new HypreDiagScale(M2_);
M2Inv_->SetPreconditioner(*M2Diag_);
}
cout << "Leaving initImplicitSolve" << endl;
}
void
HybridThermalDiffusionTDO::ImplicitSolve(const double dt,
const Vector &X, Vector &dX_dt)
{
cout << "Entering ImplicitSolve" << endl;
Vector T(X.GetData(), tsize_);
Vector q(&(X.GetData())[tsize_], qsize_);
Vector dT_dt(dX_dt.GetData(), tsize_);
Vector dq_dt(&(dX_dt.GetData())[tsize_], qsize_);
cout << 1 << endl;
cout << "Norms of T and q: " << T.Norml2() << " " << q.Norml2() << endl;
dX_dt = 0.0;
cout << 2 << endl;
T_.Distribute(T);
{
// q_.MakeRef(const_cast<ParFiniteElementSpace*>(HDiv_FESpace_),
// const_cast<Vector&>(y), 0);
// u_.MakeRef(const_cast<ParFiniteElementSpace*>(L2_FESpace_),
// const_cast<Vector&>(y), HDiv_FESpace_->GetVSize());
q_.Distribute(q);
// dqdt_.MakeRef(HDiv_FESpace_, dy_dt, 0);
// dudt_.MakeRef(L2_FESpace_, dy_dt, HDiv_FESpace_->GetVSize());
// cout << "sC size: " << sC_->Width() << ", q_ size: " << q_.Size() << ", rhs_ size: " << rhs_.Size() << endl;
cout << 3 << endl;
sC_->Mult(q_, rhs_);
dC_->Mult(*Qs_, dQs_);
rhs_ += dQs_;
rhs_ *= -1.0;
cout << 4 << endl;
// dqdt_gf_->ProjectBdrCoefficientNormal(*dqdtBdrCoef_, *bdr_attr_);
dqdt_.ProjectBdrCoefficientNormal(*dqdtBdrCoef_, *bdr_attr_);
cout << 5 << endl;
this->initA(dt);
// a_->FormLinearSystem(ess_bdr_tdofs_, *dqdt_gf_, rhs_, A_, X_, RHS_);
a_->FormLinearSystem(ess_bdr_tdofs_, dqdt_, rhs_, A_, X_, RHS_);
this->initImplicitSolve();
AInv_->Mult(RHS_, X_);
a_->RecoverFEMSolution(X_, rhs_, dqdt_);
cout << "Norm of dqdt_: " << dqdt_.Normlinf() << endl;
dq_dt = X_;
Q_perp_ = 0.0;
/*
mPerp_->Mult(dqdt_, dqdt_perp_dual_);
cout << "Norm of dqdt_perp_dual_: " << dqdt_perp_dual_.Normlinf() << endl;
Vector RHS(qsize_);
Vector X(qsize_);
dqdt_perp_dual_.ParallelAssemble(RHS);
M2Inv_->Mult(RHS, dq_dt);
dqdt_perp_.Distribute(dq_dt);
dqdt_para_ = dqdt_;
dqdt_para_ -= dqdt_perp_;
mPerp_->Mult(q_, dqdt_perp_dual_);
dqdt_perp_dual_.ParallelAssemble(RHS);
M2Inv_->Mult(RHS, X);
q1_perp_.Distribute(X);
q1_perp_.Add(dt, dqdt_perp_);
// dq_dt = X;
cout << "Norm of dqdt_perp_: " << dqdt_perp_.Normlinf() << endl;
Div_->Mult(q1_perp_, Q_perp_);
// Q_perp_ += tmp_;
Q_perp_ *= 0.0;
// dudt_ += *Qs_;
*/
}
impOp_.SetState(T_, Q_perp_, this->GetTime(), dt);
Solver & solver = impOp_.GetGradientSolver();
if (!nonLinear_)
{
solver.Mult(impOp_.GetRHS(), dT_dt);
}
else
{
newton_.SetOperator(impOp_);
newton_.SetSolver(solver);
newton_.Mult(impOp_.GetRHS(), dT_dt);
}
if (false)
{
cout << 6 << endl;
dT_.Distribute(dT_dt);
T_.Add(dt, dT_);
cout << 7 << endl;
gPara_->Update();
gPara_->Assemble();
cout << 8 << endl;
gPara_->Mult(dT_, dqdt_para_dual_);
cout << 9 << endl;
Vector X(qsize_), RHS(qsize_);
dqdt_para_dual_.ParallelAssemble(RHS);
M2Inv_->Mult(RHS, X);
dqdt_para_from_T_.Distribute(X);
Grad_->Mult(dT_, dqdt_from_T_);
cout << "Norm of dqdt_para: " << X.Norml2() << endl;
cout << 10 << endl;
dq_dt += X;
}
cout << "Norms of dT and dq: " << dT_dt.Norml2() << " " << dq_dt.Norml2() <<
endl;
solveCount_++;
}
void
HybridThermalDiffusionTDO::GetParaFluxFromTemp(const ParGridFunction &T,
ParGridFunction & q_para)
{
gPara_->Mult(T, dqdt_para_dual_);
Vector X(qsize_), RHS(qsize_);
dqdt_para_dual_.ParallelAssemble(RHS);
RHS *= -1.0;
M2Inv_->Mult(RHS, X);
q_para.Distribute(X);
}
void
HybridThermalDiffusionTDO::GetPerpFluxFromTemp(const ParGridFunction &T,
ParGridFunction & q_perp)
{
gPerp_->Mult(T, dqdt_para_dual_);
Vector X(qsize_), RHS(qsize_);
dqdt_para_dual_.ParallelAssemble(RHS);
RHS *= -1.0;
M2Inv_->Mult(RHS, X);
q_perp.Distribute(X);
}
void
HybridThermalDiffusionTDO::GetParaFluxFromFlux(const ParGridFunction &q,
ParGridFunction & q_para)
{
mPara_->Mult(q, dqdt_perp_dual_);
Vector RHS(qsize_);
Vector X(qsize_);
dqdt_perp_dual_.ParallelAssemble(RHS);
M2Inv_->Mult(RHS, X);
q_para.Distribute(X);
}
void
HybridThermalDiffusionTDO::GetPerpFluxFromFlux(const ParGridFunction &q,
ParGridFunction & q_perp)
{
mPerp_->Mult(q, dqdt_perp_dual_);
Vector RHS(qsize_);
Vector X(qsize_);
dqdt_perp_dual_.ParallelAssemble(RHS);
M2Inv_->Mult(RHS, X);
q_perp.Distribute(X);
}
ImplicitDiffOp::ImplicitDiffOp(ParFiniteElementSpace & H1_FESpace,
Coefficient & dTdtBdr, bool tdBdr,
Array<int> & bdr_attr,
Coefficient & heatCap, bool tdCp,
MatrixCoefficient & chi, bool tdChi,
MatrixCoefficient & dchi, bool tdDChi,
Coefficient & heatSource, bool tdQ,
bool nonlinear)
: Operator(H1_FESpace.GetTrueVSize()),
first_(true),
tdBdr_(tdBdr),
tdCp_(tdCp),
tdChi_(tdChi),
tdDChi_(tdDChi),
tdQ_(tdQ),
nonLinear_(nonlinear),
newTime_(true),
newTimeStep_(true),
t_(0.0),
dt_(-1.0),
ess_bdr_attr_(bdr_attr),
bdrCoef_(&dTdtBdr),
cpCoef_(&heatCap),
chiCoef_(&chi),
dChiCoef_(&dchi),
chiNLCoef_(&dynamic_cast<NLCoefficient&>(chi)),
dChiNLCoef_(&dynamic_cast<NLCoefficient&>(dchi)),
QPerpCoef_(NULL),
QCoef_(heatSource, QPerpCoef_),
dtChiCoef_(1.0, *chiCoef_),
T0_(&H1_FESpace),
T1_(&H1_FESpace),
dT_(&H1_FESpace),
gradTCoef_(&T0_),
dtGradTCoef_(-1.0, gradTCoef_),
dtdChiGradTCoef_(*dChiCoef_, dtGradTCoef_),
m0cp_(&H1_FESpace),
s0chi_(&H1_FESpace),
a0_(&H1_FESpace),
dTdt_(&H1_FESpace),
Q_(&H1_FESpace),
Qs_(&H1_FESpace),
rhs_(&H1_FESpace),
RHS_(H1_FESpace.GetTrueVSize()),
// RHS0_(0),
AInv_(NULL),
APrecond_(NULL)
{
cout << "Entering ImplicitDiffOp c'tor" << endl;
H1_FESpace.GetEssentialTrueDofs(ess_bdr_attr_, ess_bdr_tdofs_);
m0cp_.AddDomainIntegrator(new MassIntegrator(*cpCoef_));
s0chi_.AddDomainIntegrator(new DiffusionIntegrator(*chiCoef_));
a0_.AddDomainIntegrator(new MassIntegrator(*cpCoef_));
a0_.AddDomainIntegrator(new DiffusionIntegrator(dtChiCoef_));
if (nonLinear_)
{
a0_.AddDomainIntegrator(new MixedScalarWeakDivergenceIntegrator(
dtdChiGradTCoef_));
}
cout << "Qs 0" << endl;
Qs_.AddDomainIntegrator(new DomainLFIntegrator(QCoef_));
cout << "Qs 1 " << tdQ_ << endl;
if (!tdQ_) { Qs_.Assemble(); }
cout << "Leaving ImplicitDiffOp c'tor" << endl;
}
ImplicitDiffOp::~ImplicitDiffOp()
{
delete AInv_;
delete APrecond_;
}
void ImplicitDiffOp::SetState(ParGridFunction & T, ParGridFunction & Q_perp,
double t, double dt)
{
T0_ = T;
newTime_ = fabs(t - t_) > 0.0;
newTimeStep_= (fabs(1.0-dt/dt_)>1e-6);
t_ = newTime_ ? t : t_;
dt_ = newTimeStep_ ? dt : dt_;
if (tdBdr_ && (newTime_ || newTimeStep_))
{
bdrCoef_->SetTime(t_ + dt_);
}
if (newTimeStep_ || first_)
{
dtChiCoef_.SetAConst(dt_);
dtGradTCoef_.SetAConst(-dt_);
}
if ((tdCp_ && newTime_) || first_)
{
m0cp_.Update();
m0cp_.Assemble();
m0cp_.Finalize();
}
if (!tdChi_ && first_)
{
s0chi_.Assemble();
s0chi_.Finalize();
ofstream ofsS0("s0_const_initial.mat");
s0chi_.SpMat().Print(ofsS0);
a0_.Assemble();
a0_.Finalize();
}
else if (tdChi_ && newTime_ && !nonLinear_)
{
chiNLCoef_->SetTemp(T0_);
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
ofstream ofsS0("s0_lin_initial.mat");
s0chi_.SpMat().Print(ofsS0);
a0_.Update();
a0_.Assemble(0);
a0_.Finalize(0);
}
if ((tdQ_ && newTime_) || first_)
{
cout << "Assembling Q" << endl;
QCoef_.SetQPerp(Q_perp);
QCoef_.SetTime(t_ + dt_);
Qs_.Assemble();
Qs_.ParallelAssemble(RHS_);
cout << "Norm of Q: " << Qs_.Norml2() << endl;
}
first_ = false;
newTime_ = false;
newTimeStep_ = false;
}
void ImplicitDiffOp::Mult(const Vector &dT, Vector &Q) const
{
dT_.Distribute(dT);
add(T0_, dt_, dT_, T1_);
if (tdChi_ && nonLinear_)
{
chiNLCoef_->SetTemp(T1_);
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
}
else
{
cout << "Well this is a surprise..." << endl;
}
m0cp_.Mult(dT_, Q_);
s0chi_.AddMult(T1_, Q_);
Q_.ParallelAssemble(Q);
Q.SetSubVector(ess_bdr_tdofs_, 0.0);
}
Operator & ImplicitDiffOp::GetGradient(const Vector &dT) const
{
if (tdChi_)
{
if (!nonLinear_)
{
chiNLCoef_->SetTemp(T0_);
}
else
{
dT_.Distribute(dT);
add(T0_, dt_, dT_, T1_);
chiNLCoef_->SetTemp(T1_);
dChiNLCoef_->SetTemp(T1_);
gradTCoef_.SetGridFunction(&T1_);
}
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
a0_.Update();
a0_.Assemble(0);
a0_.Finalize(0);
}
if (!nonLinear_)
{
s0chi_.Mult(T0_, rhs_);
rhs_ -= Qs_;
rhs_ *= -1.0;
}
else
{
rhs_ = Qs_;
}
dTdt_.ProjectBdrCoefficient(*bdrCoef_, ess_bdr_attr_);
a0_.FormLinearSystem(ess_bdr_tdofs_, dTdt_, rhs_, A_, SOL_, RHS_);
return A_;
}
Solver & ImplicitDiffOp::GetGradientSolver() const
{
if (!nonLinear_)
{
Operator & A_op = this->GetGradient(T0_); // T0_ will be ignored
HypreParMatrix & A_hyp = dynamic_cast<HypreParMatrix &>(A_op);
if (tdChi_)
{
delete AInv_; AInv_ = NULL;
delete APrecond_; APrecond_ = NULL;
}
if ( AInv_ == NULL )
{
// A_hyp.Print("A.mat");
HyprePCG * AInv_pcg = NULL;
cout << "Building PCG" << endl;
AInv_pcg = new HyprePCG(A_hyp);
AInv_pcg->SetTol(1e-12);
AInv_pcg->SetMaxIter(200);
AInv_pcg->SetPrintLevel(0);
if ( APrecond_ == NULL )
{
cout << "Building AMG" << endl;
APrecond_ = new HypreBoomerAMG(A_hyp);
APrecond_->SetPrintLevel(0);
AInv_pcg->SetPreconditioner(*APrecond_);
}
AInv_ = AInv_pcg;
}
}
else
{
if (AInv_ == NULL)
{
/*
HypreSmoother *J_hypreSmoother = new HypreSmoother;
J_hypreSmoother->SetType(HypreSmoother::l1Jacobi);
J_hypreSmoother->SetPositiveDiagonal(true);
JPrecond_ = J_hypreSmoother;
GMRESSolver * AInv_gmres = NULL;
cout << "Building GMRES" << endl;
AInv_gmres = new GMRESSolver(T0_.ParFESpace()->GetComm());
AInv_gmres->SetRelTol(1e-12);
AInv_gmres->SetAbsTol(0.0);
AInv_gmres->SetMaxIter(20000);
AInv_gmres->SetPrintLevel(2);
AInv_gmres->SetPreconditioner(*JPrecond_);
AInv_ = AInv_gmres;
*/
HypreGMRES * AInv_gmres = NULL;
cout << "Building HypreGMRES" << endl;
AInv_gmres = new HypreGMRES(T0_.ParFESpace()->GetComm());
AInv_gmres->SetTol(1e-12);
AInv_gmres->SetMaxIter(200);
AInv_gmres->SetPrintLevel(2);
if ( APrecond_ == NULL )
{
cout << "Building AMG" << endl;
APrecond_ = new HypreBoomerAMG();
APrecond_->SetPrintLevel(0);
AInv_gmres->SetPreconditioner(*APrecond_);
}
AInv_ = AInv_gmres;
}
}
return *AInv_;
}
} // namespace thermal
void
MatrixInverseCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K.Invert();
}
void
ScaledMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K *= a_;
}
} // namespace mfem
#endif // MFEM_USE_MPI
+555
View File
@@ -0,0 +1,555 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#ifndef MFEM_FOURIER_HYBRID_SOLVER
#define MFEM_FOURIER_HYBRID_SOLVER
#include "../common/pfem_extras.hpp"
#ifdef MFEM_USE_MPI
#include <memory>
#include <iostream>
#include <fstream>
namespace mfem
{
class NLCoefficient
{
protected:
NLCoefficient() : T_(NULL) {};
NLCoefficient(GridFunctionCoefficient & T) : T_(&T) {};
GridFunctionCoefficient * T_;
public:
virtual void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
};
class ChiParaCoef : public MatrixCoefficient, public NLCoefficient
{
private:
MatrixCoefficient * bbT_;
double chi_para_;
bool nonlin_;
public:
ChiParaCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_para, bool nonlin = false)
: MatrixCoefficient(2), NLCoefficient(T), bbT_(&bbT), //T_(&T),
chi_para_(chi_para), nonlin_(nonlin)
{}
//void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
class ChiPerpCoef : public MatrixCoefficient, public NLCoefficient
{
private:
MatrixCoefficient * bbT_;
// GridFunctionCoefficient * T_;
double chi_perp_;
bool nonlin_;
public:
ChiPerpCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_perp, bool nonlin = false)
: MatrixCoefficient(2), NLCoefficient(T), bbT_(&bbT),// T_(&T),
chi_perp_(chi_perp), nonlin_(nonlin)
{}
// void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
class ChiCoef : public MatrixSumCoefficient, public NLCoefficient
{
private:
ChiPerpCoef * chiPerpCoef_;
ChiParaCoef * chiParaCoef_;
public:
ChiCoef(ChiPerpCoef & chiPerp, ChiParaCoef & chiPara)
: MatrixSumCoefficient(chiPerp, chiPara),
chiPerpCoef_(&chiPerp), chiParaCoef_(&chiPara) {}
void SetTemp(GridFunction & T)
{
NLCoefficient::SetTemp(T);
chiPerpCoef_->SetTemp(T);
chiParaCoef_->SetTemp(T);
}
using MatrixSumCoefficient::Eval;
};
class dChiParaCoef : public MatrixCoefficient, public NLCoefficient
{
private:
MatrixCoefficient * bbT_;
// GridFunctionCoefficient * T_;
double chi_para_;
public:
dChiParaCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_para)
: MatrixCoefficient(2), NLCoefficient(T), bbT_(&bbT), //T_(&T),
chi_para_(chi_para)
{}
// void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
class dChiCoef : public MatrixCoefficient, public NLCoefficient
{
private:
MatrixCoefficient * bbT_;
// GridFunctionCoefficient * T_;
double chi_perp_;
double chi_para_;
public:
dChiCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_perp, double chi_para)
: MatrixCoefficient(2), NLCoefficient(T), bbT_(&bbT),// T_(&T),
chi_perp_(chi_perp), chi_para_(chi_para)
{}
// void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
class ChiInvParaCoef : public MatrixCoefficient, public NLCoefficient
{
private:
MatrixCoefficient * bbT_;
// GridFunctionCoefficient * T_;
double chi_para_;
bool nonlin_;
public:
ChiInvParaCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_para, bool nonlin = false)
: MatrixCoefficient(2), NLCoefficient(T), bbT_(&bbT),// T_(&T),
chi_para_(chi_para), nonlin_(nonlin)
{}
// void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
class ChiInvPerpCoef : public MatrixCoefficient, public NLCoefficient
{
private:
MatrixCoefficient * bbT_;
// GridFunctionCoefficient * T_;
double chi_perp_;
bool nonlin_;
public:
ChiInvPerpCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_perp, bool nonlin = false)
: MatrixCoefficient(2), NLCoefficient(T), bbT_(&bbT),// T_(&T),
chi_perp_(chi_perp), nonlin_(nonlin)
{}
// void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
class ChiInvCoef : public MatrixSumCoefficient, public NLCoefficient
{
private:
ChiInvPerpCoef * chiInvPerpCoef_;
ChiInvParaCoef * chiInvParaCoef_;
public:
ChiInvCoef(ChiInvPerpCoef & chiInvPerp, ChiInvParaCoef & chiInvPara)
: MatrixSumCoefficient(chiInvPerp, chiInvPara),
chiInvPerpCoef_(&chiInvPerp), chiInvParaCoef_(&chiInvPara) {}
void SetTemp(GridFunction & T)
{
NLCoefficient::SetTemp(T);
chiInvPerpCoef_->SetTemp(T);
chiInvParaCoef_->SetTemp(T);
}
};
class QParaCoef : public Coefficient
{
private:
Coefficient * Q_;
GridFunctionCoefficient * Q_perp_;
public:
QParaCoef(Coefficient & Q, GridFunctionCoefficient &Q_perp)
: Q_(&Q), Q_perp_(&Q_perp)
{}
void SetQPerp(GridFunction & Q) { Q_perp_->SetGridFunction(&Q); }
double Eval(ElementTransformation &T, const IntegrationPoint &ip)
{ return Q_->Eval(T, ip) - Q_perp_->Eval(T, ip); }
};
namespace thermal
{
class ImplicitDiffOp : public Operator
{
public:
ImplicitDiffOp(ParFiniteElementSpace & H1_FESpace,
Coefficient & dTdtBdr, bool tdBdr,
Array<int> & bdr_attr,
Coefficient & heatCap, bool tdCp,
MatrixCoefficient & chi, bool tdChi,
MatrixCoefficient & dchi, bool tdDChi,
Coefficient & heatSource, bool tdQ,
bool nonlinear = false);
~ImplicitDiffOp();
void SetState(ParGridFunction & T, ParGridFunction & Q_perp,
double t, double dt);
void Mult(const Vector &x, Vector &y) const;
Operator & GetGradient(const Vector &x) const;
Solver & GetGradientSolver() const;
const Vector & GetRHS() const { return RHS_; }
private:
bool first_;
bool tdBdr_;
bool tdCp_;
bool tdChi_;
bool tdDChi_;
bool tdQ_;
bool nonLinear_;
bool newTime_;
bool newTimeStep_;
double t_;
double dt_;
Array<int> & ess_bdr_attr_;
Array<int> ess_bdr_tdofs_;
Coefficient * bdrCoef_;
Coefficient * cpCoef_;
MatrixCoefficient * chiCoef_;
MatrixCoefficient * dChiCoef_;
NLCoefficient * chiNLCoef_;
NLCoefficient * dChiNLCoef_;
// Coefficient * QCoef_;
GridFunctionCoefficient QPerpCoef_;
QParaCoef QCoef_;
ScalarMatrixProductCoefficient dtChiCoef_;
mutable ParGridFunction T0_;
mutable ParGridFunction T1_;
mutable ParGridFunction dT_;
mutable GradientGridFunctionCoefficient gradTCoef_;
ScalarVectorProductCoefficient dtGradTCoef_;
MatVecCoefficient dtdChiGradTCoef_;
ParBilinearForm m0cp_;
mutable ParBilinearForm s0chi_;
mutable ParBilinearForm a0_;
mutable HypreParMatrix A_;
mutable ParGridFunction dTdt_;
mutable ParLinearForm Q_;
mutable ParLinearForm Qs_;
mutable ParLinearForm rhs_;
mutable Vector SOL_;
mutable Vector RHS_;
// Vector RHS0_; // Dummy RHS vector which hase length zero
mutable Solver * AInv_;
mutable HypreBoomerAMG * APrecond_;
};
/**
The thermal diffusion equation can be written:
dcT/dt = Div (chi Grad T) + Q_s
where
T is the temperature.
Div is the divergence operator,
grad is the gradient operator,
chi is the thermal conductivity tensor,
c is the heat capacity,
Q_s is the heat source
Class ThermalDiffusionTDO represents the right-hand side of the above
system of ODEs.
f(t, T) = -M_0(c)^{-1}(S_0(chi)T - M_0 Q_s)
where
M_0(c) is an H_1 mass matrix
S_0(chi) is the diffusion operator
The implicit solve method will solve
(M_0(c)+dt S_0(chi))k = -S_0(chi)T + M_0 Q_s
*/
class HybridThermalDiffusionTDO : public TimeDependentOperator
{
public:
HybridThermalDiffusionTDO(ParFiniteElementSpace &H1_FES,
ParFiniteElementSpace &HCurl_FES,
ParFiniteElementSpace &HDiv_FES,
ParFiniteElementSpace &L2_FES,
VectorCoefficient & dqdtBdr,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
double chi_perp,
double chi_para,
int prob,
int coef_type,
VectorCoefficient & UnitB,
Coefficient & c, bool td_c,
Coefficient & Q, bool td_Q);
void SetTime(const double time);
/** @brief Perform the action of the operator: @a q = f(@a y, t), where
q solves the algebraic equation F(@a y, q, t) = G(@a y, t) and t is the
current time. */
virtual void Mult(const Vector &y, Vector &q) const;
/** @brief Solve the equation: @a q = f(@a y + @a dt @a q, t), for the
unknown @a q at the current time t.
For general F and G, the equation for @a q becomes:
F(@a y + @a dt @a q, @a q, t) = G(@a y + @a dt @a q, t).
The input vector @a y corresponds to time index (or cycle) n, while the
currently set time, #t, and the result vector @a q correspond to time
index n+1. The time step @a dt corresponds to the time interval between
cycles n and n+1.
This method allows for the abstract implementation of some time
integration methods, including diagonal implicit Runge-Kutta (DIRK)
methods and the backward Euler method in particular.
If not re-implemented, this method simply generates an error. */
virtual void ImplicitSolve(const double dt, const Vector &y, Vector &q);
virtual ~HybridThermalDiffusionTDO();
void SetVisItDC(VisItDataCollection & visit_dc);
void GetParaFluxFromFlux(const ParGridFunction &q, ParGridFunction & q_para);
void GetPerpFluxFromFlux(const ParGridFunction &q, ParGridFunction & q_perp);
void GetParaFluxFromTemp(const ParGridFunction &T, ParGridFunction & q_para);
void GetPerpFluxFromTemp(const ParGridFunction &T, ParGridFunction & q_perp);
private:
void init();
void initA(double dt);
void initImplicitSolve();
bool init_;
bool newTime_;
bool nonLinear_;
bool testGradient_;
int dim_;
int tsize_;
int qsize_;
mutable int multCount_;
int solveCount_;
mutable ParGridFunction T_;
mutable ParGridFunction dT_;
// mutable ParGridFunction q_;
mutable ParGridFunction Q_perp_;
GridFunctionCoefficient TCoef_;
VectorCoefficient * unitBCoef_;
OuterProductCoefficient bbTCoef_;
IdentityMatrixCoefficient ICoef_;
MatrixSumCoefficient PPerpCoef_;
ChiPerpCoef chiPerpCoef_;
ChiParaCoef chiParaCoef_;
ChiCoef chiCoef_;
dChiCoef dChiCoef_;
dChiParaCoef dChiParaCoef_;
ChiInvPerpCoef chiInvPerpCoef_;
ChiInvParaCoef chiInvParaCoef_;
ChiInvCoef chiInvCoef_;
ParFiniteElementSpace * H1_FESpace_;
ParFiniteElementSpace * HCurl_FESpace_;
ParFiniteElementSpace * HDiv_FESpace_;
ParFiniteElementSpace * L2_FESpace_;
ParBilinearForm * m2_;
ParBilinearForm * mPara_;
ParBilinearForm * mPerp_;
ParBilinearForm * sC_;
ParMixedBilinearForm * dC_;
ParBilinearForm * a_;
ParMixedBilinearForm * gPara_;
ParMixedBilinearForm * gPerp_;
ParDiscreteLinearOperator * Div_;
ParDiscreteLinearOperator * Grad_;
ParGridFunction * dqdt_gf_;
ParGridFunction * Qs_;
mutable HypreParMatrix M2_;
mutable HyprePCG * M2Inv_;
mutable HypreDiagScale * M2Diag_;
HypreParMatrix A_;
HyprePCG * AInv_;
HypreSolver * APrecond_;
// HypreParVector * T_;
mutable ParGridFunction q_;
// mutable ParGridFunction u_;
mutable ParGridFunction dqdt_;
mutable ParGridFunction dqdt_perp_;
mutable ParGridFunction dqdt_para_;
mutable ParGridFunction dqdt_from_T_;
mutable ParGridFunction dqdt_para_from_T_;
mutable ParGridFunction q1_perp_;
mutable ParLinearForm dqdt_perp_dual_;
mutable ParLinearForm dqdt_para_dual_;
// mutable ParGridFunction dudt_;
mutable Vector X_;
mutable Vector RHS_;
mutable Vector rhs_;
mutable Vector dQs_;
// mutable Vector tmp_;
Array<int> * bdr_attr_;
Array<int> ess_bdr_tdofs_;
VectorCoefficient * dqdtBdrCoef_;
bool tdQ_;
bool tdC_;
bool tdK_;
Coefficient * QCoef_;
Coefficient * CCoef_;
// Coefficient * kCoef_;
// MatrixCoefficient * KCoef_;
Coefficient * CInvCoef_;
// Coefficient * kInvCoef_;
// MatrixCoefficient * KInvCoef_;
Coefficient * dtCInvCoef_;
ImplicitDiffOp impOp_;
NewtonSolver newton_;
};
} // namespace thermal
class InverseCoefficient : public Coefficient
{
public:
InverseCoefficient(Coefficient & c) : c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return 1.0 / c_->Eval(T, ip); }
private:
Coefficient * c_;
};
class MatrixInverseCoefficient :public MatrixCoefficient
{
public:
MatrixInverseCoefficient(MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
MatrixCoefficient * M_;
};
class ScaledCoefficient : public Coefficient
{
public:
ScaledCoefficient(double a, Coefficient & c) : a_(a), c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return a_ * c_->Eval(T, ip); }
private:
double a_;
Coefficient * c_;
};
class ScaledMatrixCoefficient :public MatrixCoefficient
{
public:
ScaledMatrixCoefficient(double a, MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), a_(a), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
double a_;
MatrixCoefficient * M_;
};
} // namespace mfem
#endif // MFEM_USE_MPI
#endif // MFEM_FOURIER_HYBRID_SOLVER
+579
View File
@@ -0,0 +1,579 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
//
// -----------------------------------------------------
// Fourier Miniapp: Thermal Diffusion
// -----------------------------------------------------
//
// This miniapp solves a time dependent heat equation.
//
#include "fourier_nl_solver.hpp"
#include <memory>
#include <iostream>
#include <fstream>
using namespace std;
using namespace mfem;
using namespace mfem::thermal;
void display_banner(ostream & os);
static int prob_ = 1;
static int unit_vec_type_ = 1;
static bool non_linear_ = false;
static double theta_ = M_PI/6.0;
static double nl_exp_ = 2.5;
static double chi_perp_ = 1.0;
static double chi_para_max_ = 1.0;
static double chi_para_min_ = 1.0;
double TFunc(const Vector &x, double t)
{
if ( prob_ % 2 == 1)
{
double e = exp(-2.0 * M_PI * M_PI * t);
return sin(M_PI * x[0]) * sin(M_PI * x[1]) * (1.0 - e);
}
else
{
double a = 0.4;
double b = 0.8;
double r = pow(x[0] / a, 2) + pow(x[1] / b, 2);
double e = exp(-0.25 * t * M_PI * M_PI / (a * b) );
return cos(0.5 * M_PI * sqrt(r)) * (1.0 - e);
}
}
double QFunc(const Vector &x, double t)
{
if ( prob_ % 2 == 1)
{
if (unit_vec_type_ == 1)
return 2.0 * chi_perp_ * M_PI * M_PI *
sin(M_PI * x[0]) * sin(M_PI * x[1]);
else
{
double chi_ratio = (nl_exp_ > 0.0) ?
pow(chi_para_min_ / chi_para_max_, 1.0 / nl_exp_) : 1.0;
double cx = cos(M_PI * x[0]);
double sx = sin(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sy = sin(M_PI * x[1]);
double ct = cos(theta_);
double st = sin(theta_);
double s2t = sin(2.0 * theta_);
double u = sx * sy;
double T = chi_ratio + (1.0 - chi_ratio) * u;
return M_PI * M_PI * (chi_perp_ * (u + cx * cy * s2t) +
chi_para_max_ * (u - cx * cy * s2t) * pow(T, nl_exp_) +
chi_para_max_ * nl_exp_ * (1.0 - chi_ratio) *
(u * u - sx * sx * st * st - sy * sy * ct * ct -
u * cx * cy * s2t) * pow(T, nl_exp_ - 1.0) );
}
}
else
{
double a = 0.4;
double b = 0.8;
double r = pow(x[0] / a, 2) + pow(x[1] / b, 2);
double r4 = pow(x[0] / (a * a), 2) + pow(x[1] / (b * b), 2);
double e = exp(-0.25 * t * M_PI * M_PI / (a * b) );
if ( r == 0.0 )
return 0.25 * M_PI * M_PI *
( chi_perp_ * (1.0 - e) * ( pow(a, -2) + pow(b, -2) ) +
e / (a * b));
return 0.25 * M_PI * M_PI *
( e / (a * b) + chi_perp_ * (r4 / r) * (1.0 - e)) *
cos(0.5 * M_PI * sqrt(r)) +
0.5 * M_PI * chi_perp_ * pow(a * b, -2) * (x * x) * (1.0 - e) *
sin(0.5 * M_PI * sqrt(r)) / pow(r, 1.5);
}
}
/*
void ChiFunc(const Vector &x, DenseMatrix &M)
{
M.SetSize(2);
double cx = cos(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sx = sin(M_PI * x[0]);
double sy = sin(M_PI * x[1]);
double den = cx * cx * sy * sy + sx * sx * cy * cy;
M(0,0) = chi_ratio_ * sx * sx * cy * cy + sy * sy * cx * cx;
M(1,1) = chi_ratio_ * sy * sy * cx * cx + sx * sx * cy * cy;
M(0,1) = (1.0 - chi_ratio_) * cx * cy * sx * sy;
M(1,0) = M(0,1);
M *= 1.0 / den;
}
*/
int main(int argc, char *argv[])
{
// 1. Initialize MPI.
MPI_Session mpi(argc, argv);
int myid = mpi.WorldRank();
// print the cool banner
if (mpi.Root()) { display_banner(cout); }
// 2. Parse command-line options.
int n = -1;
int order = 1;
int irOrder = -1;
int el_type = Element::QUADRILATERAL;
int ode_solver_type = 1;
int coef_type = 0;
int vis_steps = 1;
double dt = -1.0;
double t_final = 5.0;
double tol = 1e-4;
const char *basename = "Fourier";
const char *mesh_file = "";
bool zero_start = true;
bool static_cond = false;
bool gfprint = true;
bool visit = true;
bool visualization = true;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&n, "-n", "--num-elems-1d",
"Number of elements in x and y directions. "
"Total number of elements is n^2.");
args.AddOption(&prob_, "-p", "--problem",
"Specify problem type: 1 - Square, 2 - Ellipse.");
args.AddOption(&unit_vec_type_, "-u", "--unit-vec-type",
"Specify B field unit vector type: \n"
" 1 - Square, 2 - Ellipse,\n"
" 3 - Constant (angle theta).");
args.AddOption(&coef_type, "-c", "--coef",
"Specify diffusion coefficient type: "
"0 - Constant, 1 - Linearized, 2 - Non-Linear.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree) or -1 for"
" isoparametric space.");
args.AddOption(&irOrder, "-iro", "--int-rule-order",
"Integration Rule Order.");
args.AddOption(&chi_perp_, "-chi-perp", "--chi-perpendicular",
"Chi_perp.");
args.AddOption(&chi_para_max_, "-chi-max", "--chi-para-max",
"Maximum value of chi along field lines.");
args.AddOption(&chi_para_min_, "-chi-min", "--chi-para-min",
"Minimum value of chi along field lines.");
args.AddOption(&dt, "-dt", "--time-step",
"Time step.");
args.AddOption(&t_final, "-tf", "--final-time",
"Final Time.");
args.AddOption(&tol, "-tol", "--tolerance",
"Tolerance used to determine convergence to steady state.");
args.AddOption(&el_type, "-e", "--element-type",
"Element type: 2-Triangle, 3-Quadrilateral.");
args.AddOption(&ode_solver_type, "-s", "--ode-solver",
"ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3\n\t."
"\t 22 - Mid-Point, 23 - SDIRK23, 34 - SDIRK34.");
args.AddOption(&zero_start, "-z", "--zero-start", "-no-z",
"--no-zero-start",
"Initial guess of zero or exact solution.");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&gfprint, "-print", "--print","-no-print","--no-print",
"Print results (grid functions) to disk.");
args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
"Enable or disable VisIt visualization.");
args.AddOption(&vis_steps, "-vs", "--visualization-steps",
"Visualize every n-th timestep.");
args.AddOption(&basename, "-k", "--outputfilename",
"Name of the visit dump files");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
if (irOrder < 0)
{
irOrder = std::max(4, 2 * order - 2);
}
non_linear_ = coef_type > 0;
// 3. Construct a (serial) mesh of the given size on all processors. We
// can handle triangular and quadrilateral surface meshes with the
// same code.
Mesh *mesh = (n > 0) ?
new Mesh(n, n, (Element::Type)el_type, 1) :
new Mesh(mesh_file, 1, 1);
int dim = mesh->Dimension();
// 4. This step is no longer needed
// 5. Define a parallel mesh by a partitioning of the serial mesh. Once the
// parallel mesh is defined, the serial mesh can be deleted.
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
delete mesh;
// 7. Determine the list of true (i.e. parallel conforming) essential
// boundary dofs. In this example, the boundary conditions are defined
// by marking all the boundary attributes from the mesh as essential
// (Dirichlet) and converting them to a list of true dofs.
Array<int> ess_tdof_list;
Array<int> ess_bdr(0);
if (pmesh->bdr_attributes.Size())
{
ess_bdr.SetSize(pmesh->bdr_attributes.Max());
ess_bdr = 1;
}
// The following is required for mesh refinement
// mesh->EnsureNCMesh();
// 6. Define the ODE solver used for time integration. Several implicit
// methods are available, including singly diagonal implicit Runge-Kutta
// (SDIRK).
ODESolver *ode_solver;
switch (ode_solver_type)
{
// Implicit L-stable methods
case 1: ode_solver = new BackwardEulerSolver; break;
case 2: ode_solver = new SDIRK23Solver(2); break;
case 3: ode_solver = new SDIRK33Solver; break;
// Implicit A-stable methods (not L-stable)
case 22: ode_solver = new ImplicitMidpointSolver; break;
case 23: ode_solver = new SDIRK23Solver; break;
case 34: ode_solver = new SDIRK34Solver; break;
default:
if (mpi.Root())
{
cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
}
delete mesh;
return 3;
}
// 12. Define the parallel finite element spaces. We use:
//
// H(curl) for electric field,
// H(div) for magnetic flux,
// H(div) for thermal flux,
// H(grad)/H1 for electrostatic potential,
// L2 for temperature
// L2 contains discontinuous "cell-center" finite elements, type 2 is
// "positive"
L2_FECollection L2FEC0(0, dim);
L2_FECollection L2FEC(order-1, dim);
// RT contains Raviart-Thomas "face-centered" vector finite elements with
// continuous normal component.
RT_FECollection HDivFEC(order-1, dim);
// H1 contains continuous "node-centered" Lagrange finite elements.
H1_FECollection HGradFEC(order, dim);
ParFiniteElementSpace L2FESpace0(pmesh, &L2FEC0);
ParFiniteElementSpace L2FESpace(pmesh, &L2FEC);
ParFiniteElementSpace HDivFESpace(pmesh, &HDivFEC);
ParFiniteElementSpace HGradFESpace(pmesh, &HGradFEC);
// The terminology is TrueVSize is the unique (non-redundant) number of dofs
// HYPRE_Int glob_size_l2 = L2FESpace.GlobalTrueVSize();
// HYPRE_Int glob_size_rt = HDivFESpace.GlobalTrueVSize();
HYPRE_Int glob_size_h1 = HGradFESpace.GlobalTrueVSize();
if (mpi.Root())
{
cout << "Number of Temperature unknowns: " << glob_size_h1 << endl;
}
// int Vsize_l2 = L2FESpace.GetVSize();
// int Vsize_rt = HDivFESpace.GetVSize();
// int Vsize_h1 = HGradFESpace.GetVSize();
// grid functions E, B, T, F, P, and w which is the Joule heating
ParGridFunction T_gf(&HGradFESpace);
ParGridFunction dT_gf(&HGradFESpace);
ParGridFunction Qs_gf(&HGradFESpace);
ParGridFunction errorT(&L2FESpace0);
T_gf = 0.0;
dT_gf = 1.0;
// 13. Get the boundary conditions, set up the exact solution grid functions
// These VectorCoefficients have an Eval function. Note that e_exact and
// b_exact in this case are exact analytical solutions, taking a 3-vector
// point as input and returning a 3-vector field
FunctionCoefficient TCoef(TFunc);
ConstantCoefficient zeroCoef(0.0);
ConstantCoefficient SpecificHeatCoef(1.0);
// MatrixFunctionCoefficient ConductionCoef(2, ChiFunc);
FunctionCoefficient HeatSourceCoef(QFunc);
Qs_gf.ProjectCoefficient(HeatSourceCoef);
T_gf.GridFunction::ComputeElementL2Errors(TCoef, errorT);
// 14. Initialize the Diffusion operator, the GLVis visualization and print
// the initial energies.
ThermalDiffusionTDO oper(HGradFESpace,
zeroCoef, ess_bdr,
chi_perp_,
chi_para_min_,
chi_para_max_,
prob_,
unit_vec_type_,
coef_type,
SpecificHeatCoef, false,
// ConductionCoef, false,
HeatSourceCoef, false);
// This function initializes all the fields to zero or some provided IC
// oper.Init(F);
socketstream vis_T, vis_Q, vis_errT;
char vishost[] = "localhost";
int visport = 19916;
if (visualization)
{
// Make sure all ranks have sent their 'v' solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
vis_T.precision(8);
vis_Q.precision(8);
vis_errT.precision(8);
int Wx = 0, Wy = 0; // window position
int Ww = 350, Wh = 350; // window size
int offx = Ww+10;//, offy = Wh+45; // window offsets
miniapps::VisualizeField(vis_Q, vishost, visport,
Qs_gf, "Heat Soruce", Wx, Wy, Ww, Wh);
Wx += offx;
miniapps::VisualizeField(vis_T, vishost, visport,
T_gf, "Temperature", Wx, Wy, Ww, Wh);
Wx += offx;
miniapps::VisualizeField(vis_errT, vishost, visport,
errorT, "Error in T", Wx, Wy, Ww, Wh);
}
// VisIt visualization
VisItDataCollection visit_dc(basename, pmesh);
if ( visit )
{
visit_dc.RegisterField("T", &T_gf);
visit_dc.RegisterField("Qs", &Qs_gf);
visit_dc.RegisterField("L2 Error T", &errorT);
visit_dc.SetCycle(0);
visit_dc.SetTime(0.0);
visit_dc.Save();
}
ostringstream oss_errs;
oss_errs << "fourier_nl_errs"
<< "_p" << prob_ << "_c" << coef_type
<< "_e" << (int)floor(log10(chi_para_max_/chi_perp_));
if (n > 0) { oss_errs << "_n" << n; }
oss_errs << "_o" << order << ".dat";
ofstream ofs_errs;
if (myid == 0) { ofs_errs.open(oss_errs.str().c_str()); }
// 15. Perform time-integration (looping over the time iterations, ti, with a
// time-step dt). The object oper is the MagneticDiffusionOperator which
// has a Mult() method and an ImplicitSolve() method which are used by
// the time integrators.
ode_solver->Init(oper);
double t = 0.0;
double dt_courant = 0.0;
{
double h_min, h_max, kappa_min, kappa_max;
pmesh->GetCharacteristics(h_min, h_max, kappa_min, kappa_max);
dt_courant = 1.0 * h_min * h_min / chi_para_max_;
}
if (dt < 0.0)
{
dt = dt_courant;
}
if ( myid == 0 )
{
cout << "Using time step: " << dt
<< " (Courant " << dt_courant << ")" << endl;
}
int tsize = HGradFESpace.GetTrueVSize();
Vector T0(tsize), T1(tsize), dT(tsize);
T0 = 0.0; T1 = 0.0; dT = 0.0;
bool last_step = false;
for (int ti = 1; !last_step; ti++)
{
if (t + dt >= t_final - dt/2)
{
if (myid == 0)
{
cout << "Final Time Reached" << endl;
}
last_step = true;
}
// F is the vector of dofs, t is the current time, and dt is the time step
// to advance.
T0 = T1;
ode_solver->Step(T1, t, dt);
T_gf.Distribute(T1);
TCoef.SetTime(t);
T_gf.GridFunction::ComputeElementL2Errors(TCoef, errorT);
double l2_error_T = T_gf.ComputeL2Error(TCoef);
if ( myid == 0 )
{
ofs_errs << t << '\t' << l2_error_T << endl;
cout << t << '\t' << l2_error_T << endl;
}
add(1.0, T1, -1.0, T0, dT);
dT_gf.Distribute(dT);
double maxT = T_gf.ComputeMaxError(zeroCoef);
double maxDiff = dT_gf.ComputeMaxError(zeroCoef);
if ( !last_step )
{
if ( maxT == 0.0 )
{
last_step = (maxDiff < tol) ? true:false;
}
else if ( maxDiff/maxT < tol )
{
last_step = true;
}
if (last_step && myid == 0)
{
cout << "Converged to Steady State" << endl;
}
}
/*
if (debug == 1)
{
oper.Debug(basename,t);
}
*/
if (gfprint)
{
ostringstream T_name, mesh_name;
T_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "T." << setfill('0') << setw(6) << myid;
mesh_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "mesh." << setfill('0') << setw(6) << myid;
ofstream mesh_ofs(mesh_name.str().c_str());
mesh_ofs.precision(8);
pmesh->Print(mesh_ofs);
mesh_ofs.close();
ofstream T_ofs(T_name.str().c_str());
T_ofs.precision(8);
T_gf.Save(T_ofs);
T_ofs.close();
}
if (last_step || (ti % vis_steps) == 0)
{
// Make sure all ranks have sent their 'v' solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
if (visualization)
{
int Wx = 0, Wy = 0; // window position
int Ww = 350, Wh = 350; // window size
int offx = Ww+10;//, offy = Wh+45; // window offsets
Wx += offx;
miniapps::VisualizeField(vis_T, vishost, visport,
T_gf, "Temperature", Wx, Wy, Ww, Wh);
Wx += offx;
miniapps::VisualizeField(vis_errT, vishost, visport,
errorT, "Error in T", Wx, Wy, Ww, Wh);
}
if (visit)
{
visit_dc.SetCycle(ti);
visit_dc.SetTime(t);
visit_dc.Save();
}
}
}
if (visualization)
{
vis_T.close();
vis_errT.close();
}
if (myid == 0) { ofs_errs.close(); }
double loc_T_max = T1.Normlinf();
double T_max = -1.0;
MPI_Allreduce(&loc_T_max, &T_max, 1, MPI_DOUBLE, MPI_MAX,
MPI_COMM_WORLD);
double err1 = T_gf.ComputeL2Error(TCoef);
if (myid == 0)
{
cout << "L2 Error of Solution: " << err1 << endl;
cout << "Maximum Temperature: " << T_max << endl;
cout << "| chi_eff - 1 | = " << fabs(1.0/T_max - 1) << endl;
}
// 16. Free the used memory.
delete ode_solver;
delete pmesh;
return 0;
}
void display_banner(ostream & os)
{
os << "___________ .__ " << endl
<< "\\_ _____/___ __ _________|__| ___________ " << endl
<< " | __)/ _ \\| | \\_ __ \\ |/ __ \\_ __ \\" << endl
<< " | | ( <_> ) | /| | \\/ \\ ___/| | \\/" << endl
<< " \\__ | \\____/|____/ |__| |__|\\___ >__| " << endl
<< " \\/ \\/ " << endl
<< flush;
}
+507
View File
@@ -0,0 +1,507 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#include "fourier_nl_solver.hpp"
#ifdef MFEM_USE_MPI
using namespace std;
namespace mfem
{
using namespace miniapps;
void
UnitVectorField::Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip)
{
double x[2];
Vector transip(x, 2);
T.Transform(T.GetIntPoint(), transip);
V.SetSize(2);
if ( prob_ % 2 == 1 )
{
if (unit_vec_type_ == 1)
{
double cx = cos(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sx = sin(M_PI * x[0]);
double sy = sin(M_PI * x[1]);
V[0] = -sx * cy;
V[1] = sy * cx;
}
else
{
V[0] = cos(M_PI/6.0);
V[1] = sin(M_PI/6.0);
}
}
else
{
V[0] = -a_ * a_ * x[1];
V[1] = b_ * b_ * x[0];
}
double nrm = V.Norml2();
V *= (nrm > 1e-6 * min(a_,b_)) ? (1.0/nrm) : 0.0;
}
void ChiParaCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
if (type_ == 0)
{
K *= chi_max_;
}
else
{
K *= chi_min_ * pow(1.0 + gamma_ * T_->Eval(T, ip), 2.5);
}
}
void dChiCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
K *= 2.5 * chi_min_ * gamma_ * pow(1.0 + gamma_ * T_->Eval(T, ip), 1.5);
}
namespace thermal
{
ThermalDiffusionTDO::ThermalDiffusionTDO(
ParFiniteElementSpace &H1_FESpace,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
double chi_perp,
double chi_para_min,
double chi_para_max,
int prob,
int unit_vec_type,
int coef_type,
Coefficient & c, bool td_c,
Coefficient & Q, bool td_Q)
: TimeDependentOperator(H1_FESpace.GetTrueVSize(), 0.0),
init_(false),
nonLinear_(coef_type == 2),
testGradient_(false),
multCount_(0), solveCount_(0),
T_(&H1_FESpace),
TCoef_(&T_),
unitBCoef_(prob, unit_vec_type),
ICoef_(2),
bbTCoef_(unitBCoef_, unitBCoef_),
chiPerpCoef_(ICoef_, bbTCoef_, chi_perp, -chi_perp),
chiParaCoef_(bbTCoef_, TCoef_, coef_type, chi_para_min, chi_para_max),
chiCoef_(chiPerpCoef_, chiParaCoef_),
dChiCoef_(bbTCoef_, TCoef_, chi_para_min, chi_para_max),
impOp_(H1_FESpace,
dTdtBdr, false,
bdr_attr,
c, false,
chiCoef_, coef_type > 0,
dChiCoef_, coef_type > 0,
Q, false,
coef_type == 2),
newton_(H1_FESpace.GetComm())
{
this->init();
}
ThermalDiffusionTDO::~ThermalDiffusionTDO()
{
}
void
ThermalDiffusionTDO::init()
{
cout << "Entering TDO::Init" << endl;
if ( init_ ) { return; }
newton_.SetPrintLevel(2);
newton_.SetRelTol(1e-10);
newton_.SetAbsTol(0.0);
if ( nonLinear_ && testGradient_ )
{
Vector x(impOp_.Height());
Vector dx(impOp_.Height());
T_.Distribute(x);
cout << "GetTime " << this->GetTime() << endl;
impOp_.SetState(T_, this->GetTime(), 0.1);
cout << "init 0" << endl;
newton_.SetOperator(impOp_);
cout << "init 1" << endl;
cout << "init 2" << endl;
x.Randomize(1);
x.Print(cout);
dx.Randomize(2);
dx *= 0.01;
dx.Print(cout);
cout << "init 3" << endl;
double ratio = newton_.CheckGradient(x, dx);
cout << "CheckGradient returns: " << ratio << endl;
}
init_ = true;
cout << "Leaving TDO::Init" << endl;
}
void
ThermalDiffusionTDO::SetTime(const double time)
{
this->TimeDependentOperator::SetTime(time);
newTime_ = true;
}
void
ThermalDiffusionTDO::Mult(const Vector &T, Vector &dT_dt) const
{
MFEM_ABORT("ThermalDiffusionTDO::Mult should not be called");
}
void
ThermalDiffusionTDO::ImplicitSolve(const double dt,
const Vector &T, Vector &dT_dt)
{
dT_dt = 0.0;
T_.Distribute(T);
impOp_.SetState(T_, this->GetTime(), dt);
Solver & solver = impOp_.GetGradientSolver();
if (!nonLinear_)
{
solver.Mult(impOp_.GetRHS(), dT_dt);
}
else
{
newton_.SetOperator(impOp_);
newton_.SetSolver(solver);
newton_.Mult(impOp_.GetRHS(), dT_dt);
}
solveCount_++;
}
ImplicitDiffOp::ImplicitDiffOp(ParFiniteElementSpace & H1_FESpace,
Coefficient & dTdtBdr, bool tdBdr,
Array<int> & bdr_attr,
Coefficient & heatCap, bool tdCp,
ChiCoef & chi, bool tdChi,
dChiCoef & dchi, bool tdDChi,
Coefficient & heatSource, bool tdQ,
bool nonlinear)
: Operator(H1_FESpace.GetTrueVSize()),
first_(true),
tdBdr_(tdBdr),
tdCp_(tdCp),
tdChi_(tdChi),
tdDChi_(tdDChi),
tdQ_(tdQ),
nonLinear_(nonlinear),
newTime_(true),
newTimeStep_(true),
t_(0.0),
dt_(-1.0),
ess_bdr_attr_(bdr_attr),
bdrCoef_(&dTdtBdr),
cpCoef_(&heatCap),
chiCoef_(&chi),
dChiCoef_(&dchi),
QCoef_(&heatSource),
dtChiCoef_(1.0, *chiCoef_),
T0_(&H1_FESpace),
T1_(&H1_FESpace),
dT_(&H1_FESpace),
gradTCoef_(&T0_),
dtGradTCoef_(-1.0, gradTCoef_),
dtdChiGradTCoef_(*dChiCoef_, dtGradTCoef_),
m0cp_(&H1_FESpace),
s0chi_(&H1_FESpace),
a0_(&H1_FESpace),
dTdt_(&H1_FESpace),
Q_(&H1_FESpace),
Qs_(&H1_FESpace),
rhs_(&H1_FESpace),
RHS_(H1_FESpace.GetTrueVSize()),
// RHS0_(0),
AInv_(NULL),
APrecond_(NULL)
{
H1_FESpace.GetEssentialTrueDofs(ess_bdr_attr_, ess_bdr_tdofs_);
m0cp_.AddDomainIntegrator(new MassIntegrator(*cpCoef_));
s0chi_.AddDomainIntegrator(new DiffusionIntegrator(*chiCoef_));
a0_.AddDomainIntegrator(new MassIntegrator(*cpCoef_));
a0_.AddDomainIntegrator(new DiffusionIntegrator(dtChiCoef_));
if (nonLinear_)
{
a0_.AddDomainIntegrator(new MixedScalarWeakDivergenceIntegrator(
dtdChiGradTCoef_));
}
Qs_.AddDomainIntegrator(new DomainLFIntegrator(*QCoef_));
if (!tdQ_) { Qs_.Assemble(); }
}
ImplicitDiffOp::~ImplicitDiffOp()
{
delete AInv_;
delete APrecond_;
}
void ImplicitDiffOp::SetState(ParGridFunction & T, double t, double dt)
{
T0_ = T;
newTime_ = fabs(t - t_) > 0.0;
newTimeStep_= (fabs(1.0-dt/dt_)>1e-6);
t_ = newTime_ ? t : t_;
dt_ = newTimeStep_ ? dt : dt_;
if (tdBdr_ && (newTime_ || newTimeStep_))
{
bdrCoef_->SetTime(t_ + dt_);
}
if (newTimeStep_ || first_)
{
dtChiCoef_.SetAConst(dt_);
dtGradTCoef_.SetAConst(-dt_);
}
if ((tdCp_ && newTime_) || first_)
{
m0cp_.Update();
m0cp_.Assemble();
m0cp_.Finalize();
}
if (!tdChi_ && first_)
{
s0chi_.Assemble();
s0chi_.Finalize();
ofstream ofsS0("s0_const_initial.mat");
s0chi_.SpMat().Print(ofsS0);
a0_.Assemble();
a0_.Finalize();
}
else if (tdChi_ && newTime_ && !nonLinear_)
{
chiCoef_->SetTemp(T0_);
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
ofstream ofsS0("s0_lin_initial.mat");
s0chi_.SpMat().Print(ofsS0);
a0_.Update();
a0_.Assemble(0);
a0_.Finalize(0);
}
if ((tdQ_ && newTime_) || first_)
{
cout << "Assembling Q" << endl;
QCoef_->SetTime(t_ + dt_);
Qs_.Assemble();
Qs_.ParallelAssemble(RHS_);
cout << "Norm of Q: " << Qs_.Norml2() << endl;
}
first_ = false;
newTime_ = false;
newTimeStep_ = false;
}
void ImplicitDiffOp::Mult(const Vector &dT, Vector &Q) const
{
dT_.Distribute(dT);
add(T0_, dt_, dT_, T1_);
if (tdChi_ && nonLinear_)
{
chiCoef_->SetTemp(T1_);
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
}
else
{
cout << "Well this is a surprise..." << endl;
}
m0cp_.Mult(dT_, Q_);
s0chi_.AddMult(T1_, Q_);
Q_.ParallelAssemble(Q);
Q.SetSubVector(ess_bdr_tdofs_, 0.0);
}
Operator & ImplicitDiffOp::GetGradient(const Vector &dT) const
{
if (tdChi_)
{
if (!nonLinear_)
{
chiCoef_->SetTemp(T0_);
}
else
{
dT_.Distribute(dT);
add(T0_, dt_, dT_, T1_);
chiCoef_->SetTemp(T1_);
dChiCoef_->SetTemp(T1_);
gradTCoef_.SetGridFunction(&T1_);
}
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
a0_.Update();
a0_.Assemble(0);
a0_.Finalize(0);
}
if (!nonLinear_)
{
s0chi_.Mult(T0_, rhs_);
rhs_ -= Qs_;
rhs_ *= -1.0;
}
else
{
rhs_ = Qs_;
}
dTdt_.ProjectBdrCoefficient(*bdrCoef_, ess_bdr_attr_);
a0_.FormLinearSystem(ess_bdr_tdofs_, dTdt_, rhs_, A_, SOL_, RHS_);
return A_;
}
Solver & ImplicitDiffOp::GetGradientSolver() const
{
if (!nonLinear_)
{
Operator & A_op = this->GetGradient(T0_); // T0_ will be ignored
HypreParMatrix & A_hyp = dynamic_cast<HypreParMatrix &>(A_op);
if (tdChi_)
{
delete AInv_; AInv_ = NULL;
delete APrecond_; APrecond_ = NULL;
}
if ( AInv_ == NULL )
{
// A_hyp.Print("A.mat");
HyprePCG * AInv_pcg = NULL;
cout << "Building PCG" << endl;
AInv_pcg = new HyprePCG(A_hyp);
AInv_pcg->SetTol(1e-10);
AInv_pcg->SetMaxIter(200);
AInv_pcg->SetPrintLevel(0);
if ( APrecond_ == NULL )
{
cout << "Building AMG" << endl;
APrecond_ = new HypreBoomerAMG(A_hyp);
APrecond_->SetPrintLevel(0);
AInv_pcg->SetPreconditioner(*APrecond_);
}
AInv_ = AInv_pcg;
}
}
else
{
if (AInv_ == NULL)
{
/*
HypreSmoother *J_hypreSmoother = new HypreSmoother;
J_hypreSmoother->SetType(HypreSmoother::l1Jacobi);
J_hypreSmoother->SetPositiveDiagonal(true);
JPrecond_ = J_hypreSmoother;
GMRESSolver * AInv_gmres = NULL;
cout << "Building GMRES" << endl;
AInv_gmres = new GMRESSolver(T0_.ParFESpace()->GetComm());
AInv_gmres->SetRelTol(1e-12);
AInv_gmres->SetAbsTol(0.0);
AInv_gmres->SetMaxIter(20000);
AInv_gmres->SetPrintLevel(2);
AInv_gmres->SetPreconditioner(*JPrecond_);
AInv_ = AInv_gmres;
*/
HypreGMRES * AInv_gmres = NULL;
cout << "Building HypreGMRES" << endl;
AInv_gmres = new HypreGMRES(T0_.ParFESpace()->GetComm());
AInv_gmres->SetTol(1e-12);
AInv_gmres->SetMaxIter(200);
AInv_gmres->SetPrintLevel(2);
if ( APrecond_ == NULL )
{
cout << "Building AMG" << endl;
APrecond_ = new HypreBoomerAMG();
APrecond_->SetPrintLevel(0);
AInv_gmres->SetPreconditioner(*APrecond_);
}
AInv_ = AInv_gmres;
}
}
return *AInv_;
}
} // namespace thermal
void
MatrixInverseCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K.Invert();
}
void
ScaledMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K *= a_;
}
} // namespace mfem
#endif // MFEM_USE_MPI
+367
View File
@@ -0,0 +1,367 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#ifndef MFEM_FOURIER_NL_SOLVER
#define MFEM_FOURIER_NL_SOLVER
#include "../common/pfem_extras.hpp"
#ifdef MFEM_USE_MPI
#include <memory>
#include <iostream>
#include <fstream>
namespace mfem
{
class UnitVectorField : public VectorCoefficient
{
private:
int prob_;
int unit_vec_type_;
double a_;
double b_;
public:
UnitVectorField(int prob, int unit_vec_type, double a = 0.4, double b = 0.8)
: VectorCoefficient(2), prob_(prob), unit_vec_type_(unit_vec_type),
a_(a), b_(b) {}
void Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip);
};
/*
class ChiGridFuncCoef : public MatrixCoefficient
{
private:
double chi_min_ratio_;
double chi_max_ratio_;
int prob_;
const GridFunction & T_;
public:
ChiGridFuncCoef(const GridFunction & T,
double chi_min_ratio, double chi_max_ratio, int prob = 1)
: MatrixCoefficient(2),
chi_min_ratio_(chi_min_ratio),
chi_max_ratio_(chi_max_ratio),
prob_(prob),
T_(T) {}
// void SetTemp() { T_ = &T; }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
*/
class ChiParaCoef : public MatrixCoefficient
{
private:
MatrixCoefficient * bbT_;
GridFunctionCoefficient * T_;
int type_;
double chi_min_;
double chi_max_;
double gamma_;
public:
ChiParaCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T, int type,
double chi_min, double chi_max)
: MatrixCoefficient(2), bbT_(&bbT), T_(&T), type_(type),
chi_min_(chi_min), chi_max_(chi_max),
gamma_(pow(chi_max/chi_min, 0.4) - 1.0)
{}
void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
class ChiCoef : public MatrixSumCoefficient
{
private:
ChiParaCoef * chiParaCoef_;
public:
ChiCoef(MatrixCoefficient & chiPerp, ChiParaCoef & chiPara)
: MatrixSumCoefficient(chiPerp, chiPara), chiParaCoef_(&chiPara) {}
void SetTemp(GridFunction & T) { chiParaCoef_->SetTemp(T); }
};
class dChiCoef : public MatrixCoefficient
{
private:
MatrixCoefficient * bbT_;
GridFunctionCoefficient * T_;
double chi_min_;
double chi_max_;
double gamma_;
public:
dChiCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_min, double chi_max)
: MatrixCoefficient(2), bbT_(&bbT), T_(&T),
chi_min_(chi_min), chi_max_(chi_max),
gamma_(pow(chi_max/chi_min, 0.4) - 1.0)
{}
void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
namespace thermal
{
class ImplicitDiffOp : public Operator
{
public:
ImplicitDiffOp(ParFiniteElementSpace & H1_FESpace,
Coefficient & dTdtBdr, bool tdBdr,
Array<int> & bdr_attr,
Coefficient & heatCap, bool tdCp,
ChiCoef & chi, bool tdChi,
dChiCoef & dchi, bool tdDChi,
Coefficient & heatSource, bool tdQ,
bool nonlinear = false);
~ImplicitDiffOp();
void SetState(ParGridFunction & T, double t, double dt);
void Mult(const Vector &x, Vector &y) const;
Operator & GetGradient(const Vector &x) const;
Solver & GetGradientSolver() const;
const Vector & GetRHS() const { return RHS_; }
private:
bool first_;
bool tdBdr_;
bool tdCp_;
bool tdChi_;
bool tdDChi_;
bool tdQ_;
bool nonLinear_;
bool newTime_;
bool newTimeStep_;
double t_;
double dt_;
Array<int> & ess_bdr_attr_;
Array<int> ess_bdr_tdofs_;
Coefficient * bdrCoef_;
Coefficient * cpCoef_;
ChiCoef * chiCoef_;
dChiCoef * dChiCoef_;
Coefficient * QCoef_;
ScalarMatrixProductCoefficient dtChiCoef_;
mutable ParGridFunction T0_;
mutable ParGridFunction T1_;
mutable ParGridFunction dT_;
mutable GradientGridFunctionCoefficient gradTCoef_;
ScalarVectorProductCoefficient dtGradTCoef_;
MatVecCoefficient dtdChiGradTCoef_;
ParBilinearForm m0cp_;
mutable ParBilinearForm s0chi_;
mutable ParBilinearForm a0_;
mutable HypreParMatrix A_;
mutable ParGridFunction dTdt_;
mutable ParLinearForm Q_;
mutable ParLinearForm Qs_;
mutable ParLinearForm rhs_;
mutable Vector SOL_;
mutable Vector RHS_;
// Vector RHS0_; // Dummy RHS vector which hase length zero
mutable Solver * AInv_;
mutable HypreBoomerAMG * APrecond_;
};
/**
The thermal diffusion equation can be written:
dcT/dt = Div (chi Grad T) + Q_s
where
T is the temperature.
Div is the divergence operator,
grad is the gradient operator,
chi is the thermal conductivity tensor,
c is the heat capacity,
Q_s is the heat source
Class ThermalDiffusionTDO represents the right-hand side of the above
system of ODEs.
f(t, T) = -M_0(c)^{-1}(S_0(chi)T - M_0 Q_s)
where
M_0(c) is an H_1 mass matrix
S_0(chi) is the diffusion operator
The implicit solve method will solve
(M_0(c)+dt S_0(chi))k = -S_0(chi)T + M_0 Q_s
*/
class ThermalDiffusionTDO : public TimeDependentOperator
{
public:
ThermalDiffusionTDO(ParFiniteElementSpace &H1_FES,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
double chi_perp,
double chi_para_min,
double chi_para_max,
int prob,
int unit_vec_type,
int coef_type,
Coefficient & c, bool td_c,
Coefficient & Q, bool td_Q);
void SetTime(const double time);
/** @brief Perform the action of the operator: @a q = f(@a y, t), where
q solves the algebraic equation F(@a y, q, t) = G(@a y, t) and t is the
current time. */
virtual void Mult(const Vector &y, Vector &q) const;
/** @brief Solve the equation: @a q = f(@a y + @a dt @a q, t), for the
unknown @a q at the current time t.
For general F and G, the equation for @a q becomes:
F(@a y + @a dt @a q, @a q, t) = G(@a y + @a dt @a q, t).
The input vector @a y corresponds to time index (or cycle) n, while the
currently set time, #t, and the result vector @a q correspond to time
index n+1. The time step @a dt corresponds to the time interval between
cycles n and n+1.
This method allows for the abstract implementation of some time
integration methods, including diagonal implicit Runge-Kutta (DIRK)
methods and the backward Euler method in particular.
If not re-implemented, this method simply generates an error. */
virtual void ImplicitSolve(const double dt, const Vector &y, Vector &q);
virtual ~ThermalDiffusionTDO();
private:
void init();
bool init_;
bool newTime_;
bool nonLinear_;
bool testGradient_;
mutable int multCount_;
int solveCount_;
mutable ParGridFunction T_;
GridFunctionCoefficient TCoef_;
UnitVectorField unitBCoef_;
IdentityMatrixCoefficient ICoef_;
OuterProductCoefficient bbTCoef_;
MatrixSumCoefficient chiPerpCoef_;
ChiParaCoef chiParaCoef_;
ChiCoef chiCoef_;
dChiCoef dChiCoef_;
ImplicitDiffOp impOp_;
NewtonSolver newton_;
};
} // namespace thermal
class InverseCoefficient : public Coefficient
{
public:
InverseCoefficient(Coefficient & c) : c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return 1.0 / c_->Eval(T, ip); }
private:
Coefficient * c_;
};
class MatrixInverseCoefficient :public MatrixCoefficient
{
public:
MatrixInverseCoefficient(MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
MatrixCoefficient * M_;
};
class ScaledCoefficient : public Coefficient
{
public:
ScaledCoefficient(double a, Coefficient & c) : a_(a), c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return a_ * c_->Eval(T, ip); }
private:
double a_;
Coefficient * c_;
};
class ScaledMatrixCoefficient :public MatrixCoefficient
{
public:
ScaledMatrixCoefficient(double a, MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), a_(a), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
double a_;
MatrixCoefficient * M_;
};
} // namespace mfem
#endif // MFEM_USE_MPI
#endif // MFEM_FOURIER_NL_SOLVER
+387
View File
@@ -0,0 +1,387 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#include "fourier_solver.hpp"
#ifdef MFEM_USE_MPI
using namespace std;
namespace mfem
{
using namespace miniapps;
namespace thermal
{
ThermalDiffusionOperator::ThermalDiffusionOperator(
ParFiniteElementSpace &H1_FES,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
Coefficient & c, bool td_c,
Coefficient & k, bool td_k,
Coefficient & Q, bool td_Q)
: TimeDependentOperator(H1_FES.GetVSize(), 0.0),
init_(false), //initA_(false), initAInv_(false),
multCount_(0), solveCount_(0),
H1_FESpace_(&H1_FES),
mC_(NULL), sK_(NULL), a_(NULL), dTdt_gf_(NULL), Qs_(NULL),
MCInv_(NULL), MCDiag_(NULL),
AInv_(NULL), APrecond_(NULL),
rhs_(NULL),
bdr_attr_(&bdr_attr), ess_bdr_tdofs_(0), dTdtBdrCoef_(&dTdtBdr),
tdQ_(td_Q), tdC_(td_c), tdK_(td_k),
QCoef_(&Q), CCoef_(&c), kCoef_(&k), KCoef_(NULL),
// CInvCoef_(NULL), kInvCoef_(NULL), KInvCoef_(NULL)
dtkCoef_(NULL), dtKCoef_(NULL)
{
this->init();
}
ThermalDiffusionOperator::ThermalDiffusionOperator(
ParFiniteElementSpace &H1_FES,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
Coefficient & c, bool td_c,
MatrixCoefficient & K, bool td_k,
Coefficient & Q, bool td_Q)
: TimeDependentOperator(H1_FES.GetVSize(), 0.0),
init_(false),
multCount_(0), solveCount_(0),
H1_FESpace_(&H1_FES),
mC_(NULL), sK_(NULL), a_(NULL), dTdt_gf_(NULL), Qs_(NULL),
MCInv_(NULL), MCDiag_(NULL),
AInv_(NULL), APrecond_(NULL),
rhs_(NULL),
bdr_attr_(&bdr_attr), ess_bdr_tdofs_(0), dTdtBdrCoef_(&dTdtBdr),
tdQ_(td_Q), tdC_(td_c), tdK_(td_k),
QCoef_(&Q), CCoef_(&c), kCoef_(NULL), KCoef_(&K),
// CInvCoef_(NULL), kInvCoef_(NULL), KInvCoef_(NULL)
dtkCoef_(NULL), dtKCoef_(NULL)
{
this->init();
}
ThermalDiffusionOperator::~ThermalDiffusionOperator()
{
delete a_;
delete mC_;
delete sK_;
delete dTdt_gf_;
delete Qs_;
delete MCInv_;
delete MCDiag_;
delete AInv_;
delete APrecond_;
}
void
ThermalDiffusionOperator::init()
{
if ( init_ ) { return; }
if ( mC_ == NULL )
{
mC_ = new ParBilinearForm(H1_FESpace_);
mC_->AddDomainIntegrator(new MassIntegrator(*CCoef_));
mC_->Assemble();
}
if ( sK_ == NULL )
{
sK_ = new ParBilinearForm(H1_FESpace_);
if ( kCoef_ != NULL )
{
sK_->AddDomainIntegrator(new DiffusionIntegrator(*kCoef_));
}
else if ( KCoef_ != NULL )
{
sK_->AddDomainIntegrator(new DiffusionIntegrator(*KCoef_));
}
sK_->Assemble();
}
if ( dTdt_gf_ == NULL )
{
dTdt_gf_ = new ParGridFunction(H1_FESpace_);
}
if ( Qs_ == NULL && QCoef_ != NULL )
{
Qs_ = new ParLinearForm(H1_FESpace_);
Qs_->AddDomainIntegrator(new DomainLFIntegrator(*QCoef_));
Qs_->Assemble();
rhs_ = new Vector(Qs_->Size());
}
/*
CInvCoef_ = new InverseCoefficient(*CCoef_);
if ( kCoef_ != NULL ) kInvCoef_ = new InverseCoefficient(*kCoef_);
if ( KCoef_ != NULL ) KInvCoef_ = new MatrixInverseCoefficient(*KCoef_);
*/
H1_FESpace_->GetEssentialTrueDofs(*bdr_attr_, ess_bdr_tdofs_);
init_ = true;
}
void
ThermalDiffusionOperator::SetTime(const double time)
{
this->TimeDependentOperator::SetTime(time);
dTdtBdrCoef_->SetTime(t);
if ( tdQ_ )
{
QCoef_->SetTime(t);
Qs_->Assemble();
}
if ( tdC_ )
{
CCoef_->SetTime(t);
mC_->Assemble();
}
if ( tdK_ )
{
if ( kCoef_ != NULL ) { kCoef_->SetTime(t); }
if ( KCoef_ != NULL ) { KCoef_->SetTime(t); }
sK_->Assemble();
}
if ( ( tdC_ || tdK_ ) && a_ != NULL )
{
a_->Assemble();
}
newTime_ = true;
}
/*
void
ThermalDiffusionOperator::SetHeatSource(Coefficient & Q, bool time_dep)
{
if ( ownsQ_ )
{
delete QCoef_;
}
tdQ_ = time_dep;
QCoef_ = &Q;
}
void
ThermalDiffusionOperator::SetConductivityCoefficient(Coefficient & k,
bool time_dep)
{
if ( ownsK_ )
{
delete kCoef_;
delete KCoef_;
}
tdK_ = time_dep;
kCoef_ = &k;
KCoef_ = NULL;
}
void
ThermalDiffusionOperator::SetConductivityCoefficient(MatrixCoefficient & K,
bool time_dep)
{
if ( ownsK_ )
{
delete kCoef_;
delete KCoef_;
}
tdK_ = time_dep;
kCoef_ = NULL;
KCoef_ = &K;
}
void
ThermalDiffusionOperator::SetSpecificHeatCoefficient(Coefficient & c,
bool time_dep)
{
if ( ownsC_ )
{
delete CCoef_;
}
tdC_ = time_dep;
CCoef_ = &c;
}
*/
void
ThermalDiffusionOperator::initMult() const
{
if ( tdC_ || MCInv_ == NULL || MCDiag_ == NULL )
{
if ( MCInv_ == NULL )
{
MCInv_ = new HyprePCG(MC_);
MCInv_->SetTol(1e-12);
MCInv_->SetMaxIter(200);
MCInv_->SetPrintLevel(0);
}
else
{
MCInv_->SetOperator(MC_);
}
if ( MCDiag_ == NULL )
{
MCDiag_ = new HypreDiagScale(MC_);
MCInv_->SetPreconditioner(*MCDiag_);
}
else
{
MCDiag_->SetOperator(MC_);
}
}
}
void
ThermalDiffusionOperator::Mult(const Vector &T, Vector &dT_dt) const
{
dT_dt = 0.0;
sK_->Mult(T, *rhs_);
*rhs_ -= *Qs_;
rhs_->Neg();
dTdt_gf_->ProjectBdrCoefficient(*dTdtBdrCoef_, *bdr_attr_);
mC_->FormLinearSystem(ess_bdr_tdofs_, *dTdt_gf_, *rhs_, MC_, dTdt_, RHS_);
this->initMult();
MCInv_->Mult(RHS_, dTdt_);
mC_->RecoverFEMSolution(dTdt_, *rhs_, dT_dt);
multCount_++;
}
void
ThermalDiffusionOperator::initA(double dt)
{
if ( kCoef_ != NULL )
{
dtkCoef_ = new ScaledCoefficient(dt, *kCoef_);
}
else
{
dtKCoef_ = new ScaledMatrixCoefficient(dt, *KCoef_);
}
if ( a_ == NULL)
{
a_ = new ParBilinearForm(H1_FESpace_);
a_->AddDomainIntegrator(new MassIntegrator(*CCoef_));
if ( kCoef_ != NULL)
{
a_->AddDomainIntegrator(new DiffusionIntegrator(*dtkCoef_));
}
else
{
a_->AddDomainIntegrator(new DiffusionIntegrator(*dtKCoef_));
}
a_->Assemble();
}
}
void
ThermalDiffusionOperator::initImplicitSolve()
{
if ( tdC_ || tdK_ || AInv_ == NULL || APrecond_ == NULL )
{
if ( AInv_ == NULL )
{
AInv_ = new HyprePCG(A_);
AInv_->SetTol(1e-12);
AInv_->SetMaxIter(200);
AInv_->SetPrintLevel(0);
}
else
{
AInv_->SetOperator(A_);
}
if ( APrecond_ == NULL )
{
APrecond_ = new HypreBoomerAMG(A_);
APrecond_->SetPrintLevel(0);
AInv_->SetPreconditioner(*APrecond_);
}
else
{
APrecond_->SetOperator(A_);
}
}
}
void
ThermalDiffusionOperator::ImplicitSolve(const double dt,
const Vector &T, Vector &dT_dt)
{
dT_dt = 0.0;
// cout << "sK size: " << sK_->Width() << ", T size: " << T.Size() << ", rhs_ size: " << rhs_->Size() << endl;
ostringstream ossT; ossT << "T_" << solveCount_ << ".vec";
ofstream ofsT(ossT.str().c_str());
T.Print(ofsT);
ofsT.close();
sK_->Mult(T, *rhs_);
ofstream ofsrhs("rhs.vec");
rhs_->Print(ofsrhs);
ofstream ofsQ("Q.vec");
Qs_->Print(ofsQ);
*rhs_ -= *Qs_;
*rhs_ *= -1.0;
dTdt_gf_->ProjectBdrCoefficient(*dTdtBdrCoef_, *bdr_attr_);
this->initA(dt);
a_->FormLinearSystem(ess_bdr_tdofs_, *dTdt_gf_, *rhs_, A_, dTdt_, RHS_);
A_.Print("A.mat");
ofstream ofsB("b.vec");
RHS_.Print(ofsB);
this->initImplicitSolve();
AInv_->Mult(RHS_, dTdt_);
a_->RecoverFEMSolution(dTdt_, *rhs_, dT_dt);
solveCount_++;
}
} // namespace thermal
void
MatrixInverseCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K.Invert();
}
void
ScaledMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K *= a_;
}
} // namespace mfem
#endif // MFEM_USE_MPI
+241
View File
@@ -0,0 +1,241 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#ifndef MFEM_FOURIER_SOLVER
#define MFEM_FOURIER_SOLVER
#include "../common/pfem_extras.hpp"
#ifdef MFEM_USE_MPI
#include <memory>
#include <iostream>
#include <fstream>
namespace mfem
{
namespace thermal
{
/**
The thermal diffusion equation can be written:
dcT/dt = Div (sigma Grad T) + Q_s
where
T is the temperature.
Div is the divergence operator,
grad is the gradient operator,
sigma is the thermal conductivity,
c is the heat capacity,
Q_s is the heat source
Class ThermalDiffusionOperator represents the right-hand side of the above
system of ODEs.
f(t, T) = -M_0(c)^{-1}(S_0(sigma)T - M_0 Q_s)
where
M_0(c) is an H_1 mass matrix
S_0(sigma) is the diffusion operator
The implicit solve method will solve
(M_0(c)+dt S_0(sigma))k = -S_0(sigma)T + M_0 Q_s
*/
class ThermalDiffusionOperator : public TimeDependentOperator
{
public:
ThermalDiffusionOperator(ParFiniteElementSpace &H1_FES,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
Coefficient & c, bool td_c,
Coefficient & k, bool td_k,
Coefficient & Q, bool td_Q);
ThermalDiffusionOperator(ParFiniteElementSpace &H1_FES,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
Coefficient & c, bool td_c,
MatrixCoefficient & K, bool td_k,
Coefficient & Q, bool td_Q);
void SetTime(const double time);
/*
void SetHeatSource(Coefficient & Q, bool time_dep = false);
void SetConductivityCoefficient(Coefficient & k,
bool time_dep = false);
void SetConductivityCoefficient(MatrixCoefficient & K,
bool time_dep = false);
void SetSpecificHeatCoefficient(
bool time_dep = false);
*/
/** @brief Perform the action of the operator: @a q = f(@a y, t), where
q solves the algebraic equation F(@a y, q, t) = G(@a y, t) and t is the
current time. */
virtual void Mult(const Vector &y, Vector &q) const;
/** @brief Solve the equation: @a q = f(@a y + @a dt @a q, t), for the
unknown @a q at the current time t.
For general F and G, the equation for @a q becomes:
F(@a y + @a dt @a q, @a q, t) = G(@a y + @a dt @a q, t).
The input vector @a y corresponds to time index (or cycle) n, while the
currently set time, #t, and the result vector @a q correspond to time
index n+1. The time step @a dt corresponds to the time interval between
cycles n and n+1.
This method allows for the abstract implementation of some time
integration methods, including diagonal implicit Runge-Kutta (DIRK)
methods and the backward Euler method in particular.
If not re-implemented, this method simply generates an error. */
virtual void ImplicitSolve(const double dt, const Vector &y, Vector &q);
virtual ~ThermalDiffusionOperator();
private:
void init();
void initMult() const;
void initA(double dt);
void initImplicitSolve();
bool init_;
// bool initA_;
// bool initAInv_;
bool newTime_;
mutable int multCount_;
int solveCount_;
ParFiniteElementSpace * H1_FESpace_;
ParBilinearForm * mC_;
ParBilinearForm * sK_;
ParBilinearForm * a_;
ParGridFunction * dTdt_gf_;
ParLinearForm * Qs_;
mutable HypreParMatrix MC_;
mutable HyprePCG * MCInv_;
mutable HypreDiagScale * MCDiag_;
HypreParMatrix A_;
HyprePCG * AInv_;
HypreBoomerAMG * APrecond_;
// HypreParVector * T_;
mutable Vector dTdt_;
mutable Vector RHS_;
Vector * rhs_;
Array<int> * bdr_attr_;
Array<int> ess_bdr_tdofs_;
Coefficient * dTdtBdrCoef_;
bool tdQ_;
bool tdC_;
bool tdK_;
/*
bool ownsQ_;
bool ownsC_;
bool ownsK_;
*/
Coefficient * QCoef_;
Coefficient * CCoef_;
Coefficient * kCoef_;
MatrixCoefficient * KCoef_;
// Coefficient * CInvCoef_;
// Coefficient * kInvCoef_;
// MatrixCoefficient * KInvCoef_;
Coefficient * dtkCoef_;
MatrixCoefficient * dtKCoef_;
};
} // namespace thermal
class InverseCoefficient : public Coefficient
{
public:
InverseCoefficient(Coefficient & c) : c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return 1.0 / c_->Eval(T, ip); }
private:
Coefficient * c_;
};
class MatrixInverseCoefficient :public MatrixCoefficient
{
public:
MatrixInverseCoefficient(MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
MatrixCoefficient * M_;
};
class ScaledCoefficient : public Coefficient
{
public:
ScaledCoefficient(double a, Coefficient & c) : a_(a), c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return a_ * c_->Eval(T, ip); }
private:
double a_;
Coefficient * c_;
};
class ScaledMatrixCoefficient :public MatrixCoefficient
{
public:
ScaledMatrixCoefficient(double a, MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), a_(a), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
double a_;
MatrixCoefficient * M_;
};
} // namespace mfem
#endif // MFEM_USE_MPI
#endif // MFEM_FOURIER_SOLVER
+612
View File
@@ -0,0 +1,612 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
//
// -----------------------------------------------------
// Fourier Miniapp: Thermal Diffusion
// -----------------------------------------------------
//
// This miniapp solves a time dependent heat equation.
//
#include "fourier_vanEs_solver.hpp"
#include <memory>
#include <iostream>
#include <fstream>
using namespace std;
using namespace mfem;
using namespace mfem::thermal;
void display_banner(ostream & os);
static int prob_ = 1;
static int unit_vec_type_ = 1;
static bool non_linear_ = false;
static double alpha_ = NAN;
static double theta_ = NAN;
static double gamma_ = 10.0;
static double nl_perp_exp_ = -0.5;
static double nl_para_exp_ = 2.5;
static double chi_perp_ = 1.0;
static double chi_para_ = 1.0;
static double a_ = 0.15;
static double b_ = 0.85;
static double xc_ = 0.0;
static double yc_ = 0.0;
double TFunc(const Vector &x, double t)
{
switch (prob_)
{
case 1:
return x[0] * x[1] * pow(sin(M_PI * x[0]) * sin(M_PI * x[1]), gamma_);
case 2:
return 1.0 - pow(pow(x[0] - xc_, 2) + pow(x[1] - yc_, 2), 1.5);
case 3:
return 1.0 + (a_ * x[0] + b_ * x[1]) * pow(x[0] * x[0] + x[1] * x[1], 1.5);
case 4:
return 1.0 - pow(a_ * pow(x[0] * cos(theta_) + x[1] * sin(theta_), 2) +
b_ * pow(x[0] * sin(theta_) - x[1] * cos(theta_), 2), 1.5);
default:
return 0.0;
}
}
void UnitBFunc(const Vector &x, Vector &b)
{
switch (unit_vec_type_)
{
case 2:
{
b[0] = -x[1] + yc_;
b[1] = x[0] - xc_;
}
break;
case 3:
{
b[0] = -3.0 * a_ * x[0] * x[1] -
b_ * (x[0] * x[0] + 4.0 * x[1] * x[1]);
b[1] = a_ * (4.0 * x[0] * x[0] + x[1] * x[1]) + 3.0 * b_ * x[0] * x[1];
}
break;
case 4:
{
double ct = cos(theta_);
double st = sin(theta_);
double ctst = 0.5 * sin(2.0 * theta_);
b[0] = x[1] * (a_ * st * st + b_ * ct * ct) + (a_ - b_) * x[0] * ctst;
b[1] = -x[0] * (a_ * ct * ct + b_ * st * st) - (a_ - b_) * x[1] * ctst;
}
break;
default:
b[0] = cos(alpha_);
b[1] = sin(alpha_);
}
double nrm = b.Norml2();
if ( nrm > 0.0 ) { b /= nrm; }
}
double QFunc(const Vector &x, double t)
{
switch (prob_)
{
case 1:
{
double cx = cos(M_PI * x[0]);
double sx = sin(M_PI * x[0]);
double s2x = sin(2.0 * M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sy = sin(M_PI * x[1]);
double s2y = sin(2.0 * M_PI * x[1]);
double ca = cos(alpha_);
double sa = sin(alpha_);
double s2a = sin(2.0 * alpha_);
double chi_sc = chi_perp_ * sa * sa + chi_para_ * ca * ca;
double chi_cs = chi_perp_ * ca * ca + chi_para_ * sa * sa;
double chi_s2 = (chi_para_ - chi_perp_) * s2a;
double s2gcx = s2x + M_PI * x[0] * (gamma_ * cx * cx - 1.0);
double s2gcy = s2y + M_PI * x[1] * (gamma_ * cy * cy - 1.0);
double sgcx = sx + M_PI * x[0] * gamma_ * cx;
double sgcy = sy + M_PI * x[1] * gamma_ * cy;
return -1.0 * (M_PI * gamma_ * x[0] * chi_cs * s2gcy * sx * sx +
M_PI * gamma_ * x[1] * chi_sc * s2gcx * sy * sy +
chi_s2 * sgcx * sgcy * sx * sy) *
pow(sx * sy, gamma_ - 2.0);
}
case 2:
{
return 9.0 * chi_perp_ * sqrt(pow(x[0] - xc_, 2) + pow(x[1] - yc_, 2));
}
default:
return 0.0;
}
}
void shiftUnitSquare(const Vector &x, Vector &p)
{
p[0] = x[0] - 0.5;
p[1] = x[1] - 0.5;
}
int main(int argc, char *argv[])
{
// 1. Initialize MPI.
MPI_Session mpi(argc, argv);
int myid = mpi.WorldRank();
// print the cool banner
if (mpi.Root()) { display_banner(cout); }
// 2. Parse command-line options.
int n = -1;
int order = 1;
int irOrder = -1;
int el_type = Element::QUADRILATERAL;
int ode_solver_type = 1;
int coef_type = 0;
int vis_steps = 1;
double dt = -1.0;
double t_final = 5.0;
double tol = 1e-4;
const char *basename = "Fourier";
const char *mesh_file = "";
bool zero_start = true;
bool static_cond = false;
bool gfprint = true;
bool visit = true;
bool visualization = true;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&n, "-n", "--num-elems-1d",
"Number of elements in x and y directions. "
"Total number of elements is n^2.");
args.AddOption(&prob_, "-p", "--problem",
"Specify problem type:\n"
" 1 - section 4.1, 2 - section 4.2, 3 - section 4.3.");
// args.AddOption(&unit_vec_type_, "-u", "--unit-vec-type",
// "Specify B field unit vector type: \n"
// " 1 - Constant, 2 - ,\n"
// " 3 - Constant (angle theta).");
args.AddOption(&alpha_, "-alpha", "--constant-angle",
"Angle for constant B field (in degrees)");
args.AddOption(&xc_, "-xc", "--x-center",
"x coordinate of field center");
args.AddOption(&yc_, "-yc", "--y-center",
"y coordinate of field center");
args.AddOption(&coef_type, "-c", "--coef",
"Specify diffusion coefficient type: "
"0 - Constant, 1 - Linearized, 2 - Non-Linear.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree) or -1 for"
" isoparametric space.");
args.AddOption(&irOrder, "-iro", "--int-rule-order",
"Integration Rule Order.");
args.AddOption(&chi_perp_, "-chi-perp", "--chi-perpendicular",
"Chi_perp.");
args.AddOption(&chi_para_, "-chi-para", "--chi-parallel",
"Value of chi along field lines.");
// args.AddOption(&nonlin_chi, "-nl", "--nonlin-chi",
// "-no-nl", "--no-nonlin-chi",
// "Enable or disable Nonlinear Diffusion.");
args.AddOption(&dt, "-dt", "--time-step",
"Time step.");
args.AddOption(&t_final, "-tf", "--final-time",
"Final Time.");
args.AddOption(&tol, "-tol", "--tolerance",
"Tolerance used to determine convergence to steady state.");
args.AddOption(&el_type, "-e", "--element-type",
"Element type: 2-Triangle, 3-Quadrilateral.");
args.AddOption(&ode_solver_type, "-s", "--ode-solver",
"ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3\n\t."
"\t 22 - Mid-Point, 23 - SDIRK23, 34 - SDIRK34.");
args.AddOption(&zero_start, "-z", "--zero-start", "-no-z",
"--no-zero-start",
"Initial guess of zero or exact solution.");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&gfprint, "-print", "--print","-no-print","--no-print",
"Print results (grid functions) to disk.");
args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
"Enable or disable VisIt visualization.");
args.AddOption(&vis_steps, "-vs", "--visualization-steps",
"Visualize every n-th timestep.");
args.AddOption(&basename, "-k", "--outputfilename",
"Name of the visit dump files");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
if (irOrder < 0)
{
irOrder = std::max(4, 2 * order - 2);
}
if (isnan(alpha_))
{
alpha_ = 0.0;
}
else
{
alpha_ *= M_PI / 180.0;
}
unit_vec_type_ = prob_;
non_linear_ = coef_type > 0;
// 3. Construct a (serial) mesh of the given size on all processors. We
// can handle triangular and quadrilateral surface meshes with the
// same code.
Mesh *mesh = (n > 0) ?
new Mesh(n, n, (Element::Type)el_type, 1) :
new Mesh(mesh_file, 1, 1);
int dim = mesh->Dimension();
if (prob_ > 1) { mesh->Transform(shiftUnitSquare); }
// 4. This step is no longer needed
// 5. Define a parallel mesh by a partitioning of the serial mesh. Once the
// parallel mesh is defined, the serial mesh can be deleted.
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
delete mesh;
// 7. Determine the list of true (i.e. parallel conforming) essential
// boundary dofs. In this example, the boundary conditions are defined
// by marking all the boundary attributes from the mesh as essential
// (Dirichlet) and converting them to a list of true dofs.
Array<int> ess_tdof_list;
Array<int> ess_bdr(0);
if (pmesh->bdr_attributes.Size())
{
ess_bdr.SetSize(pmesh->bdr_attributes.Max());
ess_bdr = 1;
}
// The following is required for mesh refinement
// mesh->EnsureNCMesh();
// 6. Define the ODE solver used for time integration. Several implicit
// methods are available, including singly diagonal implicit Runge-Kutta
// (SDIRK).
ODESolver *ode_solver;
switch (ode_solver_type)
{
// Implicit L-stable methods
case 1: ode_solver = new BackwardEulerSolver; break;
case 2: ode_solver = new SDIRK23Solver(2); break;
case 3: ode_solver = new SDIRK33Solver; break;
// Implicit A-stable methods (not L-stable)
case 22: ode_solver = new ImplicitMidpointSolver; break;
case 23: ode_solver = new SDIRK23Solver; break;
case 34: ode_solver = new SDIRK34Solver; break;
default:
if (mpi.Root())
{
cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
}
delete mesh;
return 3;
}
// 12. Define the parallel finite element spaces. We use:
//
// H(curl) for electric field,
// H(div) for magnetic flux,
// H(div) for thermal flux,
// H(grad)/H1 for electrostatic potential,
// L2 for temperature
// L2 contains discontinuous "cell-center" finite elements, type 2 is
// "positive"
L2_FECollection L2FEC0(0, dim);
L2_FECollection L2FEC(order-1, dim);
// RT contains Raviart-Thomas "face-centered" vector finite elements with
// continuous normal component.
RT_FECollection HDivFEC(order-1, dim);
// H1 contains continuous "node-centered" Lagrange finite elements.
H1_FECollection HGradFEC(order, dim);
ParFiniteElementSpace L2FESpace0(pmesh, &L2FEC0);
ParFiniteElementSpace L2FESpace(pmesh, &L2FEC);
ParFiniteElementSpace HDivFESpace(pmesh, &HDivFEC);
ParFiniteElementSpace HGradFESpace(pmesh, &HGradFEC);
// The terminology is TrueVSize is the unique (non-redundant) number of dofs
// HYPRE_Int glob_size_l2 = L2FESpace.GlobalTrueVSize();
// HYPRE_Int glob_size_rt = HDivFESpace.GlobalTrueVSize();
HYPRE_Int glob_size_h1 = HGradFESpace.GlobalTrueVSize();
if (mpi.Root())
{
cout << "Number of Temperature unknowns: " << glob_size_h1 << endl;
}
// int Vsize_l2 = L2FESpace.GetVSize();
// int Vsize_rt = HDivFESpace.GetVSize();
// int Vsize_h1 = HGradFESpace.GetVSize();
// grid functions E, B, T, F, P, and w which is the Joule heating
ParGridFunction T_gf(&HGradFESpace);
ParGridFunction dT_gf(&HGradFESpace);
ParGridFunction Qs_gf(&HGradFESpace);
ParGridFunction errorT(&L2FESpace0);
T_gf = 1.0;
dT_gf = 1.0;
// 13. Get the boundary conditions, set up the exact solution grid functions
// These VectorCoefficients have an Eval function. Note that e_exact and
// b_exact in this case are exact analytical solutions, taking a 3-vector
// point as input and returning a 3-vector field
FunctionCoefficient TCoef(TFunc);
ConstantCoefficient zeroCoef(0.0);
ConstantCoefficient SpecificHeatCoef(1.0);
// MatrixFunctionCoefficient ConductionCoef(2, ChiFunc);
FunctionCoefficient HeatSourceCoef(QFunc);
VectorFunctionCoefficient UnitBCoef(2, UnitBFunc);
Qs_gf.ProjectCoefficient(HeatSourceCoef);
T_gf.ProjectBdrCoefficient(TCoef, ess_bdr);
T_gf.GridFunction::ComputeElementL2Errors(TCoef, errorT);
// 14. Initialize the Diffusion operator, the GLVis visualization and print
// the initial energies.
ThermalDiffusionTDO oper(HGradFESpace,
zeroCoef, ess_bdr,
chi_perp_,
chi_para_,
prob_,
coef_type,
UnitBCoef,
SpecificHeatCoef, false,
// ConductionCoef, false,
HeatSourceCoef, false);
// This function initializes all the fields to zero or some provided IC
// oper.Init(F);
socketstream vis_T, vis_Q, vis_errT;
char vishost[] = "localhost";
int visport = 19916;
if (visualization)
{
// Make sure all ranks have sent their 'v' solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
vis_T.precision(8);
vis_Q.precision(8);
vis_errT.precision(8);
int Wx = 0, Wy = 0; // window position
int Ww = 350, Wh = 350; // window size
int offx = Ww+10;//, offy = Wh+45; // window offsets
miniapps::VisualizeField(vis_Q, vishost, visport,
Qs_gf, "Heat Soruce", Wx, Wy, Ww, Wh);
Wx += offx;
miniapps::VisualizeField(vis_T, vishost, visport,
T_gf, "Temperature", Wx, Wy, Ww, Wh);
Wx += offx;
miniapps::VisualizeField(vis_errT, vishost, visport,
errorT, "Error in T", Wx, Wy, Ww, Wh);
}
// VisIt visualization
VisItDataCollection visit_dc(basename, pmesh);
if ( visit )
{
visit_dc.RegisterField("T", &T_gf);
visit_dc.RegisterField("Qs", &Qs_gf);
visit_dc.RegisterField("L2 Error T", &errorT);
visit_dc.SetCycle(0);
visit_dc.SetTime(0.0);
visit_dc.Save();
}
ostringstream oss_errs;
oss_errs << "fourier_nl_errs"
<< "_p" << prob_ << "_c" << coef_type
<< "_e" << (int)floor(log10(chi_para_/chi_perp_));
if (n > 0) { oss_errs << "_n" << n; }
oss_errs << "_o" << order << ".dat";
ofstream ofs_errs;
if (myid == 0) { ofs_errs.open(oss_errs.str().c_str()); }
// 15. Perform time-integration (looping over the time iterations, ti, with a
// time-step dt). The object oper is the MagneticDiffusionOperator which
// has a Mult() method and an ImplicitSolve() method which are used by
// the time integrators.
ode_solver->Init(oper);
double t = 0.0;
double dt_courant = 0.0;
{
double h_min, h_max, kappa_min, kappa_max;
pmesh->GetCharacteristics(h_min, h_max, kappa_min, kappa_max);
dt_courant = 1.0 * h_min * h_min / chi_para_;
}
if (dt < 0.0)
{
dt = dt_courant;
}
if ( myid == 0 )
{
cout << "Using time step: " << dt
<< " (Courant " << dt_courant << ")" << endl;
}
int tsize = HGradFESpace.GetTrueVSize();
Vector T0(tsize), T1(tsize), dT(tsize);
T0 = 0.0; T1 = 0.0; dT = 0.0;
T_gf.ParallelProject(T1);
bool last_step = false;
for (int ti = 1; !last_step; ti++)
{
if (t + dt >= t_final - dt/2)
{
if (myid == 0)
{
cout << "Final Time Reached" << endl;
}
last_step = true;
}
// F is the vector of dofs, t is the current time, and dt is the time step
// to advance.
T0 = T1;
ode_solver->Step(T1, t, dt);
T_gf.Distribute(T1);
TCoef.SetTime(t);
T_gf.GridFunction::ComputeElementL2Errors(TCoef, errorT);
double l2_error_T = T_gf.ComputeL2Error(TCoef);
double maxT = T_gf.ComputeMaxError(zeroCoef);
if ( myid == 0 )
{
ofs_errs << t << '\t' << l2_error_T << endl;
cout << t << '\t' << l2_error_T << endl;
}
add(1.0, T1, -1.0, T0, dT);
dT_gf.Distribute(dT);
double maxDiff = dT_gf.ComputeMaxError(zeroCoef);
if ( !last_step )
{
if ( maxT == 0.0 )
{
last_step = (maxDiff < tol) ? true:false;
}
else if ( maxDiff/maxT < tol )
{
last_step = true;
}
if (last_step && myid == 0)
{
cout << "Converged to Steady State" << endl;
}
}
/*
if (debug == 1)
{
oper.Debug(basename,t);
}
*/
if (gfprint)
{
ostringstream T_name, mesh_name;
T_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "T." << setfill('0') << setw(6) << myid;
mesh_name << basename << "_" << setfill('0') << setw(6) << t << "_"
<< "mesh." << setfill('0') << setw(6) << myid;
ofstream mesh_ofs(mesh_name.str().c_str());
mesh_ofs.precision(8);
pmesh->Print(mesh_ofs);
mesh_ofs.close();
ofstream T_ofs(T_name.str().c_str());
T_ofs.precision(8);
T_gf.Save(T_ofs);
T_ofs.close();
}
if (last_step || (ti % vis_steps) == 0)
{
// Make sure all ranks have sent their 'v' solution before initiating
// another set of GLVis connections (one from each rank):
MPI_Barrier(pmesh->GetComm());
if (visualization)
{
int Wx = 0, Wy = 0; // window position
int Ww = 350, Wh = 350; // window size
int offx = Ww+10;//, offy = Wh+45; // window offsets
Wx += offx;
miniapps::VisualizeField(vis_T, vishost, visport,
T_gf, "Temperature", Wx, Wy, Ww, Wh);
Wx += offx;
miniapps::VisualizeField(vis_errT, vishost, visport,
errorT, "Error in T", Wx, Wy, Ww, Wh);
}
if (visit)
{
visit_dc.SetCycle(ti);
visit_dc.SetTime(t);
visit_dc.Save();
}
}
}
if (visualization)
{
vis_T.close();
vis_errT.close();
}
if (myid == 0) { ofs_errs.close(); }
double loc_T_max = T1.Normlinf();
double T_max = -1.0;
MPI_Allreduce(&loc_T_max, &T_max, 1, MPI_DOUBLE, MPI_MAX,
MPI_COMM_WORLD);
double err1 = T_gf.ComputeL2Error(TCoef);
if (myid == 0)
{
cout << "L2 Error of Solution: " << err1 << endl;
cout << "Maximum Temperature: " << T_max << endl;
cout << "| T - T_exact |/|max T| = " << err1 / T_max << endl;
}
// 16. Free the used memory.
delete ode_solver;
delete pmesh;
return 0;
}
void display_banner(ostream & os)
{
os << "___________ .__ " << endl
<< "\\_ _____/___ __ _________|__| ___________ " << endl
<< " | __)/ _ \\| | \\_ __ \\ |/ __ \\_ __ \\" << endl
<< " | | ( <_> ) | /| | \\/ \\ ___/| | \\/" << endl
<< " \\__ | \\____/|____/ |__| |__|\\___ >__| " << endl
<< " \\/ \\/ " << endl
<< flush;
}
+523
View File
@@ -0,0 +1,523 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#include "fourier_vanEs_solver.hpp"
#ifdef MFEM_USE_MPI
using namespace std;
namespace mfem
{
using namespace miniapps;
/*
void
UnitVectorField::Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip)
{
double x[2];
Vector transip(x, 2);
T.Transform(T.GetIntPoint(), transip);
V.SetSize(2);
if ( prob_ % 2 == 1 )
{
if (unit_vec_type_ == 1)
{
double cx = cos(M_PI * x[0]);
double cy = cos(M_PI * x[1]);
double sx = sin(M_PI * x[0]);
double sy = sin(M_PI * x[1]);
V[0] = -sx * cy;
V[1] = sy * cx;
}
else
{
V[0] = cos(M_PI/6.0);
V[1] = sin(M_PI/6.0);
}
}
else
{
V[0] = -a_ * a_ * x[1];
V[1] = b_ * b_ * x[0];
}
double nrm = V.Norml2();
V *= (nrm > 1e-6 * min(a_,b_)) ? (1.0/nrm) : 0.0;
}
*/
void ChiParaCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
if (nonlin_)
{
K *= pow(T_->Eval(T, ip), 2.5);
}
K *= chi_para_;
}
void ChiPerpCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
bbT_->Eval(K, T, ip);
K *= -1.0;
K(0,0) += 1.0;
K(1,1) += 1.0;
if (nonlin_)
{
K *= 1.0 / sqrt(T_->Eval(T, ip));
}
K *= chi_perp_;
}
void dChiCoef::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
double temp = T_->Eval(T, ip);
double perp_factor = 0.5 * chi_perp_ * pow(temp, -1.5);
double para_factor = 2.5 * chi_para_ * pow(temp, 1.5);
bbT_->Eval(K, T, ip);
K *= perp_factor + para_factor;
K(0,0) -= perp_factor;
K(1,1) -= perp_factor;
}
namespace thermal
{
ThermalDiffusionTDO::ThermalDiffusionTDO(
ParFiniteElementSpace &H1_FESpace,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
double chi_perp,
double chi_para,
int prob,
int coef_type,
VectorCoefficient & UnitB,
Coefficient & c, bool td_c,
Coefficient & Q, bool td_Q)
: TimeDependentOperator(H1_FESpace.GetTrueVSize(), 0.0),
init_(false),
nonLinear_(coef_type == 2),
testGradient_(false),
multCount_(0), solveCount_(0),
T_(&H1_FESpace),
TCoef_(&T_),
unitBCoef_(&UnitB),
// ICoef_(2),
bbTCoef_(*unitBCoef_, *unitBCoef_),
chiPerpCoef_(bbTCoef_, TCoef_, chi_perp, coef_type != 0),
chiParaCoef_(bbTCoef_, TCoef_, chi_para, coef_type != 0),
chiCoef_(chiPerpCoef_, chiParaCoef_),
dChiCoef_(bbTCoef_, TCoef_, chi_perp, chi_para),
impOp_(H1_FESpace,
dTdtBdr, false,
bdr_attr,
c, false,
chiCoef_, coef_type != 0,
dChiCoef_, coef_type != 0,
Q, false, coef_type == 2 ),
newton_(H1_FESpace.GetComm())
{
this->init();
}
ThermalDiffusionTDO::~ThermalDiffusionTDO()
{
}
void
ThermalDiffusionTDO::init()
{
cout << "Entering TDO::Init" << endl;
if ( init_ ) { return; }
newton_.SetPrintLevel(2);
newton_.SetRelTol(1e-10);
newton_.SetAbsTol(0.0);
if ( nonLinear_ && testGradient_ )
{
Vector x(impOp_.Height());
Vector dx(impOp_.Height());
T_.Distribute(x);
cout << "GetTime " << this->GetTime() << endl;
impOp_.SetState(T_, this->GetTime(), 0.1);
cout << "init 0" << endl;
newton_.SetOperator(impOp_);
cout << "init 1" << endl;
cout << "init 2" << endl;
x.Randomize(1);
x.Print(cout);
dx.Randomize(2);
dx *= 0.01;
dx.Print(cout);
cout << "init 3" << endl;
double ratio = newton_.CheckGradient(x, dx);
cout << "CheckGradient returns: " << ratio << endl;
}
init_ = true;
cout << "Leaving TDO::Init" << endl;
}
void
ThermalDiffusionTDO::SetTime(const double time)
{
this->TimeDependentOperator::SetTime(time);
newTime_ = true;
}
void
ThermalDiffusionTDO::Mult(const Vector &T, Vector &dT_dt) const
{
MFEM_ABORT("ThermalDiffusionTDO::Mult should not be called");
}
void
ThermalDiffusionTDO::ImplicitSolve(const double dt,
const Vector &T, Vector &dT_dt)
{
dT_dt = 0.0;
T_.Distribute(T);
impOp_.SetState(T_, this->GetTime(), dt);
Solver & solver = impOp_.GetGradientSolver();
if (!nonLinear_)
{
solver.Mult(impOp_.GetRHS(), dT_dt);
}
else
{
newton_.SetOperator(impOp_);
newton_.SetSolver(solver);
newton_.Mult(impOp_.GetRHS(), dT_dt);
}
solveCount_++;
}
ImplicitDiffOp::ImplicitDiffOp(ParFiniteElementSpace & H1_FESpace,
Coefficient & dTdtBdr, bool tdBdr,
Array<int> & bdr_attr,
Coefficient & heatCap, bool tdCp,
ChiCoef & chi, bool tdChi,
dChiCoef & dchi, bool tdDChi,
Coefficient & heatSource, bool tdQ,
bool nonlinear)
: Operator(H1_FESpace.GetTrueVSize()),
first_(true),
tdBdr_(tdBdr),
tdCp_(tdCp),
tdChi_(tdChi),
tdDChi_(tdDChi),
tdQ_(tdQ),
nonLinear_(nonlinear),
newTime_(true),
newTimeStep_(true),
t_(0.0),
dt_(-1.0),
ess_bdr_attr_(bdr_attr),
bdrCoef_(&dTdtBdr),
cpCoef_(&heatCap),
chiCoef_(&chi),
dChiCoef_(&dchi),
QCoef_(&heatSource),
dtChiCoef_(1.0, *chiCoef_),
T0_(&H1_FESpace),
T1_(&H1_FESpace),
dT_(&H1_FESpace),
gradTCoef_(&T0_),
dtGradTCoef_(-1.0, gradTCoef_),
dtdChiGradTCoef_(*dChiCoef_, dtGradTCoef_),
m0cp_(&H1_FESpace),
s0chi_(&H1_FESpace),
a0_(&H1_FESpace),
dTdt_(&H1_FESpace),
Q_(&H1_FESpace),
Qs_(&H1_FESpace),
rhs_(&H1_FESpace),
RHS_(H1_FESpace.GetTrueVSize()),
// RHS0_(0),
AInv_(NULL),
APrecond_(NULL)
{
H1_FESpace.GetEssentialTrueDofs(ess_bdr_attr_, ess_bdr_tdofs_);
m0cp_.AddDomainIntegrator(new MassIntegrator(*cpCoef_));
s0chi_.AddDomainIntegrator(new DiffusionIntegrator(*chiCoef_));
a0_.AddDomainIntegrator(new MassIntegrator(*cpCoef_));
a0_.AddDomainIntegrator(new DiffusionIntegrator(dtChiCoef_));
if (nonLinear_)
{
a0_.AddDomainIntegrator(new MixedScalarWeakDivergenceIntegrator(
dtdChiGradTCoef_));
}
Qs_.AddDomainIntegrator(new DomainLFIntegrator(*QCoef_));
if (!tdQ_) { Qs_.Assemble(); }
}
ImplicitDiffOp::~ImplicitDiffOp()
{
delete AInv_;
delete APrecond_;
}
void ImplicitDiffOp::SetState(ParGridFunction & T, double t, double dt)
{
T0_ = T;
newTime_ = fabs(t - t_) > 0.0;
newTimeStep_= (fabs(1.0-dt/dt_)>1e-6);
t_ = newTime_ ? t : t_;
dt_ = newTimeStep_ ? dt : dt_;
if (tdBdr_ && (newTime_ || newTimeStep_))
{
bdrCoef_->SetTime(t_ + dt_);
}
if (newTimeStep_ || first_)
{
dtChiCoef_.SetAConst(dt_);
dtGradTCoef_.SetAConst(-dt_);
}
if ((tdCp_ && newTime_) || first_)
{
m0cp_.Update();
m0cp_.Assemble();
m0cp_.Finalize();
}
if (!tdChi_ && first_)
{
s0chi_.Assemble();
s0chi_.Finalize();
ofstream ofsS0("s0_const_initial.mat");
s0chi_.SpMat().Print(ofsS0);
a0_.Assemble();
a0_.Finalize();
}
else if (tdChi_ && newTime_ && !nonLinear_)
{
chiCoef_->SetTemp(T0_);
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
ofstream ofsS0("s0_lin_initial.mat");
s0chi_.SpMat().Print(ofsS0);
a0_.Update();
a0_.Assemble(0);
a0_.Finalize(0);
}
if ((tdQ_ && newTime_) || first_)
{
cout << "Assembling Q" << endl;
QCoef_->SetTime(t_ + dt_);
Qs_.Assemble();
Qs_.ParallelAssemble(RHS_);
cout << "Norm of Q: " << Qs_.Norml2() << endl;
}
first_ = false;
newTime_ = false;
newTimeStep_ = false;
}
void ImplicitDiffOp::Mult(const Vector &dT, Vector &Q) const
{
dT_.Distribute(dT);
add(T0_, dt_, dT_, T1_);
if (tdChi_ && nonLinear_)
{
chiCoef_->SetTemp(T1_);
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
}
else
{
cout << "Well this is a surprise..." << endl;
}
m0cp_.Mult(dT_, Q_);
s0chi_.AddMult(T1_, Q_);
Q_.ParallelAssemble(Q);
Q.SetSubVector(ess_bdr_tdofs_, 0.0);
}
Operator & ImplicitDiffOp::GetGradient(const Vector &dT) const
{
if (tdChi_)
{
if (!nonLinear_)
{
chiCoef_->SetTemp(T0_);
}
else
{
dT_.Distribute(dT);
add(T0_, dt_, dT_, T1_);
chiCoef_->SetTemp(T1_);
dChiCoef_->SetTemp(T1_);
gradTCoef_.SetGridFunction(&T1_);
}
s0chi_.Update();
s0chi_.Assemble(0);
s0chi_.Finalize(0);
a0_.Update();
a0_.Assemble(0);
a0_.Finalize(0);
}
if (!nonLinear_)
{
s0chi_.Mult(T0_, rhs_);
rhs_ -= Qs_;
rhs_ *= -1.0;
}
else
{
rhs_ = Qs_;
}
dTdt_.ProjectBdrCoefficient(*bdrCoef_, ess_bdr_attr_);
a0_.FormLinearSystem(ess_bdr_tdofs_, dTdt_, rhs_, A_, SOL_, RHS_);
return A_;
}
Solver & ImplicitDiffOp::GetGradientSolver() const
{
if (!nonLinear_)
{
Operator & A_op = this->GetGradient(T0_); // T0_ will be ignored
HypreParMatrix & A_hyp = dynamic_cast<HypreParMatrix &>(A_op);
if (tdChi_)
{
delete AInv_; AInv_ = NULL;
delete APrecond_; APrecond_ = NULL;
}
if ( AInv_ == NULL )
{
// A_hyp.Print("A.mat");
HyprePCG * AInv_pcg = NULL;
cout << "Building PCG" << endl;
AInv_pcg = new HyprePCG(A_hyp);
AInv_pcg->SetTol(1e-10);
AInv_pcg->SetMaxIter(200);
AInv_pcg->SetPrintLevel(0);
if ( APrecond_ == NULL )
{
cout << "Building AMG" << endl;
APrecond_ = new HypreBoomerAMG(A_hyp);
APrecond_->SetPrintLevel(0);
AInv_pcg->SetPreconditioner(*APrecond_);
}
AInv_ = AInv_pcg;
}
}
else
{
if (AInv_ == NULL)
{
/*
HypreSmoother *J_hypreSmoother = new HypreSmoother;
J_hypreSmoother->SetType(HypreSmoother::l1Jacobi);
J_hypreSmoother->SetPositiveDiagonal(true);
JPrecond_ = J_hypreSmoother;
GMRESSolver * AInv_gmres = NULL;
cout << "Building GMRES" << endl;
AInv_gmres = new GMRESSolver(T0_.ParFESpace()->GetComm());
AInv_gmres->SetRelTol(1e-12);
AInv_gmres->SetAbsTol(0.0);
AInv_gmres->SetMaxIter(20000);
AInv_gmres->SetPrintLevel(2);
AInv_gmres->SetPreconditioner(*JPrecond_);
AInv_ = AInv_gmres;
*/
HypreGMRES * AInv_gmres = NULL;
cout << "Building HypreGMRES" << endl;
AInv_gmres = new HypreGMRES(T0_.ParFESpace()->GetComm());
AInv_gmres->SetTol(1e-12);
AInv_gmres->SetMaxIter(200);
AInv_gmres->SetPrintLevel(2);
if ( APrecond_ == NULL )
{
cout << "Building AMG" << endl;
APrecond_ = new HypreBoomerAMG();
APrecond_->SetPrintLevel(0);
AInv_gmres->SetPreconditioner(*APrecond_);
}
AInv_ = AInv_gmres;
}
}
return *AInv_;
}
} // namespace thermal
void
MatrixInverseCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K.Invert();
}
void
ScaledMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip)
{
M_->Eval(K, T, ip); K *= a_;
}
} // namespace mfem
#endif // MFEM_USE_MPI
+362
View File
@@ -0,0 +1,362 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#ifndef MFEM_FOURIER_NL_SOLVER
#define MFEM_FOURIER_NL_SOLVER
#include "../common/pfem_extras.hpp"
#ifdef MFEM_USE_MPI
#include <memory>
#include <iostream>
#include <fstream>
namespace mfem
{
/*
class UnitVectorField : public VectorCoefficient
{
private:
int prob_;
int unit_vec_type_;
double a_;
double b_;
public:
UnitVectorField(int prob, int unit_vec_type, double a = 0.4, double b = 0.8)
: VectorCoefficient(2), prob_(prob), unit_vec_type_(unit_vec_type),
a_(a), b_(b) {}
void Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip);
};
*/
class ChiParaCoef : public MatrixCoefficient
{
private:
MatrixCoefficient * bbT_;
GridFunctionCoefficient * T_;
double chi_para_;
bool nonlin_;
public:
ChiParaCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_para, bool nonlin = false)
: MatrixCoefficient(2), bbT_(&bbT), T_(&T),
chi_para_(chi_para), nonlin_(nonlin)
{}
void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
class ChiPerpCoef : public MatrixCoefficient
{
private:
MatrixCoefficient * bbT_;
GridFunctionCoefficient * T_;
double chi_perp_;
bool nonlin_;
public:
ChiPerpCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_perp, bool nonlin = false)
: MatrixCoefficient(2), bbT_(&bbT), T_(&T),
chi_perp_(chi_perp), nonlin_(nonlin)
{}
void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
class ChiCoef : public MatrixSumCoefficient
{
private:
ChiPerpCoef * chiPerpCoef_;
ChiParaCoef * chiParaCoef_;
public:
ChiCoef(ChiPerpCoef & chiPerp, ChiParaCoef & chiPara)
: MatrixSumCoefficient(chiPerp, chiPara),
chiPerpCoef_(&chiPerp), chiParaCoef_(&chiPara) {}
void SetTemp(GridFunction & T)
{ chiPerpCoef_->SetTemp(T); chiParaCoef_->SetTemp(T); }
};
class dChiCoef : public MatrixCoefficient
{
private:
MatrixCoefficient * bbT_;
GridFunctionCoefficient * T_;
double chi_perp_;
double chi_para_;
public:
dChiCoef(MatrixCoefficient &bbT, GridFunctionCoefficient &T,
double chi_perp, double chi_para)
: MatrixCoefficient(2), bbT_(&bbT), T_(&T),
chi_perp_(chi_perp), chi_para_(chi_para)
{}
void SetTemp(GridFunction & T) { T_->SetGridFunction(&T); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
};
namespace thermal
{
class ImplicitDiffOp : public Operator
{
public:
ImplicitDiffOp(ParFiniteElementSpace & H1_FESpace,
Coefficient & dTdtBdr, bool tdBdr,
Array<int> & bdr_attr,
Coefficient & heatCap, bool tdCp,
ChiCoef & chi, bool tdChi,
dChiCoef & dchi, bool tdDChi,
Coefficient & heatSource, bool tdQ,
bool nonlinear = false);
~ImplicitDiffOp();
void SetState(ParGridFunction & T, double t, double dt);
void Mult(const Vector &x, Vector &y) const;
Operator & GetGradient(const Vector &x) const;
Solver & GetGradientSolver() const;
const Vector & GetRHS() const { return RHS_; }
private:
bool first_;
bool tdBdr_;
bool tdCp_;
bool tdChi_;
bool tdDChi_;
bool tdQ_;
bool nonLinear_;
bool newTime_;
bool newTimeStep_;
double t_;
double dt_;
Array<int> & ess_bdr_attr_;
Array<int> ess_bdr_tdofs_;
Coefficient * bdrCoef_;
Coefficient * cpCoef_;
ChiCoef * chiCoef_;
dChiCoef * dChiCoef_;
Coefficient * QCoef_;
ScalarMatrixProductCoefficient dtChiCoef_;
mutable ParGridFunction T0_;
mutable ParGridFunction T1_;
mutable ParGridFunction dT_;
mutable GradientGridFunctionCoefficient gradTCoef_;
ScalarVectorProductCoefficient dtGradTCoef_;
MatVecCoefficient dtdChiGradTCoef_;
ParBilinearForm m0cp_;
mutable ParBilinearForm s0chi_;
mutable ParBilinearForm a0_;
mutable HypreParMatrix A_;
mutable ParGridFunction dTdt_;
mutable ParLinearForm Q_;
mutable ParLinearForm Qs_;
mutable ParLinearForm rhs_;
mutable Vector SOL_;
mutable Vector RHS_;
// Vector RHS0_; // Dummy RHS vector which hase length zero
mutable Solver * AInv_;
mutable HypreBoomerAMG * APrecond_;
};
/**
The thermal diffusion equation can be written:
dcT/dt = Div (chi Grad T) + Q_s
where
T is the temperature.
Div is the divergence operator,
grad is the gradient operator,
chi is the thermal conductivity tensor,
c is the heat capacity,
Q_s is the heat source
Class ThermalDiffusionTDO represents the right-hand side of the above
system of ODEs.
f(t, T) = -M_0(c)^{-1}(S_0(chi)T - M_0 Q_s)
where
M_0(c) is an H_1 mass matrix
S_0(chi) is the diffusion operator
The implicit solve method will solve
(M_0(c)+dt S_0(chi))k = -S_0(chi)T + M_0 Q_s
*/
class ThermalDiffusionTDO : public TimeDependentOperator
{
public:
ThermalDiffusionTDO(ParFiniteElementSpace &H1_FES,
Coefficient & dTdtBdr,
Array<int> & bdr_attr,
double chi_perp,
double chi_para,
int prob,
int coef_type,
VectorCoefficient & UnitB,
Coefficient & c, bool td_c,
Coefficient & Q, bool td_Q);
void SetTime(const double time);
/** @brief Perform the action of the operator: @a q = f(@a y, t), where
q solves the algebraic equation F(@a y, q, t) = G(@a y, t) and t is the
current time. */
virtual void Mult(const Vector &y, Vector &q) const;
/** @brief Solve the equation: @a q = f(@a y + @a dt @a q, t), for the
unknown @a q at the current time t.
For general F and G, the equation for @a q becomes:
F(@a y + @a dt @a q, @a q, t) = G(@a y + @a dt @a q, t).
The input vector @a y corresponds to time index (or cycle) n, while the
currently set time, #t, and the result vector @a q correspond to time
index n+1. The time step @a dt corresponds to the time interval between
cycles n and n+1.
This method allows for the abstract implementation of some time
integration methods, including diagonal implicit Runge-Kutta (DIRK)
methods and the backward Euler method in particular.
If not re-implemented, this method simply generates an error. */
virtual void ImplicitSolve(const double dt, const Vector &y, Vector &q);
virtual ~ThermalDiffusionTDO();
private:
void init();
bool init_;
bool newTime_;
bool nonLinear_;
bool testGradient_;
mutable int multCount_;
int solveCount_;
mutable ParGridFunction T_;
GridFunctionCoefficient TCoef_;
VectorCoefficient * unitBCoef_;
// IdentityMatrixCoefficient ICoef_;
OuterProductCoefficient bbTCoef_;
ChiPerpCoef chiPerpCoef_;
ChiParaCoef chiParaCoef_;
ChiCoef chiCoef_;
dChiCoef dChiCoef_;
ImplicitDiffOp impOp_;
NewtonSolver newton_;
};
} // namespace thermal
class InverseCoefficient : public Coefficient
{
public:
InverseCoefficient(Coefficient & c) : c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return 1.0 / c_->Eval(T, ip); }
private:
Coefficient * c_;
};
class MatrixInverseCoefficient :public MatrixCoefficient
{
public:
MatrixInverseCoefficient(MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
MatrixCoefficient * M_;
};
class ScaledCoefficient : public Coefficient
{
public:
ScaledCoefficient(double a, Coefficient & c) : a_(a), c_(&c) {}
void SetTime(double t) { time = t; c_->SetTime(t); }
double Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{ return a_ * c_->Eval(T, ip); }
private:
double a_;
Coefficient * c_;
};
class ScaledMatrixCoefficient :public MatrixCoefficient
{
public:
ScaledMatrixCoefficient(double a, MatrixCoefficient & M)
: MatrixCoefficient(M.GetWidth()), a_(a), M_(&M) {}
void SetTime(double t) { time = t; M_->SetTime(t); }
void Eval(DenseMatrix &K, ElementTransformation &T,
const IntegrationPoint &ip);
private:
double a_;
MatrixCoefficient * M_;
};
} // namespace mfem
#endif // MFEM_USE_MPI
#endif // MFEM_FOURIER_NL_SOLVER
+93
View File
@@ -0,0 +1,93 @@
# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
# See file COPYRIGHT for details.
#
# This file is part of the MFEM library. For more information and source code
# availability see http://mfem.org.
#
# MFEM is free software; you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License (as published by the Free
# Software Foundation) version 2.1 dated February 1999.
# Use the MFEM build directory
MFEM_DIR ?= ../..
MFEM_BUILD_DIR ?= ../..
SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/thermal/,)
CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
# Use the MFEM install directory
# MFEM_INSTALL_DIR = ../../mfem
# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
MFEM_LIB_FILE = mfem_is_not_built
-include $(CONFIG_MK)
SEQ_MINIAPPS =
PAR_MINIAPPS = fourier fourier_nl fourier_vanEs fourier_hybrid \
fourier_flux fourier_nl_flux \
fourier_refine fourier_flux_refine ex1p_nl
ifeq ($(MFEM_USE_MPI),NO)
MINIAPPS = $(SEQ_MINIAPPS)
else
MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
endif
.SUFFIXES:
.SUFFIXES: .o .cpp .mk
.PHONY: all clean clean-build clean-exec
.PRECIOUS: %.o
COMMON_O=../common/pfem_extras.o
# Remove built-in rules
%: %.cpp
%.o: %.cpp
all: $(MINIAPPS)
# Rules for building the miniapps
%: $(SRC)%.cpp %_solver.o $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $@_solver.o $(COMMON_O) $(MFEM_LIBS)
fourier_refine: fourier_refine.cpp fourier_solver.o $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ fourier_solver.o $(COMMON_O) $(MFEM_LIBS)
fourier_flux_refine: fourier_flux_refine.cpp fourier_flux_solver.o $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ fourier_flux_solver.o $(COMMON_O) $(MFEM_LIBS)
curve_mesh: curve_mesh.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(COMMON_O) $(MFEM_LIBS)
ncd2mesh: ncd2mesh.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(COMMON_O) $(MFEM_LIBS)
ex1p_nl: ex1p_nl.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(COMMON_O) $(MFEM_LIBS)
# Rules for compiling miniapp dependencies
$(COMMON_O) $(addsuffix _solver.o,$(MINIAPPS)): \
%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
MFEM_TESTS = MINIAPPS
include $(MFEM_TEST_MK)
# Testing: Specific execution options
RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
fourier-test-par: fourier
@$(call mfem-test,$<, $(RUN_MPI), Thermal miniapp,\
)
# Testing: "test" target and mfem-test* variables are defined in config/test.mk
# Generate an error message if the MFEM library is not built and exit
$(MFEM_LIB_FILE):
$(error The MFEM library is not built)
clean: clean-build clean-exec
clean-build:
rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
rm -rf *.dSYM *.TVD.*breakpoints
clean-exec:
@rm -rf Fourier_*