Add ability to linearly ramp up input power

pecos · Feb 3, 2023 · 91a7e15 · 91a7e15
1 parent d220e9c
commit 91a7e15
Show file tree

Hide file tree

Showing 3 changed files with 18 additions and 5 deletions.
diff --git a/src/cycle_avg_joule_coupling.cpp b/src/cycle_avg_joule_coupling.cpp
@@ -36,8 +36,12 @@
 #include "quasimagnetostatic.hpp"
 
 CycleAvgJouleCoupling::CycleAvgJouleCoupling(MPI_Session &mpi, string &inputFileName, TPS::Tps *tps, int max_out,
-                                             bool axisym, double input_power)
-    : mpi_(mpi), em_opt_(), max_outer_iters_(max_out), input_power_(input_power) {
+                                             bool axisym, double input_power, double initial_input_power)
+    : mpi_(mpi),
+      em_opt_(),
+      max_outer_iters_(max_out),
+      input_power_(input_power),
+      initial_input_power_(initial_input_power) {
   if (axisym) {
     qmsa_solver_ = new QuasiMagnetostaticSolverAxiSym(mpi, em_opt_, tps);
   } else {
@@ -242,6 +246,11 @@ void CycleAvgJouleCoupling::solve() {
   int curr_iter = flow_solver_->getCurrentIterations();
   flow_solver_->setMaximumIterations(curr_iter + increment);
 
+  double delta_power = 0;
+  if (input_power_ > 0) {
+    delta_power = (input_power_ - initial_input_power_) / max_outer_iters_;
+  }
+
   for (int outer_iters = 0; outer_iters < max_outer_iters_; outer_iters++) {
     // EM
     interpConductivityFromFlowToEM();
@@ -252,7 +261,8 @@ void CycleAvgJouleCoupling::solve() {
     }
 
     if (input_power_ > 0) {
-      const double ratio = input_power_ / tot_jh;
+      const double target_power = initial_input_power_ + (outer_iters + 1) * delta_power;
+      const double ratio = target_power / tot_jh;
       qmsa_solver_->scaleJouleHeating(ratio);
       const double upd_jh = qmsa_solver_->totalJouleHeating();
       if (mpi_.Root()) {

diff --git a/src/cycle_avg_joule_coupling.hpp b/src/cycle_avg_joule_coupling.hpp
@@ -65,10 +65,11 @@ class CycleAvgJouleCoupling : public TPS::Solver {
   int n_flow_interp_nodes_;
 
   double input_power_;
+  double initial_input_power_;
 
  public:
   CycleAvgJouleCoupling(MPI_Session &mpi, string &inputFileName, TPS::Tps *tps, int max_out, bool axisym,
-                        double input_power = -1.);
+                        double input_power = -1., double initial_input_power = -1.);
   ~CycleAvgJouleCoupling();
 
   void initializeInterpolationData();

diff --git a/src/tps.cpp b/src/tps.cpp
@@ -242,7 +242,9 @@ void Tps::chooseSolver() {
     getInput("solver/axisymmetric", axisym, false);
     double input_power;
     getInput("solver/input-power", input_power, -1.);
-    solver_ = new CycleAvgJouleCoupling(mpi_, iFile_, this, max_out, axisym, input_power);
+    double initial_input_power;
+    getInput("solver/initial-input-power", initial_input_power, input_power);
+    solver_ = new CycleAvgJouleCoupling(mpi_, iFile_, this, max_out, axisym, input_power, initial_input_power);
   } else if (input_solver_type_ == "coupled") {
     isFlowEMCoupledMode_ = true;
     grvy_printf(GRVY_ERROR, "\nSlow your roll.  Solid high-five for whoever implements this coupled solver mode!\n");