AdvDiffSolver::GetFlowField and AdvDiffFlowPastArray class. still nee…

…ds debugging.

include/advdiff_solver.hpp

@@ -27,6 +27,9 @@ friend class ParameterizedProblem;
       flow solver to obtain the prescribed velocity field. both StokesSolver / SteadyNSSolver can be used.
    */
    StokesSolver *stokes_solver = NULL;
+   bool load_flow = false;
+   bool save_flow = false;
+   std::string flow_file = "";
 
 public:
    AdvDiffSolver();
@@ -46,6 +49,9 @@ friend class ParameterizedProblem;
 
 protected:
    void SetMUMPSSolver() override;
+
+private:
+   void GetFlowField(ParameterizedProblem *flow_problem);
 };
 
 #endif

include/parameterized_problem.hpp

@@ -92,9 +92,24 @@ extern double rforce_f;
 void tip_force(const Vector &x, Vector &u);
 }
 
+namespace advdiff_problem
+{
+
+extern bool analytic_flow;
+
+namespace advdiff_flow_past_array
+{
+   extern double q0, dq, qoffset;
+   extern Vector qk;
+
+   double qbdr(const Vector &x);
+}  //  namespace advdiff_flow_past_array
+
+}  // namespace advdiff_problem
+
 }
 
-enum BoundaryType
+enum class BoundaryType
 { 
    ZERO,
    DIRICHLET,
@@ -169,8 +184,8 @@ class PoissonComponent : public PoissonProblem
 public:
    PoissonComponent();
    virtual ~PoissonComponent() {};
-   virtual void SetParams(const std::string &key, const double &value) override;
-   virtual void SetParams(const Array<int> &indexes, const Vector &values) override;
+   void SetParams(const std::string &key, const double &value) override;
+   void SetParams(const Array<int> &indexes, const Vector &values) override;
 
 private:
    void SetBattr();
@@ -206,15 +221,17 @@ class StokesComponent : public StokesProblem
 
 class StokesFlowPastArray : public StokesComponent
 {
+friend class AdvDiffFlowPastArray;
+
 public:
    StokesFlowPastArray();
 
-   virtual void SetParams(const std::string &key, const double &value);
-   virtual void SetParams(const Array<int> &indexes, const Vector &values);
+   virtual void SetParams(const std::string &key, const double &value) override;
+   virtual void SetParams(const Array<int> &indexes, const Vector &values) override;
 
-private:
+protected:
    Vector *u0;
-   void SetBattr();
+   virtual void SetBattr();
 };
 
 class LinElastProblem : public ParameterizedProblem
@@ -251,6 +268,43 @@ class LinElastForceCantilever : public LinElastProblem
    LinElastForceCantilever();
 };
 
+namespace function_factory
+{
+
+namespace advdiff_problem
+{
+/*
+   flow_problem will be passed down to StokesSolver/SteadyNSSolver for obtaining velocity field.
+   It must be set appropriately within each AdvDiffSolver problems.
+*/
+extern StokesProblem *flow_problem;
+
+}  // namespace advdiff_problem
+
+}  // namespace function_factory
+
+class AdvDiffFlowPastArray : public StokesFlowPastArray
+{
+protected:
+   /*
+      flow_problem shares the same pointers with this class.
+      Thus every parameter set by this class is reflected to StokesFlowPastArrayProblem as well.
+      flow_problem will be passed down to StokesSolver/SteadyNSSolver for obtaining velocity field.
+   */
+   StokesFlowPastArray *flow_problem = NULL;
+
+public:
+   AdvDiffFlowPastArray();
+   virtual ~AdvDiffFlowPastArray();
+
+protected:
+   void SetBattr() override
+   {
+      StokesFlowPastArray::SetBattr();
+      flow_problem->SetBattr();
+   }
+};
+
 ParameterizedProblem* InitParameterizedProblem();
 
 #endif

src/advdiff_solver.cpp

@@ -19,11 +19,17 @@ AdvDiffSolver::AdvDiffSolver()
    Pe = config.GetOption<double>("adv-diff/peclet_number", 0.0);
 
    flow_coeffs.SetSize(numSub);
+
+   save_flow = config.GetOption<bool>("adv-diff/save_flow", false);
+   load_flow = config.GetOption<bool>("adv-diff/load_flow", false);
+   if (save_flow || load_flow)
+      flow_file = config.GetRequiredOption<std::string>("adv-diff/flow_file");
 }
 
 AdvDiffSolver::~AdvDiffSolver()
 {
    DeletePointers(flow_coeffs);
+   delete stokes_solver;
 }
 
 void AdvDiffSolver::BuildDomainOperators()
@@ -160,9 +166,10 @@ void AdvDiffSolver::SetFlowAtSubdomain(std::function<void(const Vector &, Vector
 
 void AdvDiffSolver::SetParameterizedProblem(ParameterizedProblem *problem)
 {
-   PoissonSolver::SetParameterizedProblem(problem);
+   if (!function_factory::advdiff_problem::analytic_flow)
+      GetFlowField(function_factory::advdiff_problem::flow_problem);
 
-   // TODO(kevin): add flow field setup.
+   PoissonSolver::SetParameterizedProblem(problem);
 }
 
 void AdvDiffSolver::SetMUMPSSolver()
@@ -171,4 +178,40 @@ void AdvDiffSolver::SetMUMPSSolver()
    mumps = new MUMPSSolver();
    mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
    mumps->SetOperator(*globalMat_hypre);
+}
+
+void AdvDiffSolver::GetFlowField(ParameterizedProblem *flow_problem)
+{
+   assert(flow_problem);
+   mfem_warning("AdvDiffSolver: Obtaining flow field. This may take a while depending on the domain size.\n");
+
+   stokes_solver = new StokesSolver;
+   stokes_solver->InitVariables();
+   stokes_solver->InitVisualization();
+
+   if (load_flow && FileExists(flow_file))
+      stokes_solver->LoadSolution(flow_file);
+   else
+   {
+      stokes_solver->SetParameterizedProblem(flow_problem);
+      // currently only support FOM.
+      stokes_solver->BuildOperators();
+      stokes_solver->SetupBCOperators();
+      stokes_solver->Assemble();
+      stokes_solver->Solve();
+   }
+
+   if (save_flow)
+      stokes_solver->SaveSolution(flow_file);
+
+   DeletePointers(flow_coeffs);
+   const int stokes_numvar = stokes_solver->GetNumVar();
+   for (int m = 0; m < numSub; m++)
+      flow_coeffs[m] = new VectorGridFunctionCoefficient(stokes_solver->GetGridFunction(m * stokes_numvar));
+
+   /*
+      VectorGridFunctionCoefficient does not own the grid function,
+      and it requires the grid function for its lifetime.
+      Thus stokes_solver will be deleted at ~AdvDiffSolver().
+   */
 }

src/main_workflow.cpp

@@ -9,6 +9,7 @@
 #include "linelast_solver.hpp"
 #include "stokes_solver.hpp"
 #include "steady_ns_solver.hpp"
+#include "advdiff_solver.hpp"
 #include "etc.hpp"
 #include <fstream>
 #include <iostream>
@@ -55,6 +56,7 @@ MultiBlockSolver* InitSolver()
    else if (solver_type == "stokes")   { solver = new StokesSolver; }
    else if (solver_type == "steady-ns")   { solver = new SteadyNSSolver; }
    else if (solver_type == "linelast")   { solver = new LinElastSolver; }
+   else if (solver_type == "adv-diff")   { solver = new AdvDiffSolver; }
    else
    {
       printf("Unknown MultiBlockSolver %s!\n", solver_type.c_str());

src/parameterized_problem.cpp

@@ -200,6 +200,32 @@ void tip_force(const Vector &x, Vector &f)
 
 }  // namespace linelast_force
 
+namespace advdiff_problem
+{
+
+bool analytic_flow;
+StokesProblem *flow_problem;
+
+namespace advdiff_flow_past_array
+{
+
+double q0, dq, qoffset;
+Vector qk;
+
+double qbdr(const Vector &x)
+{
+   assert(qk.Size() >= x.Size());
+   double tmp = 0.0;
+   for (int d = 0; d < x.Size(); d++)
+      tmp += advdiff_flow_past_array::qk(d) * x(d);
+   tmp += advdiff_flow_past_array::qoffset;
+   return advdiff_flow_past_array::q0 + advdiff_flow_past_array::dq * sin(2.0 * pi * tmp);
+}
+
+}  // namespace advdiff_flow_past_array
+
+}  // namespace advdiff_problem
+
 }  // namespace function_factory
 
 ParameterizedProblem::ParameterizedProblem()
@@ -303,6 +329,10 @@ ParameterizedProblem* InitParameterizedProblem()
    {
       problem = new LinElastForceCantilever();
    }
+   else if (problem_name == "advdiff_flow_past_array")
+   {
+      problem = new AdvDiffFlowPastArray();
+   }
    else
    {
       mfem_error("Unknown parameterized problem name!\n");
@@ -762,3 +792,51 @@ LinElastForceCantilever::LinElastForceCantilever()
    general_vector_ptr.SetSize(1);
    general_vector_ptr[0] = NULL;
 }
+
+/*
+   AdvDiffFlowPastArray
+*/
+
+AdvDiffFlowPastArray::AdvDiffFlowPastArray()
+   : StokesFlowPastArray(), flow_problem(new StokesFlowPastArray)
+{
+   function_factory::advdiff_problem::analytic_flow = false;
+   function_factory::advdiff_problem::flow_problem = flow_problem;
+
+   bdr_type.SetSize(5);
+   bdr_type = BoundaryType::DIRICHLET;
+   bdr_type[4] = BoundaryType::NEUMANN;
+
+   scalar_rhs_ptr = NULL;
+   scalar_bdr_ptr.SetSize(5);
+   scalar_bdr_ptr = &(function_factory::advdiff_problem::advdiff_flow_past_array::qbdr);
+
+   // q0 + dq + qoffset + qk(3)
+   param_num += 1 + 1 + 1 + 3;
+   const int p0 = flow_problem->GetNumParams();
+
+   param_map["q0"] = p0;
+   param_map["dq"] = p0 + 1;
+   param_map["qoffset"] = p0 + 2;
+   param_map["qk_x"] = p0 + 3;
+   param_map["qk_y"] = p0 + 4;
+   param_map["qk_z"] = p0 + 5;
+
+   // default values.
+   function_factory::advdiff_problem::advdiff_flow_past_array::q0 = 1.0;
+   function_factory::advdiff_problem::advdiff_flow_past_array::dq = 0.1;
+   function_factory::advdiff_problem::advdiff_flow_past_array::qoffset = 0.0;
+   function_factory::advdiff_problem::advdiff_flow_past_array::qk.SetSize(3);
+   function_factory::advdiff_problem::advdiff_flow_past_array::qk = 0.0;
+
+   param_ptr.Append(&(function_factory::advdiff_problem::advdiff_flow_past_array::q0));
+   param_ptr.Append(&(function_factory::advdiff_problem::advdiff_flow_past_array::dq));
+   param_ptr.Append(&(function_factory::advdiff_problem::advdiff_flow_past_array::qoffset));
+   for (int j = 0; j < 3; j++)
+      param_ptr.Append(&(function_factory::advdiff_problem::advdiff_flow_past_array::qk(j)));
+}
+
+AdvDiffFlowPastArray::~AdvDiffFlowPastArray()
+{
+   delete flow_problem;
+}

src/poisson_solver.cpp

@@ -554,6 +554,9 @@ void PoissonSolver::SanityCheckOnCoeffs()
 
 void PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem)
 {
+   /* set up boundary types */
+   MultiBlockSolver::SetParameterizedProblem(problem);
+
    // clean up rhs for parametrized problem.
    if (rhs_coeffs.Size() > 0)
    {

src/stokes_solver.cpp

@@ -1078,6 +1078,9 @@ void StokesSolver::SanityCheckOnCoeffs()
 
 void StokesSolver::SetParameterizedProblem(ParameterizedProblem *problem)
 {
+   /* set up boundary types */
+   MultiBlockSolver::SetParameterizedProblem(problem);
+
    nu = function_factory::stokes_problem::nu;
    delete nu_coeff;
    nu_coeff = new ConstantCoefficient(nu);

test/CMakeLists.txt

@@ -33,6 +33,7 @@ add_executable(test_workflow test_workflow.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 file(COPY inputs/test.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
 file(COPY inputs/stokes.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
 file(COPY inputs/steady_ns.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
+file(COPY inputs/advdiff.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
 file(COPY meshes/test.2x2.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes)
 file(COPY meshes/test.2x1.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes)