Skip to content

Commit

Permalink
WIP (#75): Generalize drive voltage input
Browse files Browse the repository at this point in the history 
Allow complex input for driving port voltages.
  • Loading branch information
@trevilo
trevilo committed Oct 25, 2021
1 parent d43aaad commit 4907b2d
Show file tree
Hide file tree
Showing 4 changed files with 69 additions and 28 deletions.
11 changes: 10 additions & 1 deletion src/em_options.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,8 @@ class ElectromagneticOptions {
int port1;
double Vstat0_real;
double Vstat1_real;
double Vstat0_imag;
double Vstat1_imag;

ElectromagneticOptions()
:
Expand All @@ -84,7 +86,8 @@ class ElectromagneticOptions {
conductor_domains(0), neumann_bc_attr(0),
nd_conductivity(1e6), nd_frequency(0.001),
port0(0), port1(1),
Vstat0_real(1.0), Vstat1_real(0.0)
Vstat0_real(1.0), Vstat1_real(0.0),
Vstat0_imag(0.0), Vstat1_imag(0.0)
{ }

void AddElectromagneticOptions(mfem::OptionsParser &args) {
Expand Down Expand Up @@ -131,6 +134,10 @@ class ElectromagneticOptions {
"Voltage (real) at port 0 (SEQS solver only)");
args.AddOption(&Vstat1_real, "-Vs1r", "--Vstat1real",
"Voltage (real) at port 1 (SEQS solver only)");
args.AddOption(&Vstat0_imag, "-Vs0i", "--Vstat0imag",
"Voltage (imaginary) at port 0 (SEQS solver only)");
args.AddOption(&Vstat1_imag, "-Vs1i", "--Vstat1imag",
"Voltage (imaginary) at port 1 (SEQS solver only)");
}

void print(std::ostream &out) {
Expand Down Expand Up @@ -170,6 +177,8 @@ class ElectromagneticOptions {
out << " nd_frequency = " << nd_frequency << std::endl;
out << " port0 = " << port0 << std::endl;
out << " port1 = " << port1 << std::endl;
out << " V0 = " << Vstat0_real << " + " << Vstat0_imag << " *1j" << std::endl;
out << " V1 = " << Vstat1_real << " + " << Vstat1_imag << " *1j" << std::endl;
}
out << std::endl;
}
Expand Down
65 changes: 43 additions & 22 deletions src/seqs_maxwell_frequency.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,8 +72,11 @@ SeqsMaxwellFrequencySolver::SeqsMaxwellFrequencySolver(MPI_Session &mpi, Electro
psi_real_ = NULL;
psi_imag_ = NULL;

V0_ = NULL;
V1_ = NULL;
V0_real_ = NULL;
V1_real_ = NULL;

V0_imag_ = NULL;
V1_imag_ = NULL;

phi_tot_real_ = NULL;
phi_tot_imag_ = NULL;
Expand All @@ -100,8 +103,10 @@ SeqsMaxwellFrequencySolver::~SeqsMaxwellFrequencySolver() {
delete A_real_;
delete phi_tot_real_;
delete phi_tot_imag_;
delete V1_;
delete V0_;
delete V1_imag_;
delete V0_imag_;
delete V1_real_;
delete V0_real_;
delete psi_imag_;
delete psi_real_;
delete phi_imag_;
Expand Down Expand Up @@ -299,13 +304,21 @@ void SeqsMaxwellFrequencySolver::SolveSEQS() {

ConstantCoefficient one(1.0);

V0_ = new ParGridFunction(pspace_);
*V0_ = 0.0;
V0_->ProjectBdrCoefficient(one, port_0_);
V0_real_ = new ParGridFunction(pspace_);
*V0_real_ = 0.0;
V0_real_->ProjectBdrCoefficient(one, port_0_);

V0_imag_ = new ParGridFunction(pspace_);
*V0_imag_ = 0.0;
V0_imag_->ProjectBdrCoefficient(one, port_0_);

V1_real_ = new ParGridFunction(pspace_);
*V1_real_ = 0.0;
V1_real_->ProjectBdrCoefficient(one, port_1_);

V1_ = new ParGridFunction(pspace_);
*V1_ = 0.0;
V1_->ProjectBdrCoefficient(one, port_1_);
V1_imag_ = new ParGridFunction(pspace_);
*V1_imag_ = 0.0;
V1_imag_->ProjectBdrCoefficient(one, port_1_);

if (verbose) grvy_printf(ginfo, "... Initializing right hand side.\n");
ParLinearForm *bp_real = new ParLinearForm(pspace_);
Expand All @@ -323,24 +336,30 @@ void SeqsMaxwellFrequencySolver::SolveSEQS() {
Kpp_real->AddDomainIntegrator(new DiffusionIntegrator(*rel_sig_));
Kpp_real->Assemble();
Kpp_real->Finalize();
Kpp_real->AddMult(*V0_, *bp_real, -em_opt_.Vstat0_real);
Kpp_real->AddMult(*V1_, *bp_real, -em_opt_.Vstat1_real);
Kpp_real->AddMult(*V0_real_, *bp_real, -em_opt_.Vstat0_real);
Kpp_real->AddMult(*V1_real_, *bp_real, -em_opt_.Vstat1_real);
Kpp_real->AddMult(*V0_imag_, *bp_imag, -em_opt_.Vstat0_imag);
Kpp_real->AddMult(*V1_imag_, *bp_imag, -em_opt_.Vstat1_imag);
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, -em_opt_.Vstat0_real);
Kpp_imag->AddMult(*V1_, *bp_imag, -em_opt_.Vstat1_real);
Kpp_imag->AddMult(*V0_real_, *bp_imag, -em_opt_.Vstat0_real);
Kpp_imag->AddMult(*V1_real_, *bp_imag, -em_opt_.Vstat1_real);
Kpp_imag->AddMult(*V0_imag_, *bp_real, em_opt_.Vstat0_imag);
Kpp_imag->AddMult(*V1_imag_, *bp_real, em_opt_.Vstat1_imag);
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, -em_opt_.Vstat0_real);
Kss_real->AddMult(*V1_, *bs_real, -em_opt_.Vstat1_real);
Kss_real->AddMult(*V0_real_, *bs_real, -em_opt_.Vstat0_real);
Kss_real->AddMult(*V1_real_, *bs_real, -em_opt_.Vstat1_real);
Kss_real->AddMult(*V0_imag_, *bs_imag, -em_opt_.Vstat0_imag);
Kss_real->AddMult(*V1_imag_, *bs_imag, -em_opt_.Vstat1_imag);
delete Kss_real;
}

Expand Down Expand Up @@ -578,13 +597,15 @@ void SeqsMaxwellFrequencySolver::SolveSEQS() {
// potential solve
phi_tot_real_ = new ParGridFunction(pspace_);
*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_);
add(*phi_tot_real_ , em_opt_.Vstat0_real, *V0_real_, *phi_tot_real_);
add(*phi_tot_real_ , em_opt_.Vstat1_real, *V1_real_, *phi_tot_real_);
*phi_tot_real_ += *phi_real_;
*phi_tot_real_ += *psi_real_;

phi_tot_imag_ = new ParGridFunction(pspace_);
*phi_tot_imag_ = 0.0;
add(*phi_tot_imag_ , em_opt_.Vstat0_imag, *V0_imag_, *phi_tot_imag_);
add(*phi_tot_imag_ , em_opt_.Vstat1_imag, *V1_imag_, *phi_tot_imag_);
*phi_tot_imag_ += *phi_imag_;
*phi_tot_imag_ += *psi_imag_;

Expand Down Expand Up @@ -937,8 +958,8 @@ void SeqsMaxwellFrequencySolver::EvaluateCurrent() {
delete Kpp_imag;

// call operator() for linear forms
const double I_stat_real = (*bp_real)(*V0_);
const double I_stat_imag = (*bp_imag)(*V0_);
const double I_stat_real = (*bp_real)(*V0_real_);
const double I_stat_imag = (*bp_imag)(*V0_real_);

if (verbose) grvy_printf(ginfo, "I_stat_real = %.8e\n", I_stat_real);
if (verbose) grvy_printf(ginfo, "I_stat_imag = %.8e\n", I_stat_imag);
Expand Down Expand Up @@ -985,8 +1006,8 @@ void SeqsMaxwellFrequencySolver::EvaluateCurrent() {
Knn_imag->AddMult(*n_imag_, *bp_real, -1.0);
delete Knn_imag;

const double I_ind_real = (*bp_real)(*V0_);
const double I_ind_imag = (*bp_imag)(*V0_);
const double I_ind_real = (*bp_real)(*V0_real_);
const double I_ind_imag = (*bp_imag)(*V0_real_);

delete bp_imag;
delete bp_real;
Expand Down
7 changes: 5 additions & 2 deletions src/seqs_maxwell_frequency.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,8 +89,11 @@ class SeqsMaxwellFrequencySolver {
mfem::ParGridFunction *psi_real_;
mfem::ParGridFunction *psi_imag_;

mfem::ParGridFunction *V0_;
mfem::ParGridFunction *V1_;
mfem::ParGridFunction *V0_real_;
mfem::ParGridFunction *V1_real_;

mfem::ParGridFunction *V0_imag_;
mfem::ParGridFunction *V1_imag_;

mfem::ParGridFunction *phi_tot_real_;
mfem::ParGridFunction *phi_tot_imag_;
Expand Down
14 changes: 11 additions & 3 deletions test/seqs.bar.test
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,19 +7,27 @@ setup() {
MESH_FILE=meshes/simple_bar.msh
CMD_LINE="--em-only -seqs -m $MESH_FILE --rtol 1e-12 --atol 1e-15 --maxiter 200"
CMD_LINE="$CMD_LINE -cd 1 -nbc 4 -p0 3 -p1 5"
CMD_LINE="$CMD_LINE -s 1e8 -f 0.0001 -Vs0r 1.0 -Vs1r -1.0"
CMD_LINE_REAL="$CMD_LINE -s 1e8 -f 0.0001 -Vs0r 1.0 -Vs1r -1.0"
CMD_LINE_IMAG="$CMD_LINE -s 1e8 -f 0.0001 -Vs0r 0.0 -Vs0i 1.0 -Vs1r 0.0 -Vs1i -1.0"
}

@test "[$TEST] verify tps SEQS Maxwell solver runs with np = 1" {
test -s $MESH_FILE
run mpirun -np 1 ../src/tps $CMD_LINE
run mpirun -np 1 ../src/tps $CMD_LINE_REAL
[[ "${status}" -eq 0 ]]
[[ "${output}" =~ "I_tot_real = 1.54692270e-02" ]]
}

@test "[$TEST] verify tps SEQS Maxwell solver runs with np = 4" {
test -s $MESH_FILE
run mpirun -np 4 ../src/tps $CMD_LINE
run mpirun -np 4 ../src/tps $CMD_LINE_REAL
[[ "${status}" -eq 0 ]]
[[ "${output}" =~ "I_tot_real = 1.54692270e-02" ]]
}

@test "[$TEST] verify tps SEQS Maxwell solver runs with imaginary drive voltage" {
test -s $MESH_FILE
run mpirun -np 4 ../src/tps $CMD_LINE_IMAG
[[ "${status}" -eq 0 ]]
[[ "${output}" =~ "I_tot_imag = 1.54692270e-02" ]]
}

0 comments on commit 4907b2d

Please sign in to comment.