WIP (#75): Starting to generalize input opts for driving voltage

pecos · Oct 25, 2021 · 0e265cb · 0e265cb
1 parent 7ff67db
commit 0e265cb
Show file tree

Hide file tree

Showing 2 changed files with 29 additions and 10 deletions.
diff --git a/src/em_options.hpp b/src/em_options.hpp
@@ -68,6 +68,10 @@ class ElectromagneticOptions {
   mfem::Array<int> neumann_bc_attr;    /**< List of boundary attributes corresponding to Neumann boundary */
   double nd_conductivity;  /**< Non-dimensional conductivity: \frac{\sigma_0}{\omega \epsilon_0} */
   double nd_frequency;     /**< Non-dimensional frequency: \frac{\omega \ell}{c} */
+  int port0;
+  int port1;
+  double Vstat0_real;
+  double Vstat1_real;
 
   ElectromagneticOptions()
     :
@@ -78,7 +82,9 @@ class ElectromagneticOptions {
     yinterp_min(0.0), yinterp_max(1.0),
     top_only(false), bot_only(false),
     conductor_domains(0), neumann_bc_attr(0),
-    nd_conductivity(1e6), nd_frequency(0.001)
+    nd_conductivity(1e6), nd_frequency(0.001),
+    port0(0), port1(1),
+    Vstat0_real(1.0), Vstat1_real(0.0)
   { }
 
   void AddElectromagneticOptions(mfem::OptionsParser &args) {
@@ -117,6 +123,14 @@ class ElectromagneticOptions {
                    "Non-dimensional conductivity, sigma_0/(omega*epsilon_0) (SEQS solver only)");
     args.AddOption(&nd_frequency, "-f", "--eta",
                    "Non-dimensional (angular) frequency, (omega*ell)/c (SEQS solver only)");
+    args.AddOption(&port0, "-p0", "--port0",
+                   "Boundary attribute of port 0 (SEQS solver only)");
+    args.AddOption(&port1, "-p1", "--port1",
+                   "Boundary attribute of port 1 (SEQS solver only)");
+    args.AddOption(&Vstat0_real, "-Vs0r", "--Vstat0real",
+                   "Voltage (real) at port 0 (SEQS solver only)");
+    args.AddOption(&Vstat1_real, "-Vs1r", "--Vstat1real",
+                   "Voltage (real) at port 1 (SEQS solver only)");
   }
 
   void print(std::ostream &out) {
@@ -154,6 +168,8 @@ class ElectromagneticOptions {
       out << ") " << std::endl;
       out << "    nd_conductivity = " << nd_conductivity << std::endl;
       out << "    nd_frequency    = " << nd_frequency << std::endl;
+      out << "    port0           = " << port0 << std::endl;
+      out << "    port1           = " << port1 << std::endl;
     }
     out << std::endl;
   }

diff --git a/src/seqs_maxwell_frequency.cpp b/src/seqs_maxwell_frequency.cpp
@@ -212,14 +212,14 @@ void SeqsMaxwellFrequencySolver::Initialize() {
 
   pspace_->GetEssentialTrueDofs(ess_bdr, h1_ess_tdof_list_);
 
-  // TODO(trevilo): Make these user-specified inputs
+  // Set the boundary indicators for the ports
   port_0_.SetSize(pmesh_->bdr_attributes.Max());
   port_0_ = 0;
-  port_0_[2] = 1;
+  port_0_[em_opt_.port0-1] = 1;
 
   port_1_.SetSize(pmesh_->bdr_attributes.Max());
   port_1_ = 0;
-  port_1_[1] = 1;
+  port_1_[em_opt_.port1-1] = 1;
 
   // Finally, get list of essential dofs for psi, which includes those
   // in the conductor as well as the phi Dirichlet dofs
@@ -297,8 +297,6 @@ void SeqsMaxwellFrequencySolver::SolveSEQS() {
   *psi_real_ = 0.0;
   *psi_imag_ = 0.0;
 
-  // Offset functions (e.g. \Phi_1 from Ostrowski and Hiptmair (eqn 5.2))
-  // TODO(trevilo): Assumed real and = 1 here... generalize!
   ConstantCoefficient one(1.0);
 
   V0_ = new ParGridFunction(pspace_);
@@ -325,21 +323,24 @@ void SeqsMaxwellFrequencySolver::SolveSEQS() {
     Kpp_real->AddDomainIntegrator(new DiffusionIntegrator(*rel_sig_));
     Kpp_real->Assemble();
     Kpp_real->Finalize();
-    Kpp_real->AddMult(*V0_, *bp_real, -1.0);
+    Kpp_real->AddMult(*V0_, *bp_real, -em_opt_.Vstat0_real);
+    Kpp_real->AddMult(*V1_, *bp_real, -em_opt_.Vstat1_real);
     delete Kpp_real;
 
     ParBilinearForm *Kpp_imag = new ParBilinearForm(pspace_);
     Kpp_imag->AddDomainIntegrator(new DiffusionIntegrator(*one_over_sigma_));
     Kpp_imag->Assemble();
     Kpp_imag->Finalize();
-    Kpp_imag->AddMult(*V0_, *bp_imag, -1.0);
+    Kpp_imag->AddMult(*V0_, *bp_imag, -em_opt_.Vstat0_real);
+    Kpp_imag->AddMult(*V1_, *bp_imag, -em_opt_.Vstat1_real);
     delete Kpp_imag;
 
     ParBilinearForm *Kss_real = new ParBilinearForm(pspace_);
     Kss_real->AddDomainIntegrator(new DiffusionIntegrator(*rel_eps_nc_));
     Kss_real->Assemble();
     Kss_real->Finalize();
-    Kss_real->AddMult(*V0_, *bs_real, -1.0);
+    Kss_real->AddMult(*V0_, *bs_real, -em_opt_.Vstat0_real);
+    Kss_real->AddMult(*V1_, *bs_real, -em_opt_.Vstat1_real);
     delete Kss_real;
   }
 
@@ -576,7 +577,9 @@ void SeqsMaxwellFrequencySolver::SolveSEQS() {
   // single function, which will be used by the magnetic vector
   // potential solve
   phi_tot_real_ = new ParGridFunction(pspace_);
-  *phi_tot_real_  = *V0_;
+  *phi_tot_real_ = 0.0;
+  add(*phi_tot_real_ , em_opt_.Vstat0_real, *V0_, *phi_tot_real_);
+  add(*phi_tot_real_ , em_opt_.Vstat1_real, *V1_, *phi_tot_real_);
   *phi_tot_real_ += *phi_real_;
   *phi_tot_real_ += *psi_real_;