Skip to content

Commit

Permalink
Moved global variables in the operator.
Browse files Browse the repository at this point in the history
  • Loading branch information
@vladotomov
vladotomov committed Nov 8, 2018
1 parent 25d2f73 commit 3798a94
Show file tree
Hide file tree
Showing 8 changed files with 56 additions and 65 deletions.
3 changes: 2 additions & 1 deletion laghos.cpp
View file
Expand Up @@ -418,7 +418,8 @@ int main(int argc, char *argv[])

LagrangianHydroOperator oper(S.Size(), H1FESpace, L2FESpace,
ess_tdofs, rho, source, cfl, mat_gf_coeff,
visc, p_assembly, cg_tol, cg_max_iter);
visc, p_assembly, cg_tol, cg_max_iter,
H1FEC.GetBasisType());

socketstream vis_rho, vis_v, vis_e;
char vishost[] = "localhost";
Expand Down
3 changes: 0 additions & 3 deletions laghos_assembly.cpp
View file
Expand Up @@ -24,9 +24,6 @@ namespace mfem
namespace hydrodynamics
{

const Tensors1D *tensors1D = NULL;
const FastEvaluator *evaluator = NULL;

Tensors1D::Tensors1D(int H1order, int L2order, int nqp1D, bool bernstein_v)
: HQshape1D(H1order + 1, nqp1D),
HQgrad1D(H1order + 1, nqp1D),
Expand Down
29 changes: 17 additions & 12 deletions laghos_assembly.hpp
View file
Expand Up @@ -70,23 +70,22 @@ struct Tensors1D

Tensors1D(int H1order, int L2order, int nqp1D, bool bernstein_v);
};
extern const Tensors1D *tensors1D;

class FastEvaluator
{
const int dim;
FiniteElementSpace &H1FESpace;
Tensors1D *tensors1D;

public:
FastEvaluator(FiniteElementSpace &h1fes)
: dim(h1fes.GetMesh()->Dimension()), H1FESpace(h1fes) { }
FastEvaluator(FiniteElementSpace &h1fes, Tensors1D *t1D)
: dim(h1fes.GetMesh()->Dimension()), H1FESpace(h1fes), tensors1D(t1D) { }

void GetL2Values(const Vector &vecL2, Vector &vecQP) const;
// The input vec is an H1 function with dim components, over a zone.
// The output is J_ij = d(vec_i) / d(x_j) with ij = 1 .. dim.
void GetVectorGrad(const DenseMatrix &vec, DenseTensor &J) const;
};
extern const FastEvaluator *evaluator;

// This class is used only for visualization. It assembles (rho, phi) in each
// zone, which is used by LagrangianHydroOperator::ComputeDensity to do an L2
Expand Down Expand Up @@ -129,6 +128,7 @@ class ForcePAOperator : public Operator

QuadratureData *quad_data;
FiniteElementSpace &H1FESpace, &L2FESpace;
Tensors1D *tensors1D;

// Force matrix action on quadrilateral elements in 2D.
void MultQuad(const Vector &vecL2, Vector &vecH1) const;
Expand All @@ -142,9 +142,11 @@ class ForcePAOperator : public Operator

public:
ForcePAOperator(QuadratureData *quad_data_,
FiniteElementSpace &h1fes, FiniteElementSpace &l2fes)
FiniteElementSpace &h1fes, FiniteElementSpace &l2fes,
Tensors1D *t1D)
: dim(h1fes.GetMesh()->Dimension()), nzones(h1fes.GetMesh()->GetNE()),
quad_data(quad_data_), H1FESpace(h1fes), L2FESpace(l2fes) { }
quad_data(quad_data_), H1FESpace(h1fes), L2FESpace(l2fes),
tensors1D(t1D) { }

virtual void Mult(const Vector &vecL2, Vector &vecH1) const;
virtual void MultTranspose(const Vector &vecH1, Vector &vecL2) const;
Expand All @@ -160,18 +162,19 @@ class MassPAOperator : public Operator

QuadratureData *quad_data;
FiniteElementSpace &FESpace;
Tensors1D *tensors1D;

// Mass matrix action on quadrilateral elements in 2D.
void MultQuad(const Vector &x, Vector &y) const;
// Mass matrix action on hexahedral elements in 3D.
void MultHex(const Vector &x, Vector &y) const;

public:
MassPAOperator(QuadratureData *quad_data_, FiniteElementSpace &fes)
MassPAOperator(QuadratureData *quad_data_, FiniteElementSpace &fes,
Tensors1D *t1D)
: Operator(fes.GetVSize()),
dim(fes.GetMesh()->Dimension()), nzones(fes.GetMesh()->GetNE()),
quad_data(quad_data_), FESpace(fes)
{ }
quad_data(quad_data_), FESpace(fes), tensors1D(t1D) { }

// Mass matrix action.
virtual void Mult(const Vector &x, Vector &y) const;
Expand Down Expand Up @@ -223,18 +226,20 @@ class LocalMassPAOperator : public Operator
int zone_id;

QuadratureData *quad_data;
Tensors1D *tensors1D;

// Mass matrix action on a quadrilateral element in 2D.
void MultQuad(const Vector &x, Vector &y) const;
// Mass matrix action on a hexahedral element in 3D.
void MultHex(const Vector &x, Vector &y) const;

public:
LocalMassPAOperator(QuadratureData *quad_data_, FiniteElementSpace &fes)
LocalMassPAOperator(QuadratureData *quad_data_, FiniteElementSpace &fes,
Tensors1D *t1D)
: Operator(fes.GetFE(0)->GetDof()),
dim(fes.GetMesh()->Dimension()), zone_id(0),
quad_data(quad_data_)
{ }
quad_data(quad_data_), tensors1D(t1D) { }

void SetZoneId(int zid) { zone_id = zid; }

virtual void Mult(const Vector &x, Vector &y) const;
Expand Down
34 changes: 13 additions & 21 deletions laghos_solver.cpp
View file
Expand Up @@ -83,7 +83,8 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size,
int source_type_, double cfl_,
Coefficient *material_,
bool visc, bool pa,
double cgt, int cgiter)
double cgt, int cgiter,
int h1_basis_type)
: TimeDependentOperator(size),
H1FESpace(h1_fes), L2FESpace(l2_fes),
ess_tdofs(essential_tdofs),
Expand All @@ -100,11 +101,16 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size,
3*h1_fes.GetOrder(0) + l2_fes.GetOrder(0) - 1)),
quad_data(dim, nzones, integ_rule.GetNPoints()),
quad_data_is_current(false), forcemat_is_assembled(false),
Force(&l2_fes, &h1_fes), ForcePA(&quad_data, h1_fes, l2_fes),
VMassPA(&quad_data, H1FESpace), VMassPA_prec(H1FESpace),
locEMassPA(&quad_data, l2_fes),
tensors1D(H1FESpace.GetFE(0)->GetOrder(), L2FESpace.GetFE(0)->GetOrder(),
int(floor(0.7 + pow(integ_rule.GetNPoints(), 1.0 / dim))),
h1_basis_type == BasisType::Positive),
evaluator(H1FESpace, &tensors1D),
Force(&l2_fes, &h1_fes), ForcePA(&quad_data, h1_fes, l2_fes, &tensors1D),
VMassPA(&quad_data, H1FESpace, &tensors1D), VMassPA_prec(H1FESpace),
locEMassPA(&quad_data, l2_fes, &tensors1D),
locCG(), timer()
{

GridFunctionCoefficient rho_coeff(&rho0);

// Standard local assembly and inversion for energy mass matrices.
Expand Down Expand Up @@ -170,15 +176,6 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size,

if (p_assembly)
{
// Compute the global 1D reference tensors.
const H1_FECollection *h1_fec =
dynamic_cast<const H1_FECollection *>(H1FESpace.FEColl());
tensors1D = new Tensors1D(H1FESpace.GetFE(0)->GetOrder(),
L2FESpace.GetFE(0)->GetOrder(),
int(floor(0.7 + pow(nqp, 1.0 / dim))),
h1_fec->GetBasisType() == BasisType::Positive);
evaluator = new FastEvaluator(H1FESpace);

// Setup the preconditioner of the velocity mass operator.
Vector d;
(dim == 2) ? VMassPA.ComputeDiagonal2D(d) : VMassPA.ComputeDiagonal3D(d);
Expand Down Expand Up @@ -479,11 +476,6 @@ void LagrangianHydroOperator::PrintTimingData(bool IamRoot, int steps) const
}
}

LagrangianHydroOperator::~LagrangianHydroOperator()
{
delete tensors1D;
}

// Smooth transition between 0 and 1 for x in [-eps, eps].
inline double smooth_step_01(double x, double eps)
{
Expand Down Expand Up @@ -548,12 +540,12 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const
// Energy values at quadrature point.
L2FESpace.GetElementDofs(z_id, L2dofs);
e.GetSubVector(L2dofs, e_loc);
evaluator->GetL2Values(e_loc, e_vals);
evaluator.GetL2Values(e_loc, e_vals);

// All reference->physical Jacobians at the quadrature points.
H1FESpace.GetElementVDofs(z_id, H1dofs);
x.GetSubVector(H1dofs, vector_vals);
evaluator->GetVectorGrad(vecvalMat, Jpr_b[z]);
evaluator.GetVectorGrad(vecvalMat, Jpr_b[z]);
}
else { e.GetValues(z_id, integ_rule, e_vals); }
for (int q = 0; q < nqp; q++)
Expand Down Expand Up @@ -585,7 +577,7 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const
// All reference->physical Jacobians at the quadrature points.
H1FESpace.GetElementVDofs(z_id, H1dofs);
v.GetSubVector(H1dofs, vector_vals);
evaluator->GetVectorGrad(vecvalMat, grad_v_ref);
evaluator.GetVectorGrad(vecvalMat, grad_v_ref);
}
for (int q = 0; q < nqp; q++)
{
Expand Down
8 changes: 5 additions & 3 deletions laghos_solver.hpp
View file
Expand Up @@ -84,6 +84,10 @@ class LagrangianHydroOperator : public TimeDependentOperator
mutable QuadratureData quad_data;
mutable bool quad_data_is_current, forcemat_is_assembled;

// Structures used to perform partial assembly.
Tensors1D tensors1D;
FastEvaluator evaluator;

// Force matrix that combines the kinematic and thermodynamic spaces. It is
// assembled in each time step and then it is used to compute the final
// right-hand sides for momentum and specific internal energy.
Expand Down Expand Up @@ -123,7 +127,7 @@ class LagrangianHydroOperator : public TimeDependentOperator
Array<int> &essential_tdofs, ParGridFunction &rho0,
int source_type_, double cfl_,
Coefficient *material_, bool visc, bool pa,
double cgt, int cgiter);
double cgt, int cgiter, int h1_basis_type);

// Solve for dx_dt, dv_dt and de_dt.
virtual void Mult(const Vector &S, Vector &dS_dt) const;
Expand All @@ -147,8 +151,6 @@ class LagrangianHydroOperator : public TimeDependentOperator
void PrintTimingData(bool IamRoot, int steps) const;

int GetH1VSize() const { return H1FESpace.GetVSize(); }

~LagrangianHydroOperator();
};

class TaylorCoefficient : public Coefficient
Expand Down
3 changes: 2 additions & 1 deletion serial/laghos.cpp
View file
Expand Up @@ -271,7 +271,8 @@ int main(int argc, char *argv[])

LagrangianHydroOperator oper(S.Size(), H1FESpace, L2FESpace,
ess_vdofs, rho, source, cfl, mat_gf_coeff,
visc, p_assembly, cg_tol, cg_max_iter);
visc, p_assembly, cg_tol, cg_max_iter,
H1FEC.GetBasisType());

socketstream vis_rho, vis_v, vis_e;
char vishost[] = "localhost";
Expand Down
33 changes: 12 additions & 21 deletions serial/laghos_solver.cpp
View file
Expand Up @@ -68,7 +68,8 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size,
int source_type_, double cfl_,
Coefficient *material_,
bool visc, bool pa,
double cgt, int cgiter)
double cgt, int cgiter,
int h1_basis_type)
: TimeDependentOperator(size),
H1FESpace(h1_fes), L2FESpace(l2_fes),
ess_tdofs(essential_tdofs),
Expand All @@ -85,9 +86,13 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size,
3*h1_fes.GetOrder(0) + l2_fes.GetOrder(0) - 1)),
quad_data(dim, nzones, integ_rule.GetNPoints()),
quad_data_is_current(false), forcemat_is_assembled(false),
Force(&l2_fes, &h1_fes), ForcePA(&quad_data, h1_fes, l2_fes),
VMassPA(&quad_data, H1FESpace), VMassPA_prec(H1FESpace),
locEMassPA(&quad_data, l2_fes),
tensors1D(H1FESpace.GetFE(0)->GetOrder(), L2FESpace.GetFE(0)->GetOrder(),
int(floor(0.7 + pow(integ_rule.GetNPoints(), 1.0 / dim))),
h1_basis_type == BasisType::Positive),
evaluator(H1FESpace, &tensors1D),
Force(&l2_fes, &h1_fes), ForcePA(&quad_data, h1_fes, l2_fes, &tensors1D),
VMassPA(&quad_data, H1FESpace, &tensors1D), VMassPA_prec(H1FESpace),
locEMassPA(&quad_data, l2_fes, &tensors1D),
locCG(), timer()
{
GridFunctionCoefficient rho_coeff(&rho0);
Expand Down Expand Up @@ -159,15 +164,6 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size,

if (p_assembly)
{
// Compute the global 1D reference tensors.
const H1_FECollection *h1_fec =
dynamic_cast<const H1_FECollection *>(H1FESpace.FEColl());
tensors1D = new Tensors1D(H1FESpace.GetFE(0)->GetOrder(),
L2FESpace.GetFE(0)->GetOrder(),
int(floor(0.7 + pow(nqp, 1.0 / dim))),
h1_fec->GetBasisType() == BasisType::Positive);
evaluator = new FastEvaluator(H1FESpace);

// Setup the preconditioner of the velocity mass operator.
Vector d;
(dim == 2) ? VMassPA.ComputeDiagonal2D(d) : VMassPA.ComputeDiagonal3D(d);
Expand Down Expand Up @@ -441,11 +437,6 @@ void LagrangianHydroOperator::PrintTimingData(int steps) const
<< 1e-6 * steps * (H1size + L2size) / runtime[4] << endl;
}

LagrangianHydroOperator::~LagrangianHydroOperator()
{
delete tensors1D;
}

// Smooth transition between 0 and 1 for x in [-eps, eps].
inline double smooth_step_01(double x, double eps)
{
Expand Down Expand Up @@ -510,12 +501,12 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const
// Energy values at quadrature point.
L2FESpace.GetElementDofs(z_id, L2dofs);
e.GetSubVector(L2dofs, e_loc);
evaluator->GetL2Values(e_loc, e_vals);
evaluator.GetL2Values(e_loc, e_vals);

// All reference->physical Jacobians at the quadrature points.
H1FESpace.GetElementVDofs(z_id, H1dofs);
x.GetSubVector(H1dofs, vector_vals);
evaluator->GetVectorGrad(vecvalMat, Jpr_b[z]);
evaluator.GetVectorGrad(vecvalMat, Jpr_b[z]);
}
else { e.GetValues(z_id, integ_rule, e_vals); }
for (int q = 0; q < nqp; q++)
Expand Down Expand Up @@ -547,7 +538,7 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const
// All reference->physical Jacobians at the quadrature points.
H1FESpace.GetElementVDofs(z_id, H1dofs);
v.GetSubVector(H1dofs, vector_vals);
evaluator->GetVectorGrad(vecvalMat, grad_v_ref);
evaluator.GetVectorGrad(vecvalMat, grad_v_ref);
}
for (int q = 0; q < nqp; q++)
{
Expand Down
8 changes: 5 additions & 3 deletions serial/laghos_solver.hpp
View file
Expand Up @@ -82,6 +82,10 @@ class LagrangianHydroOperator : public TimeDependentOperator
mutable QuadratureData quad_data;
mutable bool quad_data_is_current, forcemat_is_assembled;

// Structures used to perform partial assembly.
Tensors1D tensors1D;
FastEvaluator evaluator;

// Force matrix that combines the kinematic and thermodynamic spaces. It is
// assembled in each time step and then it is used to compute the final
// right-hand sides for momentum and specific internal energy.
Expand Down Expand Up @@ -121,7 +125,7 @@ class LagrangianHydroOperator : public TimeDependentOperator
Array<int> &essential_tdofs, GridFunction &rho0,
int source_type_, double cfl_,
Coefficient *material_, bool visc, bool pa,
double cgt, int cgiter);
double cgt, int cgiter, int h1_basis_type);

// Solve for dx_dt, dv_dt and de_dt.
virtual void Mult(const Vector &S, Vector &dS_dt) const;
Expand All @@ -145,8 +149,6 @@ class LagrangianHydroOperator : public TimeDependentOperator
void PrintTimingData(int steps) const;

int GetH1VSize() const { return H1FESpace.GetVSize(); }

~LagrangianHydroOperator();
};

class TaylorCoefficient : public Coefficient
Expand Down

0 comments on commit 3798a94

Please sign in to comment.