From 32b4f4baa4043596e5c35e7ed96a9caca6a8bdfb Mon Sep 17 00:00:00 2001 From: Arturo Vargas Date: Mon, 20 Dec 2021 15:11:13 -0800 Subject: [PATCH 1/7] sketch of MF discrete upwind --- remhos.cpp | 17 +++-- remhos_lo.cpp | 189 ++++++++++++++++++++++++++++++++++++++++++++++++-- remhos_lo.hpp | 47 +++++++++++++ 3 files changed, 242 insertions(+), 11 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 055878d..c9db842 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -390,16 +390,22 @@ int main(int argc, char *argv[]) M_HO.AddDomainIntegrator(new MassIntegrator); ParBilinearForm k(&pfes); + ParBilinearForm k_du(&pfes); ParBilinearForm K_HO(&pfes); + ConvectionIntegrator *conv_int = nullptr; if (exec_mode == 0) { + conv_int = new ConvectionIntegrator(velocity, -1.0); k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0)); + k_du.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0)); K_HO.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0)); } else if (exec_mode == 1) { - k.AddDomainIntegrator(new ConvectionIntegrator(v_coef)); - K_HO.AddDomainIntegrator(new ConvectionIntegrator(v_coef)); + conv_int = new ConvectionIntegrator(v_coef); + k.AddDomainIntegrator(conv_int); + k_du.AddDomainIntegrator(conv_int); + K_HO.AddDomainIntegrator(conv_int); } if (ho_type == HOSolverType::CG || @@ -430,11 +436,14 @@ int main(int argc, char *argv[]) K_HO.SetAssemblyLevel(AssemblyLevel::PARTIAL); } + k_du.SetAssemblyLevel(AssemblyLevel::PARTIAL); + if (next_gen_full) { K_HO.SetAssemblyLevel(AssemblyLevel::FULL); } + k_du.Assemble(); M_HO.Assemble(); K_HO.Assemble(0); @@ -613,8 +622,8 @@ int main(int argc, char *argv[]) if (lo_type == LOSolverType::DiscrUpwind) { lo_smap = SparseMatrix_Build_smap(k.SpMat()); - lo_solver = new DiscreteUpwind(pfes, k.SpMat(), lo_smap, - lumpedM, asmbl, time_dep); + lo_solver = new PADiscreteUpwind(pfes, k_du, conv_int, k.SpMat(), lo_smap, + lumpedM, asmbl, time_dep); } else if (lo_type == LOSolverType::DiscrUpwindPrec) { diff --git a/remhos_lo.cpp b/remhos_lo.cpp index ace6bc6..395e5af 100644 --- a/remhos_lo.cpp +++ b/remhos_lo.cpp @@ -22,6 +22,13 @@ using namespace std; namespace mfem { +const DofToQuad *get_maps(ParFiniteElementSpace &pfes, Assembly &asmbly) +{ + const FiniteElement *el_trace = + pfes.GetTraceElement(0, pfes.GetParMesh()->GetFaceBaseGeometry(0)); + return &el_trace->GetDofToQuad(*asmbly.lom.irF, DofToQuad::TENSOR); +} + DiscreteUpwind::DiscreteUpwind(ParFiniteElementSpace &space, const SparseMatrix &adv, const Array &adv_smap, const Vector &Mlump, @@ -34,6 +41,22 @@ DiscreteUpwind::DiscreteUpwind(ParFiniteElementSpace &space, ComputeDiscreteUpwindMatrix(); } +PADiscreteUpwind::PADiscreteUpwind(ParFiniteElementSpace &space, + ParBilinearForm &Kbf, ConvectionIntegrator *Conv_, + const SparseMatrix &adv, + const Array &adv_smap, const Vector &Mlump, + Assembly &asmbly, bool updateD) + : LOSolver(space), + k_pbilinear(Kbf), Conv(Conv_), K(adv), D(), K_smap(adv_smap), M_lumped(Mlump), + assembly(asmbly), update_D(updateD), + quad1D(get_maps(pfes, assembly)->nqpt), + dofs1D(get_maps(pfes, assembly)->ndof), + face_dofs((pfes.GetMesh()->Dimension() ==2) ? quad1D : quad1D * quad1D) +{ + D = K; + ComputeDiscreteUpwindMatrix(); +} + void DiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const { const int ndof = pfes.GetFE(0)->GetDof(); @@ -67,6 +90,139 @@ void DiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const for (int i = 0; i < s; i++) { du(i) /= M_lumped(i); } } +void PADiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const +{ + //mfem_error("Here! \n"); + const int ndof = pfes.GetFE(0)->GetDof(); + const int NE = pfes.GetMesh()->GetNE(); + Vector alpha(ndof); alpha = 0.0; + + // Recompute D due to mesh changes (K changes) in remap mode. + if (update_D) { ComputeDiscreteUpwindMatrix(); } + + // Discretization and monotonicity terms. + D.Mult(u, du); + //K.Mult(u, du); + + Vector y(du.Size()); y = 0.0; + + //Clearly incorrect... + //k_pbilinear.Mult(u, y); + + Vector ConvMats(ndof * ndof * NE); + Vector AlgDiffMats(ndof * ndof * NE); AlgDiffMats = 0.0; + + //Given in row major format + Conv->AssembleEA(pfes, ConvMats, false); + + ComputeAlgebraicDiffusion(ConvMats, AlgDiffMats); + + AddBlkMult(ConvMats, u, y); + + AddBlkMult(AlgDiffMats, u, y); + + + y-=du; + double error = y.Norml2(); + if (error > 1e-12) + { + std::cout<<"Error too high "<GetDof(); + const int NE = pfes.GetMesh()->GetNE(); + auto conv_blk = Reshape(ConvMats.Read(), ndof, ndof, NE); + auto alg_blk = Reshape(AlgDiff.ReadWrite(), ndof, ndof, NE); + + //Matrices are assumed to be in row major format + for (int e=0; eGetDof(); + const int NE = pfes.GetMesh()->GetNE(); + + auto X = Reshape(x.Read(), ndof, NE); + auto Y = Reshape(y.Write(), ndof, NE); + auto Me = Reshape(Mat.Read(), ndof, ndof, NE); + + //Takes row major format + for (int e=0; eGetFaceBaseGeometry(0)); - return &el_trace->GetDofToQuad(*asmbly.lom.irF, DofToQuad::TENSOR); -} - //==== //Residual Distribution // diff --git a/remhos_lo.hpp b/remhos_lo.hpp index f0133c1..cf042cf 100644 --- a/remhos_lo.hpp +++ b/remhos_lo.hpp @@ -58,6 +58,53 @@ class DiscreteUpwind : public LOSolver virtual void CalcLOSolution(const Vector &u, Vector &du) const; }; +class PADiscreteUpwind : public LOSolver +{ +protected: + ParBilinearForm &k_pbilinear; + ConvectionIntegrator *Conv; + const SparseMatrix &K; + mutable SparseMatrix D; + const Array &K_smap; + const Vector &M_lumped; + Assembly &assembly; + const bool update_D; + + void ComputeDiscreteUpwindMatrix() const; + + // Data at quadrature points //for face terms + const int quad1D, dofs1D, face_dofs; + mutable Array D_int, D_bdry; + mutable Array IntVelocity, BdryVelocity; + +public: + PADiscreteUpwind(ParFiniteElementSpace &space, + ParBilinearForm &Kbf, ConvectionIntegrator *Conv_, + const SparseMatrix &adv, + const Array &adv_smap, const Vector &Mlump, + Assembly &asmbly, bool updateD); + + void SampleVelocity(FaceType type) const; + + void SetupPA(FaceType type) const; + + void SetupPA2D(FaceType) const; + + void SetupPA3D(FaceType) const; + + void ApplyFaceTerms(const Vector &x, Vector &y, FaceType type) const; + + void ApplyFaceTerms2D(const Vector &x, Vector &y, FaceType type) const; + + void ApplyFaceTerms3D(const Vector &x, Vector &y, FaceType type) const; + + void ComputeAlgebraicDiffusion(Vector &ConvMats, Vector &AlgDiff) const; + + void AddBlkMult(const Vector &Mat, const Vector &x, Vector &y) const; + + virtual void CalcLOSolution(const Vector &u, Vector &du) const; +}; + class ResidualDistribution : public LOSolver { protected: From 53757e5d2cbd7622fd286cf2c3b11592d9508a73 Mon Sep 17 00:00:00 2001 From: Arturo Vargas Date: Mon, 20 Dec 2021 15:57:49 -0800 Subject: [PATCH 2/7] add trace class for faceterms --- remhos_lo.cpp | 75 ++++++++++++++++++++++++++------------------------- remhos_lo.hpp | 66 ++++++++++++++++++++++----------------------- 2 files changed, 72 insertions(+), 69 deletions(-) diff --git a/remhos_lo.cpp b/remhos_lo.cpp index 395e5af..98619e9 100644 --- a/remhos_lo.cpp +++ b/remhos_lo.cpp @@ -429,21 +429,24 @@ PAResidualDistribution::PAResidualDistribution(ParFiniteElementSpace &space, const Vector &Mlump, bool subcell, bool timedep) : ResidualDistribution(space, Kbf, asmbly, Mlump, subcell, timedep), - quad1D(get_maps(pfes, assembly)->nqpt), - dofs1D(get_maps(pfes, assembly)->ndof), - face_dofs((pfes.GetMesh()->Dimension() ==2) ? quad1D : quad1D * quad1D) + PADGTraceLOSolver(space, get_maps(pfes, trace_assembly)->nqpt, + get_maps(pfes, trace_assembly)->ndof, + (pfes.GetMesh()->Dimension() ==2) ? + get_maps(pfes, trace_assembly)->nqpt : + get_maps(pfes, trace_assembly)->nqpt * get_maps(pfes, trace_assembly)->nqpt, + assembly) { } // Taken from DGTraceIntegrator::SetupPA L:145 -void PAResidualDistribution::SampleVelocity(FaceType type) const +void PADGTraceLOSolver::SampleVelocity(FaceType type) const { - const int nf = pfes.GetNFbyType(type); + const int nf = trace_pfes.GetNFbyType(type); if (nf == 0) { return; } - const IntegrationRule *ir = assembly.lom.irF; + const IntegrationRule *ir = trace_assembly.lom.irF; - Mesh *mesh = pfes.GetMesh(); + Mesh *mesh = trace_pfes.GetMesh(); const int dim = mesh->Dimension(); const int nq = ir->GetNPoints(); @@ -464,7 +467,7 @@ void PAResidualDistribution::SampleVelocity(FaceType type) const Vector Vq(dim); int f_idx = 0; - for (int f = 0; f < pfes.GetNF(); ++f) + for (int f = 0; f < trace_pfes.GetNF(); ++f) { int e1, e2; int inf1, inf2; @@ -484,7 +487,7 @@ void PAResidualDistribution::SampleVelocity(FaceType type) const int iq = ToLexOrdering(dim, my_face_id, quad1D, q); T.SetAllIntPoints(&ir->IntPoint(q)); const IntegrationPoint &eip1 = T.GetElement1IntPoint(); - assembly.lom.coef->Eval(Vq, *T.Elem1, eip1); + trace_assembly.lom.coef->Eval(Vq, *T.Elem1, eip1); for (int i = 0; i < dim; ++i) { C(i,iq,f_idx) = Vq(i); @@ -495,9 +498,9 @@ void PAResidualDistribution::SampleVelocity(FaceType type) const } } -void PAResidualDistribution::SetupPA(FaceType type) const +void PADGTraceLOSolver::SetupPA(FaceType type) const { - const FiniteElementSpace *fes = assembly.GetFes(); + const FiniteElementSpace *fes = trace_assembly.GetFes(); int nf = fes->GetNFbyType(type); if (nf == 0) {return;} @@ -509,13 +512,13 @@ void PAResidualDistribution::SetupPA(FaceType type) const if (dim == 3) { return SetupPA3D(type); } } -void PAResidualDistribution::SetupPA2D(FaceType type) const +void PADGTraceLOSolver::SetupPA2D(FaceType type) const { - const int nf = pfes.GetNFbyType(type); - Mesh *mesh = pfes.GetMesh(); + const int nf = trace_pfes.GetNFbyType(type); + Mesh *mesh = trace_pfes.GetMesh(); const int dim = mesh->Dimension(); - const IntegrationRule *ir = assembly.lom.irF; + const IntegrationRule *ir = trace_assembly.lom.irF; const FaceGeometricFactors *geom = mesh->GetFaceGeometricFactors(*ir, @@ -532,7 +535,7 @@ void PAResidualDistribution::SetupPA2D(FaceType type) const auto vel = mfem::Reshape(vel_ptr, dim, ir->GetNPoints(), nf); - const int execMode = (int) assembly.GetExecMode(); + const int execMode = (int) trace_assembly.GetExecMode(); if (type == FaceType::Interior) { @@ -600,13 +603,13 @@ void PAResidualDistribution::SetupPA2D(FaceType type) const }//boundary } -void PAResidualDistribution::SetupPA3D(FaceType type) const +void PADGTraceLOSolver::SetupPA3D(FaceType type) const { - const int nf = pfes.GetNFbyType(type); - Mesh *mesh = pfes.GetMesh(); + const int nf = trace_pfes.GetNFbyType(type); + Mesh *mesh = trace_pfes.GetMesh(); const int dim = mesh->Dimension(); - const IntegrationRule *ir = assembly.lom.irF; + const IntegrationRule *ir = trace_assembly.lom.irF; const FaceGeometricFactors *geom = mesh->GetFaceGeometricFactors(*ir, @@ -623,7 +626,7 @@ void PAResidualDistribution::SetupPA3D(FaceType type) const auto vel = mfem::Reshape(vel_ptr, dim, quad1D, quad1D, nf); - const int execMode = (int) assembly.GetExecMode(); + const int execMode = (int) trace_assembly.GetExecMode(); if (type == FaceType::Interior) { @@ -699,27 +702,27 @@ void PAResidualDistribution::SetupPA3D(FaceType type) const }//bdry } -void PAResidualDistribution::ApplyFaceTerms(const Vector &x, Vector &y, +void PADGTraceLOSolver::ApplyFaceTerms(const Vector &x, Vector &y, FaceType type) const { - Mesh *mesh = pfes.GetMesh(); + Mesh *mesh = trace_pfes.GetMesh(); int dim = mesh->Dimension(); if (dim == 2) { return ApplyFaceTerms2D(x, y, type); } if (dim == 3) { return ApplyFaceTerms3D(x, y, type); } } -void PAResidualDistribution::ApplyFaceTerms2D(const Vector &x, Vector &y, +void PADGTraceLOSolver::ApplyFaceTerms2D(const Vector &x, Vector &y, FaceType type) const { const int Q1D = quad1D; const int D1D = dofs1D; - const int nf = pfes.GetNFbyType(type); + const int nf = trace_pfes.GetNFbyType(type); const FaceRestriction * face_restrict_lex = nullptr; - const IntegrationRule *ir = assembly.lom.irF; + const IntegrationRule *ir = trace_assembly.lom.irF; const FiniteElement &el_trace = - *pfes.GetTraceElement(0, pfes.GetMesh()->GetFaceBaseGeometry(0)); + *trace_pfes.GetTraceElement(0, trace_pfes.GetMesh()->GetFaceBaseGeometry(0)); const DofToQuad *maps = &el_trace.GetDofToQuad(*ir, DofToQuad::TENSOR); auto B = mfem::Reshape(maps->B.Read(), quad1D, dofs1D); @@ -727,7 +730,7 @@ void PAResidualDistribution::ApplyFaceTerms2D(const Vector &x, Vector &y, if (type == FaceType::Interior) { face_restrict_lex = - pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type); + trace_pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type); Vector x_loc(face_restrict_lex->Height()); Vector y_loc(face_restrict_lex->Height()); @@ -794,7 +797,7 @@ void PAResidualDistribution::ApplyFaceTerms2D(const Vector &x, Vector &y, { face_restrict_lex = - pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type, + trace_pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type, L2FaceValues::SingleValued); Vector x_loc(face_restrict_lex->Height()); @@ -850,17 +853,17 @@ void PAResidualDistribution::ApplyFaceTerms2D(const Vector &x, Vector &y, } } -void PAResidualDistribution::ApplyFaceTerms3D(const Vector &x, Vector &y, - FaceType type) const +void PADGTraceLOSolver::ApplyFaceTerms3D(const Vector &x, Vector &y, + FaceType type) const { const int Q1D = quad1D; const int D1D = dofs1D; - const int nf = pfes.GetNFbyType(type); + const int nf = trace_pfes.GetNFbyType(type); const FaceRestriction * face_restrict_lex = nullptr; - const IntegrationRule *ir = assembly.lom.irF; + const IntegrationRule *ir = trace_assembly.lom.irF; const FiniteElement &el_trace = - *pfes.GetTraceElement(0, pfes.GetMesh()->GetFaceBaseGeometry(0)); + *trace_pfes.GetTraceElement(0, trace_pfes.GetMesh()->GetFaceBaseGeometry(0)); const DofToQuad *maps = &el_trace.GetDofToQuad(*ir, DofToQuad::TENSOR); auto B = mfem::Reshape(maps->B.Read(), quad1D, dofs1D); @@ -868,7 +871,7 @@ void PAResidualDistribution::ApplyFaceTerms3D(const Vector &x, Vector &y, if (type == FaceType::Interior) { face_restrict_lex = - pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type); + trace_pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type); Vector x_loc(face_restrict_lex->Height()); Vector y_loc(face_restrict_lex->Height()); @@ -968,7 +971,7 @@ void PAResidualDistribution::ApplyFaceTerms3D(const Vector &x, Vector &y, { face_restrict_lex = - pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type, + trace_pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type, L2FaceValues::SingleValued); Vector x_loc(face_restrict_lex->Height()); diff --git a/remhos_lo.hpp b/remhos_lo.hpp index cf042cf..d167e6f 100644 --- a/remhos_lo.hpp +++ b/remhos_lo.hpp @@ -38,6 +38,38 @@ class LOSolver class Assembly; +class PADGTraceLOSolver +{ +protected: + // Data at quadrature points + ParFiniteElementSpace &trace_pfes; + const int quad1D, dofs1D, face_dofs; + mutable Array D_int, D_bdry; + mutable Array IntVelocity, BdryVelocity; + Assembly &trace_assembly; + +public: + PADGTraceLOSolver(ParFiniteElementSpace &pfes_, const int q1D, + const int d1D, const int fdofs, + Assembly &assembly_) : + trace_pfes(pfes_), quad1D(q1D), dofs1D(d1D), face_dofs(fdofs), + trace_assembly(assembly_) {} + + void SampleVelocity(FaceType type) const; + + void SetupPA(FaceType type) const; + + void SetupPA2D(FaceType) const; + + void SetupPA3D(FaceType) const; + + void ApplyFaceTerms(const Vector &x, Vector &y, FaceType type) const; + + void ApplyFaceTerms2D(const Vector &x, Vector &y, FaceType type) const; + + void ApplyFaceTerms3D(const Vector &x, Vector &y, FaceType type) const; +}; + class DiscreteUpwind : public LOSolver { protected: @@ -84,20 +116,6 @@ class PADiscreteUpwind : public LOSolver const Array &adv_smap, const Vector &Mlump, Assembly &asmbly, bool updateD); - void SampleVelocity(FaceType type) const; - - void SetupPA(FaceType type) const; - - void SetupPA2D(FaceType) const; - - void SetupPA3D(FaceType) const; - - void ApplyFaceTerms(const Vector &x, Vector &y, FaceType type) const; - - void ApplyFaceTerms2D(const Vector &x, Vector &y, FaceType type) const; - - void ApplyFaceTerms3D(const Vector &x, Vector &y, FaceType type) const; - void ComputeAlgebraicDiffusion(Vector &ConvMats, Vector &AlgDiff) const; void AddBlkMult(const Vector &Mat, const Vector &x, Vector &y) const; @@ -123,33 +141,15 @@ class ResidualDistribution : public LOSolver }; //PA based Residual Distribution -class PAResidualDistribution : public ResidualDistribution +class PAResidualDistribution : public ResidualDistribution , public PADGTraceLOSolver { protected: - // Data at quadrature points - const int quad1D, dofs1D, face_dofs; - mutable Array D_int, D_bdry; - mutable Array IntVelocity, BdryVelocity; public: PAResidualDistribution(ParFiniteElementSpace &space, ParBilinearForm &Kbf, Assembly &asmbly, const Vector &Mlump, bool subcell, bool timedep); - void SampleVelocity(FaceType type) const; - - void SetupPA(FaceType type) const; - - void SetupPA2D(FaceType) const; - - void SetupPA3D(FaceType) const; - - void ApplyFaceTerms(const Vector &x, Vector &y, FaceType type) const; - - void ApplyFaceTerms2D(const Vector &x, Vector &y, FaceType type) const; - - void ApplyFaceTerms3D(const Vector &x, Vector &y, FaceType type) const; - virtual void CalcLOSolution(const Vector &u, Vector &du) const; }; From 42f3bc0a7737803e484147d5173f6b42cde0f1e2 Mon Sep 17 00:00:00 2001 From: Arturo Vargas Date: Mon, 20 Dec 2021 17:45:53 -0800 Subject: [PATCH 3/7] functional code -- need to clean up --- remhos.cpp | 17 ++++++++ remhos_lo.cpp | 106 ++++++++++++++++++++++---------------------------- remhos_lo.hpp | 27 +++++-------- 3 files changed, 74 insertions(+), 76 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index c9db842..a28c8f2 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -624,6 +624,16 @@ int main(int argc, char *argv[]) lo_smap = SparseMatrix_Build_smap(k.SpMat()); lo_solver = new PADiscreteUpwind(pfes, k_du, conv_int, k.SpMat(), lo_smap, lumpedM, asmbl, time_dep); + + if (exec_mode == 0) + { + const PADiscreteUpwind *lo_ptr = + dynamic_cast(lo_solver); + lo_ptr->SampleVelocity(FaceType::Interior); + lo_ptr->SampleVelocity(FaceType::Boundary); + lo_ptr->SetupPA(FaceType::Interior); + lo_ptr->SetupPA(FaceType::Boundary); + } } else if (lo_type == LOSolverType::DiscrUpwindPrec) { @@ -1219,6 +1229,13 @@ void AdvectionOperator::Mult(const Vector &X, Vector &Y) const RD_ptr->SetupPA(FaceType::Interior); RD_ptr->SetupPA(FaceType::Boundary); } + else if (auto lo_ptr = dynamic_cast(lo_solver)) + { + lo_ptr->SampleVelocity(FaceType::Interior); + lo_ptr->SampleVelocity(FaceType::Boundary); + lo_ptr->SetupPA(FaceType::Interior); + lo_ptr->SetupPA(FaceType::Boundary); + } else { for (int k = 0; k < ne; k++) diff --git a/remhos_lo.cpp b/remhos_lo.cpp index 98619e9..0196377 100644 --- a/remhos_lo.cpp +++ b/remhos_lo.cpp @@ -45,16 +45,17 @@ PADiscreteUpwind::PADiscreteUpwind(ParFiniteElementSpace &space, ParBilinearForm &Kbf, ConvectionIntegrator *Conv_, const SparseMatrix &adv, const Array &adv_smap, const Vector &Mlump, - Assembly &asmbly, bool updateD) + Assembly &asmbly, bool timedep) : LOSolver(space), + PADGTraceLOSolver(space, get_maps(pfes, asmbly)->nqpt, + get_maps(pfes, asmbly)->ndof, + (pfes.GetMesh()->Dimension() ==2) ? + get_maps(pfes, asmbly)->nqpt : + get_maps(pfes, asmbly)->nqpt * get_maps(pfes, asmbly)->nqpt, + asmbly), k_pbilinear(Kbf), Conv(Conv_), K(adv), D(), K_smap(adv_smap), M_lumped(Mlump), - assembly(asmbly), update_D(updateD), - quad1D(get_maps(pfes, assembly)->nqpt), - dofs1D(get_maps(pfes, assembly)->ndof), - face_dofs((pfes.GetMesh()->Dimension() ==2) ? quad1D : quad1D * quad1D) + time_dep(timedep) { - D = K; - ComputeDiscreteUpwindMatrix(); } void DiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const @@ -93,18 +94,15 @@ void DiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const void PADiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const { //mfem_error("Here! \n"); - const int ndof = pfes.GetFE(0)->GetDof(); - const int NE = pfes.GetMesh()->GetNE(); + const int ndof = trace_pfes.GetFE(0)->GetDof(); + const int NE = trace_pfes.GetMesh()->GetNE(); Vector alpha(ndof); alpha = 0.0; - - // Recompute D due to mesh changes (K changes) in remap mode. - if (update_D) { ComputeDiscreteUpwindMatrix(); } + du = 0.0; // Discretization and monotonicity terms. - D.Mult(u, du); - //K.Mult(u, du); + //D.Mult(u, du); - Vector y(du.Size()); y = 0.0; + //Vector y(du.Size()); y = 0.0; //Clearly incorrect... //k_pbilinear.Mult(u, y); @@ -117,11 +115,19 @@ void PADiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const ComputeAlgebraicDiffusion(ConvMats, AlgDiffMats); - AddBlkMult(ConvMats, u, y); + //K* = K + D + + //Mult K - AddBlkMult(AlgDiffMats, u, y); + AddBlkMult(ConvMats, u, du); + //Mult D + AddBlkMult(AlgDiffMats, u, du); + + //std::cout< 1e-12) @@ -129,24 +135,32 @@ void PADiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const std::cout<<"Error too high "<nqpt, - get_maps(pfes, trace_assembly)->ndof, - (pfes.GetMesh()->Dimension() ==2) ? - get_maps(pfes, trace_assembly)->nqpt : - get_maps(pfes, trace_assembly)->nqpt * get_maps(pfes, trace_assembly)->nqpt, - assembly) + get_maps(pfes, trace_assembly)->ndof, + (pfes.GetMesh()->Dimension() ==2) ? + get_maps(pfes, trace_assembly)->nqpt : + get_maps(pfes, trace_assembly)->nqpt * get_maps(pfes, trace_assembly)->nqpt, + assembly) { } @@ -703,7 +691,7 @@ void PADGTraceLOSolver::SetupPA3D(FaceType type) const } void PADGTraceLOSolver::ApplyFaceTerms(const Vector &x, Vector &y, - FaceType type) const + FaceType type) const { Mesh *mesh = trace_pfes.GetMesh(); int dim = mesh->Dimension(); @@ -713,7 +701,7 @@ void PADGTraceLOSolver::ApplyFaceTerms(const Vector &x, Vector &y, } void PADGTraceLOSolver::ApplyFaceTerms2D(const Vector &x, Vector &y, - FaceType type) const + FaceType type) const { const int Q1D = quad1D; const int D1D = dofs1D; @@ -798,7 +786,7 @@ void PADGTraceLOSolver::ApplyFaceTerms2D(const Vector &x, Vector &y, face_restrict_lex = trace_pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type, - L2FaceValues::SingleValued); + L2FaceValues::SingleValued); Vector x_loc(face_restrict_lex->Height()); Vector y_loc(face_restrict_lex->Height()); @@ -854,7 +842,7 @@ void PADGTraceLOSolver::ApplyFaceTerms2D(const Vector &x, Vector &y, } void PADGTraceLOSolver::ApplyFaceTerms3D(const Vector &x, Vector &y, - FaceType type) const + FaceType type) const { const int Q1D = quad1D; const int D1D = dofs1D; @@ -972,7 +960,7 @@ void PADGTraceLOSolver::ApplyFaceTerms3D(const Vector &x, Vector &y, face_restrict_lex = trace_pfes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,type, - L2FaceValues::SingleValued); + L2FaceValues::SingleValued); Vector x_loc(face_restrict_lex->Height()); Vector y_loc(face_restrict_lex->Height()); diff --git a/remhos_lo.hpp b/remhos_lo.hpp index d167e6f..4d88e37 100644 --- a/remhos_lo.hpp +++ b/remhos_lo.hpp @@ -49,11 +49,11 @@ class PADGTraceLOSolver Assembly &trace_assembly; public: - PADGTraceLOSolver(ParFiniteElementSpace &pfes_, const int q1D, - const int d1D, const int fdofs, - Assembly &assembly_) : - trace_pfes(pfes_), quad1D(q1D), dofs1D(d1D), face_dofs(fdofs), - trace_assembly(assembly_) {} + PADGTraceLOSolver(ParFiniteElementSpace &pfes_, const int q1D, + const int d1D, const int fdofs, + Assembly &assembly_) : + trace_pfes(pfes_), quad1D(q1D), dofs1D(d1D), face_dofs(fdofs), + trace_assembly(assembly_) {} void SampleVelocity(FaceType type) const; @@ -90,7 +90,7 @@ class DiscreteUpwind : public LOSolver virtual void CalcLOSolution(const Vector &u, Vector &du) const; }; -class PADiscreteUpwind : public LOSolver +class PADiscreteUpwind : public LOSolver, public PADGTraceLOSolver { protected: ParBilinearForm &k_pbilinear; @@ -99,22 +99,14 @@ class PADiscreteUpwind : public LOSolver mutable SparseMatrix D; const Array &K_smap; const Vector &M_lumped; - Assembly &assembly; - const bool update_D; - - void ComputeDiscreteUpwindMatrix() const; - - // Data at quadrature points //for face terms - const int quad1D, dofs1D, face_dofs; - mutable Array D_int, D_bdry; - mutable Array IntVelocity, BdryVelocity; + const bool time_dep; public: PADiscreteUpwind(ParFiniteElementSpace &space, ParBilinearForm &Kbf, ConvectionIntegrator *Conv_, const SparseMatrix &adv, const Array &adv_smap, const Vector &Mlump, - Assembly &asmbly, bool updateD); + Assembly &asmbly, bool timedep); void ComputeAlgebraicDiffusion(Vector &ConvMats, Vector &AlgDiff) const; @@ -141,7 +133,8 @@ class ResidualDistribution : public LOSolver }; //PA based Residual Distribution -class PAResidualDistribution : public ResidualDistribution , public PADGTraceLOSolver +class PAResidualDistribution : public ResidualDistribution, + public PADGTraceLOSolver { protected: From be1d046912f17859dbd28fca1f9252d67d328afd Mon Sep 17 00:00:00 2001 From: Arturo Vargas Date: Mon, 20 Dec 2021 18:04:36 -0800 Subject: [PATCH 4/7] clean up pass -- need to remove commented code now --- remhos.cpp | 30 ++++++++++++------- remhos_lo.cpp | 82 +++++++++++++++++---------------------------------- remhos_lo.hpp | 21 ++++++++----- 3 files changed, 59 insertions(+), 74 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index a28c8f2..f17c911 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -621,18 +621,26 @@ int main(int argc, char *argv[]) const bool time_dep = (exec_mode == 0) ? false : true; if (lo_type == LOSolverType::DiscrUpwind) { - lo_smap = SparseMatrix_Build_smap(k.SpMat()); - lo_solver = new PADiscreteUpwind(pfes, k_du, conv_int, k.SpMat(), lo_smap, - lumpedM, asmbl, time_dep); - - if (exec_mode == 0) + if (pa) { - const PADiscreteUpwind *lo_ptr = - dynamic_cast(lo_solver); - lo_ptr->SampleVelocity(FaceType::Interior); - lo_ptr->SampleVelocity(FaceType::Boundary); - lo_ptr->SetupPA(FaceType::Interior); - lo_ptr->SetupPA(FaceType::Boundary); + lo_solver = new PADiscreteUpwind(pfes, conv_int, + lumpedM, asmbl, time_dep); + if (exec_mode == 0) + { + const PADiscreteUpwind *lo_ptr = + dynamic_cast(lo_solver); + lo_ptr->AssembleBlkOperators(); + lo_ptr->SampleVelocity(FaceType::Interior); + lo_ptr->SampleVelocity(FaceType::Boundary); + lo_ptr->SetupPA(FaceType::Interior); + lo_ptr->SetupPA(FaceType::Boundary); + } + } + else + { + lo_smap = SparseMatrix_Build_smap(k.SpMat()); + lo_solver = new DiscreteUpwind(pfes, k.SpMat(), lo_smap, + lumpedM, asmbl, time_dep); } } else if (lo_type == LOSolverType::DiscrUpwindPrec) diff --git a/remhos_lo.cpp b/remhos_lo.cpp index 0196377..4813bac 100644 --- a/remhos_lo.cpp +++ b/remhos_lo.cpp @@ -42,9 +42,9 @@ DiscreteUpwind::DiscreteUpwind(ParFiniteElementSpace &space, } PADiscreteUpwind::PADiscreteUpwind(ParFiniteElementSpace &space, - ParBilinearForm &Kbf, ConvectionIntegrator *Conv_, - const SparseMatrix &adv, - const Array &adv_smap, const Vector &Mlump, + /*ParBilinearForm &Kbf,*/ ConvectionIntegrator *Conv_, + /*const SparseMatrix &adv,*/ + /*const Array &adv_smap,*/ const Vector &Mlump, Assembly &asmbly, bool timedep) : LOSolver(space), PADGTraceLOSolver(space, get_maps(pfes, asmbly)->nqpt, @@ -53,7 +53,8 @@ PADiscreteUpwind::PADiscreteUpwind(ParFiniteElementSpace &space, get_maps(pfes, asmbly)->nqpt : get_maps(pfes, asmbly)->nqpt * get_maps(pfes, asmbly)->nqpt, asmbly), - k_pbilinear(Kbf), Conv(Conv_), K(adv), D(), K_smap(adv_smap), M_lumped(Mlump), + /*k_pbilinear(Kbf),*/ Conv(Conv_), /*K(adv), D(), K_smap(adv_smap), */ M_lumped( + Mlump), time_dep(timedep) { } @@ -93,79 +94,50 @@ void DiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const void PADiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const { - //mfem_error("Here! \n"); - const int ndof = trace_pfes.GetFE(0)->GetDof(); - const int NE = trace_pfes.GetMesh()->GetNE(); - Vector alpha(ndof); alpha = 0.0; du = 0.0; - // Discretization and monotonicity terms. - //D.Mult(u, du); - - //Vector y(du.Size()); y = 0.0; - - //Clearly incorrect... - //k_pbilinear.Mult(u, y); - - Vector ConvMats(ndof * ndof * NE); - Vector AlgDiffMats(ndof * ndof * NE); AlgDiffMats = 0.0; - - //Given in row major format - Conv->AssembleEA(pfes, ConvMats, false); - - ComputeAlgebraicDiffusion(ConvMats, AlgDiffMats); + //Construct K, D + if (false) + { + //Given in row major format + AssembleBlkOperators(); + } - //K* = K + D + //Apply K* = K + D //Mult K - AddBlkMult(ConvMats, u, du); //Mult D AddBlkMult(AlgDiffMats, u, du); - - //std::cout< 1e-12) - { - std::cout<<"Error too high "<GetDof(); + const int NE = trace_pfes.GetMesh()->GetNE(); + + ConvMats.SetSize(ndof * ndof * NE); + AlgDiffMats.SetSize(ndof * ndof * NE); AlgDiffMats = 0.0; + + Conv->AssembleEA(pfes, ConvMats, false); + ComputeAlgebraicDiffusion(ConvMats, AlgDiffMats); +} + void PADiscreteUpwind::ComputeAlgebraicDiffusion(Vector &ConvMats, Vector &AlgDiff) const { diff --git a/remhos_lo.hpp b/remhos_lo.hpp index 4d88e37..df96e6b 100644 --- a/remhos_lo.hpp +++ b/remhos_lo.hpp @@ -93,21 +93,26 @@ class DiscreteUpwind : public LOSolver class PADiscreteUpwind : public LOSolver, public PADGTraceLOSolver { protected: - ParBilinearForm &k_pbilinear; - ConvectionIntegrator *Conv; - const SparseMatrix &K; - mutable SparseMatrix D; - const Array &K_smap; + //ParBilinearForm &k_pbilinear; + ConvectionIntegrator *Conv = nullptr; + //const SparseMatrix &K; + //mutable SparseMatrix D; + //const Array &K_smap; const Vector &M_lumped; const bool time_dep; + mutable Vector ConvMats; + mutable Vector AlgDiffMats; public: PADiscreteUpwind(ParFiniteElementSpace &space, - ParBilinearForm &Kbf, ConvectionIntegrator *Conv_, - const SparseMatrix &adv, - const Array &adv_smap, const Vector &Mlump, + /*ParBilinearForm &Kbf,*/ ConvectionIntegrator *Conv_, + /* const SparseMatrix &adv,*/ + /*const Array &adv_smap, */ const Vector &Mlump, Assembly &asmbly, bool timedep); + /*Block Operators are in row major format*/ + void AssembleBlkOperators() const; + void ComputeAlgebraicDiffusion(Vector &ConvMats, Vector &AlgDiff) const; void AddBlkMult(const Vector &Mat, const Vector &x, Vector &y) const; From bcddb191f5842e8e832b651ef394903c25b657be Mon Sep 17 00:00:00 2001 From: Arturo Vargas Date: Mon, 20 Dec 2021 18:08:22 -0800 Subject: [PATCH 5/7] ready for review --- remhos.cpp | 6 ------ remhos_lo.cpp | 9 +++------ remhos_lo.hpp | 9 ++------- 3 files changed, 5 insertions(+), 19 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index f17c911..811968a 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -390,21 +390,18 @@ int main(int argc, char *argv[]) M_HO.AddDomainIntegrator(new MassIntegrator); ParBilinearForm k(&pfes); - ParBilinearForm k_du(&pfes); ParBilinearForm K_HO(&pfes); ConvectionIntegrator *conv_int = nullptr; if (exec_mode == 0) { conv_int = new ConvectionIntegrator(velocity, -1.0); k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0)); - k_du.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0)); K_HO.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0)); } else if (exec_mode == 1) { conv_int = new ConvectionIntegrator(v_coef); k.AddDomainIntegrator(conv_int); - k_du.AddDomainIntegrator(conv_int); K_HO.AddDomainIntegrator(conv_int); } @@ -436,14 +433,11 @@ int main(int argc, char *argv[]) K_HO.SetAssemblyLevel(AssemblyLevel::PARTIAL); } - k_du.SetAssemblyLevel(AssemblyLevel::PARTIAL); - if (next_gen_full) { K_HO.SetAssemblyLevel(AssemblyLevel::FULL); } - k_du.Assemble(); M_HO.Assemble(); K_HO.Assemble(0); diff --git a/remhos_lo.cpp b/remhos_lo.cpp index 4813bac..ab6c7e7 100644 --- a/remhos_lo.cpp +++ b/remhos_lo.cpp @@ -42,9 +42,8 @@ DiscreteUpwind::DiscreteUpwind(ParFiniteElementSpace &space, } PADiscreteUpwind::PADiscreteUpwind(ParFiniteElementSpace &space, - /*ParBilinearForm &Kbf,*/ ConvectionIntegrator *Conv_, - /*const SparseMatrix &adv,*/ - /*const Array &adv_smap,*/ const Vector &Mlump, + ConvectionIntegrator *Conv_, + const Vector &Mlump, Assembly &asmbly, bool timedep) : LOSolver(space), PADGTraceLOSolver(space, get_maps(pfes, asmbly)->nqpt, @@ -53,9 +52,7 @@ PADiscreteUpwind::PADiscreteUpwind(ParFiniteElementSpace &space, get_maps(pfes, asmbly)->nqpt : get_maps(pfes, asmbly)->nqpt * get_maps(pfes, asmbly)->nqpt, asmbly), - /*k_pbilinear(Kbf),*/ Conv(Conv_), /*K(adv), D(), K_smap(adv_smap), */ M_lumped( - Mlump), - time_dep(timedep) + Conv(Conv_), M_lumped(Mlump), time_dep(timedep) { } diff --git a/remhos_lo.hpp b/remhos_lo.hpp index df96e6b..791b564 100644 --- a/remhos_lo.hpp +++ b/remhos_lo.hpp @@ -93,11 +93,7 @@ class DiscreteUpwind : public LOSolver class PADiscreteUpwind : public LOSolver, public PADGTraceLOSolver { protected: - //ParBilinearForm &k_pbilinear; ConvectionIntegrator *Conv = nullptr; - //const SparseMatrix &K; - //mutable SparseMatrix D; - //const Array &K_smap; const Vector &M_lumped; const bool time_dep; mutable Vector ConvMats; @@ -105,9 +101,8 @@ class PADiscreteUpwind : public LOSolver, public PADGTraceLOSolver public: PADiscreteUpwind(ParFiniteElementSpace &space, - /*ParBilinearForm &Kbf,*/ ConvectionIntegrator *Conv_, - /* const SparseMatrix &adv,*/ - /*const Array &adv_smap, */ const Vector &Mlump, + ConvectionIntegrator *Conv_, + const Vector &Mlump, Assembly &asmbly, bool timedep); /*Block Operators are in row major format*/ From 01116604f39255f1b7e605b049f609a501c7a0f3 Mon Sep 17 00:00:00 2001 From: Arturo Vargas Date: Mon, 20 Dec 2021 18:14:35 -0800 Subject: [PATCH 6/7] clean up pass --- remhos_lo.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/remhos_lo.cpp b/remhos_lo.cpp index ab6c7e7..cf6cf30 100644 --- a/remhos_lo.cpp +++ b/remhos_lo.cpp @@ -156,10 +156,6 @@ void PADiscreteUpwind::ComputeAlgebraicDiffusion(Vector &ConvMats, if (r == c) { continue; } - //for(int k=0; k Date: Mon, 20 Dec 2021 19:51:28 -0800 Subject: [PATCH 7/7] fix segfault bugs --- remhos.cpp | 6 ++++-- remhos_lo.cpp | 7 +------ 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 811968a..46d99a0 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -401,8 +401,8 @@ int main(int argc, char *argv[]) else if (exec_mode == 1) { conv_int = new ConvectionIntegrator(v_coef); - k.AddDomainIntegrator(conv_int); - K_HO.AddDomainIntegrator(conv_int); + k.AddDomainIntegrator(new ConvectionIntegrator(v_coef)); + K_HO.AddDomainIntegrator(new ConvectionIntegrator(v_coef)); } if (ho_type == HOSolverType::CG || @@ -1233,6 +1233,8 @@ void AdvectionOperator::Mult(const Vector &X, Vector &Y) const } else if (auto lo_ptr = dynamic_cast(lo_solver)) { + //Construct K, D in K* = K + D + lo_ptr->AssembleBlkOperators(); lo_ptr->SampleVelocity(FaceType::Interior); lo_ptr->SampleVelocity(FaceType::Boundary); lo_ptr->SetupPA(FaceType::Interior); diff --git a/remhos_lo.cpp b/remhos_lo.cpp index cf6cf30..5867144 100644 --- a/remhos_lo.cpp +++ b/remhos_lo.cpp @@ -93,12 +93,7 @@ void PADiscreteUpwind::CalcLOSolution(const Vector &u, Vector &du) const { du = 0.0; - //Construct K, D - if (false) - { - //Given in row major format - AssembleBlkOperators(); - } + //Apply K* = K + D