Skip to content

Commit

Permalink
WIP (#75): Add E, B, and I calculations
Browse files Browse the repository at this point in the history 
First, compute E and B, since these are the physically meaningful
quantities.  Specifically, as a postprocessing step, we compute

B = curl(A)
E = -grad(phi) - j*omega*A

where phi and A denote the "full" or "total" potentials, which are
split up by Ostrowski and Hiptmair.

Only E and B are written to the paraview output now.

Regarding the current, it is split into two parts, the "static" and
"induced" parts, as defined by Ostrowski and Hiptmair, and evaluated
accordingly.
  • Loading branch information
@trevilo
trevilo committed Oct 21, 2021
1 parent b31ec76 commit 7ff67db
Show file tree
Hide file tree
Showing 2 changed files with 188 additions and 5 deletions.
182 changes: 177 additions & 5 deletions src/seqs_maxwell_frequency.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,9 +82,18 @@ SeqsMaxwellFrequencySolver::SeqsMaxwellFrequencySolver(MPI_Session &mpi, Electro
A_imag_ = NULL;
n_real_ = NULL;
n_imag_ = NULL;

B_real_ = NULL;
B_imag_ = NULL;
E_real_ = NULL;
E_imag_ = NULL;
}

SeqsMaxwellFrequencySolver::~SeqsMaxwellFrequencySolver() {
delete E_imag_;
delete E_real_;
delete B_imag_;
delete B_real_;
delete n_imag_;
delete n_real_;
delete A_imag_;
Expand Down Expand Up @@ -268,6 +277,8 @@ void SeqsMaxwellFrequencySolver::Initialize() {
void SeqsMaxwellFrequencySolver::Solve() {
SolveSEQS();
SolveStabMaxwell();
EvaluateEandB();
EvaluateCurrent();
WriteParaview();
}

Expand Down Expand Up @@ -787,7 +798,7 @@ void SeqsMaxwellFrequencySolver::SolveStabMaxwell() {
BDP.SetDiagonalBlock(3, Pnni);
BDP.owns_blocks = 1;

// 13. Solve the thing (finally!)
// Solve the thing (finally!)
FGMRESSolver solver(MPI_COMM_WORLD);
solver.SetPreconditioner(BDP);
solver.SetOperator(*maxwellOp);
Expand All @@ -805,6 +816,163 @@ void SeqsMaxwellFrequencySolver::SolveStabMaxwell() {
}
}

void SeqsMaxwellFrequencySolver::EvaluateEandB() {
bool verbose = mpi_.Root();
if (verbose) grvy_printf(ginfo, "Evaluating the magnetic field.\n");

assert(A_real_ != NULL);
assert(A_imag_ != NULL);

// First, we evaluate the B field, which is the curl of the magnetic
// vector potential. The full magnetic vector potential is given by
// A + grad(n). Thus, B = curl(A + grad(n)) = curl(A).
B_real_ = new ParGridFunction(Bspace_);
*B_real_ = 0;

B_imag_ = new ParGridFunction(Bspace_);
*B_imag_ = 0;

ParDiscreteLinearOperator *curl = new ParDiscreteLinearOperator(Aspace_, Bspace_);
curl->AddDomainInterpolator(new CurlInterpolator);
curl->Assemble();
curl->Finalize();

curl->Mult(*A_real_, *B_real_);
curl->Mult(*A_imag_, *B_imag_);

delete curl;

if (verbose) grvy_printf(ginfo, "Evaluating the electric field.\n");
assert(phi_tot_real_ != NULL);
assert(phi_tot_imag_ != NULL);
assert(n_real_ != NULL);
assert(n_imag_ != NULL);

// Second, we evaluate the E field, which is the opposite of the
// gradient of the electric potential minus the time derivative of
// the magnetic vector potential:
//
// E = -grad(phi_tot) - dA_tot/dt,
//
// where phi_tot and A_tot are the "total" potentials (i.e. phi_tot
// = phi+psi+V_0+V_1) and (A_tot = A + grad(n)). In the frequency
// domain, we have
//
// E = -grad(phi_tot) - j*omega*A_tot.
//
// After non-dimensionalizing, this becomes
//
// E = -grad(phi_tot) - j*A_tot,
//
// such that
//
// E_real = -grad(phi_tot_real) + A_tot_imag
// E_imag = -grad(phi_tot_imag) - A_tot_real
E_real_ = new ParGridFunction(Aspace_);
*E_real_ = 0;

E_imag_ = new ParGridFunction(Aspace_);
*E_imag_ = 0;

ParDiscreteLinearOperator *grad = new ParDiscreteLinearOperator(pspace_, Aspace_);
grad->AddDomainInterpolator(new GradientInterpolator);
grad->Assemble();
grad->Finalize();

grad->AddMult(*phi_tot_real_, *E_real_, -1.0);
grad->AddMult(*phi_tot_imag_, *E_imag_, -1.0);

*E_real_ += *A_imag_;
*E_imag_ -= *A_real_;

grad->AddMult(*n_imag_, *E_real_, +1.0);
grad->AddMult(*n_real_, *E_imag_, -1.0);

delete grad;
}

void SeqsMaxwellFrequencySolver::EvaluateCurrent() {
bool verbose = mpi_.Root();
if (verbose) grvy_printf(ginfo, "Evaluating the current.\n");

ParLinearForm *bp_real = new ParLinearForm(pspace_);
*bp_real = 0.0;

ParLinearForm *bp_imag = new ParLinearForm(pspace_);
*bp_imag = 0.0;

// First, the "static" part
ParBilinearForm *Kpp_real = new ParBilinearForm(pspace_);
Kpp_real->AddDomainIntegrator(new DiffusionIntegrator(*rel_sig_));
Kpp_real->Assemble();
Kpp_real->Finalize();
Kpp_real->AddMult(*phi_tot_real_, *bp_real, 1.0);
Kpp_real->AddMult(*phi_tot_imag_, *bp_imag, 1.0);
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(*phi_tot_imag_, *bp_real, -1.0);
Kpp_imag->AddMult(*phi_tot_real_, *bp_imag, 1.0);
delete Kpp_imag;

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

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);

// these linear forms are re-used, so zero them out
*bp_real = 0.0;
*bp_imag = 0.0;


// Second, the "induced" part
ParMixedBilinearForm *Kna_imag = new ParMixedBilinearForm(Aspace_, pspace_);
Kna_imag->AddDomainIntegrator(new VectorFEWeakDivergenceIntegrator(*rel_sig_)); // should be -sigma
Kna_imag->Assemble();
Kna_imag->Finalize();

Kna_imag->AddMult(*A_real_, *bp_imag, -1.0); // sign flip accounts for -sigma (see above)
Kna_imag->AddMult(*A_imag_, *bp_real, 1.0); // sign flip accounts for -sigma (see above)

ParMixedBilinearForm *Kna_real = new ParMixedBilinearForm(Aspace_, pspace_);
Kna_real->AddDomainIntegrator(new VectorFEWeakDivergenceIntegrator(*one_over_sigma_)); // should be -1/sigma
Kna_real->Assemble();
Kna_real->Finalize();

Kna_real->AddMult(*A_real_, *bp_real, 1.0); // sign flip accounts for -1/sigma (see above)
Kna_real->AddMult(*A_imag_, *bp_imag, 1.0); // sign flip accounts for -1/sigma (see above)

ParBilinearForm *Knn_real = new ParBilinearForm(pspace_);
Knn_real->AddDomainIntegrator(new DiffusionIntegrator(*one_over_sigma_));
Knn_real->Assemble();
Knn_real->Finalize();

Knn_real->AddMult(*n_real_, *bp_real, -1.0);
Knn_real->AddMult(*n_imag_, *bp_imag, -1.0);

ParBilinearForm *Knn_imag = new ParBilinearForm(pspace_);
Knn_imag->AddDomainIntegrator(new DiffusionIntegrator(*rel_sig_));
Knn_imag->Assemble();
Knn_imag->Finalize();

Knn_imag->AddMult(*n_real_, *bp_imag, 1.0);
Knn_imag->AddMult(*n_imag_, *bp_real, -1.0);

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

if (verbose) grvy_printf(ginfo, "I_ind_real = %.8e\n", I_ind_real);
if (verbose) grvy_printf(ginfo, "I_ind_imag = %.8e\n", I_ind_imag);

if (verbose) grvy_printf(ginfo, "I_tot_real = %.8e\n", I_stat_real + I_ind_real);
if (verbose) grvy_printf(ginfo, "I_tot_imag = %.8e\n", I_stat_imag + I_ind_imag);
}

void SeqsMaxwellFrequencySolver::WriteParaview() {
bool verbose = mpi_.Root();
if (verbose) grvy_printf(ginfo, "Writing Maxwell solution to paraview output.\n");
Expand All @@ -815,9 +983,13 @@ void SeqsMaxwellFrequencySolver::WriteParaview() {
paraview_dc.SetDataFormat(VTKFormat::BINARY);
paraview_dc.SetHighOrderOutput(true);
paraview_dc.SetTime(0.0);
paraview_dc.RegisterField("phi_tot_real", phi_tot_real_);
paraview_dc.RegisterField("phi_tot_imag", phi_tot_imag_);
paraview_dc.RegisterField("A_real", A_real_);
paraview_dc.RegisterField("A_imag", A_imag_);
// paraview_dc.RegisterField("phi_tot_real", phi_tot_real_);
// paraview_dc.RegisterField("phi_tot_imag", phi_tot_imag_);
// paraview_dc.RegisterField("A_real", A_real_);
// paraview_dc.RegisterField("A_imag", A_imag_);
paraview_dc.RegisterField("E_real", E_real_);
paraview_dc.RegisterField("E_imag", E_imag_);
paraview_dc.RegisterField("B_real", B_real_);
paraview_dc.RegisterField("B_imag", B_imag_);
paraview_dc.Save();
}
11 changes: 11 additions & 0 deletions src/seqs_maxwell_frequency.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,12 +100,23 @@ class SeqsMaxwellFrequencySolver {
mfem::ParGridFunction *n_real_;
mfem::ParGridFunction *n_imag_;

mfem::ParGridFunction *B_real_;
mfem::ParGridFunction *B_imag_;
mfem::ParGridFunction *E_real_;
mfem::ParGridFunction *E_imag_;

// Solve for the electric potential
void SolveSEQS();

// Solve for the magnetic vector potential
void SolveStabMaxwell();

// Evaluate electric and magnetic fields
void EvaluateEandB();

// Evaluate the current
void EvaluateCurrent();

// Write paraview data
void WriteParaview();

Expand Down

0 comments on commit 7ff67db

Please sign in to comment.