Skip to content

Adaptive time-stepping for transient solver using SUNDIALS #292

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 41 commits into from
Dec 3, 2024
Merged

Adaptive time-stepping for transient solver using SUNDIALS #292

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
41 commits
Select commit Hold shift + click to select a range
f0ef4f2
Add SUNDIALS dependency and build MFEM with SUNDIALS
simlapointe Oct 14, 2024
e2edda8
Add weak curl operator for first-order transient formulation
simlapointe Oct 14, 2024
0c4d8e1
Add new transient solver config options for adaptive time-stepping
simlapointe Oct 14, 2024
d437156
Add first-order transient formulation with SUNDIALS ARKode schemes
simlapointe Oct 14, 2024
0c10b67
Change first order ODE system formulation, add CVODE integrator, remo…
simlapointe Oct 21, 2024
983cc9c
Use explicit RHS (y' = M^-1 f(y,t)) for both ARKODE and CVODE
simlapointe Oct 25, 2024
cff74b9
Remove weakCurl operator, no longer needed
simlapointe Oct 30, 2024
e5b1e1d
Remove unnecessary ODE integrator prints
simlapointe Oct 30, 2024
2f10218
Merge branch 'main' into simlapointe/transient-adapt-dt
simlapointe Oct 30, 2024
c4b4a01
Remove 2nd order ODE formulation
simlapointe Nov 5, 2024
937fb47
Update reference for transient solver tests
simlapointe Nov 5, 2024
b9ce371
Clean up
simlapointe Nov 5, 2024
11092e9
Add SUNDIALS to spack install instructions
simlapointe Nov 5, 2024
11ad777
Change SUNDIALS version for spack installer
simlapointe Nov 5, 2024
aea8d7f
Fix formatting issues
simlapointe Nov 5, 2024
1bb2a95
Merge branch 'main' into simlapointe/transient-adapt-dt
simlapointe Nov 6, 2024
cd7cfc4
Use MFEM's SDIRK23 solver as the fixed step Runge-Kutta option, along…
simlapointe Nov 7, 2024
c86a4c6
Remove SUNDIALS MAGMA dependency
simlapointe Nov 8, 2024
100d2f4
Fix formatting issues
simlapointe Nov 8, 2024
8ed5f1e
Change Vector management to avoid GPU memory issue
simlapointe Nov 11, 2024
f1ebc29
Fix formatting issues
simlapointe Nov 11, 2024
cc65463
Fix typo
simlapointe Nov 11, 2024
2c59a02
Make rhs/RHS lower/upper case consistent
simlapointe Nov 13, 2024
47a04e9
Use MakeRef instead of Get/SetSubVector to split ODE system vectors
simlapointe Nov 13, 2024
137ec8b
Fix formatting issues
simlapointe Nov 13, 2024
df38dbb
Simplify transient solver config file verification
simlapointe Nov 13, 2024
1fd974c
Merge branch 'main' into simlapointe/transient-adapt-dt
simlapointe Nov 13, 2024
9ecd6d1
fix formatting issues
simlapointe Nov 13, 2024
db8c2a5
Add dB/dt to ODE system instead of using trapezoidal integration
simlapointe Nov 19, 2024
1b0b6b2
Remove unused En vector
simlapointe Nov 19, 2024
56844da
Fix formatting issues
simlapointe Nov 19, 2024
7469ea9
Merge branch 'main' into simlapointe/transient-adapt-dt
simlapointe Nov 19, 2024
d615337
Update coaxial example reference
simlapointe Nov 19, 2024
d36f954
Update time domain formulation reference doc
simlapointe Nov 19, 2024
443b3b8
Add SUNDIALS to configuration options doc
simlapointe Nov 19, 2024
863aee4
Fix formatting issues
simlapointe Nov 19, 2024
4eab1b5
Doc fixes
simlapointe Dec 2, 2024
e6034e5
Remove unnecessary sundials options
simlapointe Dec 2, 2024
abde1ec
Merge main
simlapointe Dec 2, 2024
2c2b9a8
Incorporate PR feedback
simlapointe Dec 2, 2024
d3d18d7
Fix formatting
simlapointe Dec 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
120 changes: 48 additions & 72 deletions palace/models/timeoperator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,9 +39,6 @@ class TimeDependentFirstOrderOperator : public mfem::TimeDependentOperator
mutable Vector RHS;
int size_E;

// Indices of the first and second halves of the solution vector, sol = [Edot, E]
mfem::Array<int> idx1, idx2;

// Bindings to SpaceOperator functions to get the system matrix and preconditioner, and
// construct the linear solver.
std::function<void(double dt)> ConfigureLinearSolver;
Expand All @@ -56,15 +53,6 @@ class TimeDependentFirstOrderOperator : public mfem::TimeDependentOperator
// Get dimensions of E and Edot vectors.
size_E = space_op.GetNDSpace().GetTrueVSize();

// Set array indices
idx1.SetSize(size_E);
idx2.SetSize(size_E);
for (int i = 0; i < size_E; i++)
{
idx1[i] = i;
idx2[i] = i + size_E;
}

// Construct the system matrices defining the linear operator. PEC boundaries are
// handled simply by setting diagonal entries of the mass matrix for the corresponding
// dofs. Because the Dirichlet BC is always homogeneous, no special elimination is
Expand Down Expand Up @@ -120,10 +108,12 @@ class TimeDependentFirstOrderOperator : public mfem::TimeDependentOperator
u2.UseDevice(true);
rhs1.UseDevice(true);
rhs2.UseDevice(true);
u.GetSubVector(idx1, u1);
u.GetSubVector(idx2, u2);
rhs.GetSubVector(idx1, rhs1);
rhs.GetSubVector(idx2, rhs2);
u.Read();
u1.MakeRef(const_cast<Vector &>(u), 0, size_E);
u2.MakeRef(const_cast<Vector &>(u), size_E, size_E);
rhs.ReadWrite();
rhs1.MakeRef(rhs, 0, size_E);
rhs2.MakeRef(rhs, size_E, size_E);

// u1 = Edot, u2 = E
// rhs1 = -(K * u2 + C * u1) - J(t)
Expand All @@ -136,9 +126,6 @@ class TimeDependentFirstOrderOperator : public mfem::TimeDependentOperator
linalg::AXPBYPCZ(-1.0, rhs1, dJ_coef(t), NegJ, 0.0, rhs1);

rhs2 = u1;

rhs.SetSubVector(idx1, rhs1);
rhs.SetSubVector(idx2, rhs2);
}

// Solve M du = rhs
Expand All @@ -153,21 +140,20 @@ class TimeDependentFirstOrderOperator : public mfem::TimeDependentOperator
}
FormRHS(u, RHS);

Vector du1, du2, rhs1, rhs2;
Vector du1, du2, RHS1, RHS2;
du1.UseDevice(true);
du2.UseDevice(true);
rhs1.UseDevice(true);
rhs2.UseDevice(true);
du.GetSubVector(idx1, du1);
du.GetSubVector(idx2, du2);
RHS.GetSubVector(idx1, rhs1);
RHS.GetSubVector(idx2, rhs2);

kspM->Mult(rhs1, du1);
du2 = rhs2;

du.SetSubVector(idx1, du1);
du.SetSubVector(idx2, du2);
RHS1.UseDevice(true);
RHS2.UseDevice(true);
du.ReadWrite();
du1.MakeRef(du, 0, size_E);
du2.MakeRef(du, size_E, size_E);
RHS.ReadWrite();
RHS1.MakeRef(RHS, 0, size_E);
RHS2.MakeRef(RHS, size_E, size_E);

kspM->Mult(RHS1, du1);
du2 = RHS2;
}

void ImplicitSolve(double dt, const Vector &u, Vector &k) override
Expand All @@ -185,25 +171,24 @@ class TimeDependentFirstOrderOperator : public mfem::TimeDependentOperator
Mpi::Print("\n");
FormRHS(u, RHS);

Vector k1, k2, rhs1, rhs2;
Vector k1, k2, RHS1, RHS2;
k1.UseDevice(true);
k2.UseDevice(true);
rhs1.UseDevice(true);
rhs2.UseDevice(true);
k.GetSubVector(idx1, k1);
k.GetSubVector(idx2, k2);
RHS.GetSubVector(idx1, rhs1);
RHS.GetSubVector(idx2, rhs2);

// A k1 = rhs1 - dt K rhs2
K->AddMult(rhs2, rhs1, -dt);
kspA->Mult(rhs1, k1);
RHS1.UseDevice(true);
RHS2.UseDevice(true);
k.ReadWrite();
k1.MakeRef(k, 0, size_E);
k2.MakeRef(k, size_E, size_E);
RHS.ReadWrite();
RHS1.MakeRef(RHS, 0, size_E);
RHS2.MakeRef(RHS, size_E, size_E);

// A k1 = RHS1 - dt K RHS2
K->AddMult(RHS2, RHS1, -dt);
kspA->Mult(RHS1, k1);

// k2 = rhs2 + dt k1
linalg::AXPBYPCZ(1.0, rhs2, dt, k1, 0.0, k2);

k.SetSubVector(idx1, k1);
k.SetSubVector(idx2, k2);
linalg::AXPBYPCZ(1.0, RHS2, dt, k1, 0.0, k2);
}

void ExplicitMult(const Vector &u, Vector &v) const override { Mult(u, v); }
Expand All @@ -230,28 +215,29 @@ class TimeDependentFirstOrderOperator : public mfem::TimeDependentOperator
// Solve (Mass - dt Jacobian) x = Mass b
int SUNImplicitSolve(const Vector &b, Vector &x, double tol) override
{
Vector b1, b2, x1, x2, rhs;
Vector b1, b2, x1, x2, RHS1;
b1.UseDevice(true);
b2.UseDevice(true);
x1.UseDevice(true);
x2.UseDevice(true);
rhs.UseDevice(true);
b.GetSubVector(idx1, b1);
b.GetSubVector(idx2, b2);
x.GetSubVector(idx1, x1);
x.GetSubVector(idx2, x2);
RHS.GetSubVector(idx1, rhs);
RHS1.UseDevice(true);
b.Read();
b1.MakeRef(const_cast<Vector &>(b), 0, size_E);
b2.MakeRef(const_cast<Vector &>(b), size_E, size_E);
x.ReadWrite();
x1.MakeRef(x, 0, size_E);
x2.MakeRef(x, size_E, size_E);
RHS.ReadWrite();
RHS1.MakeRef(RHS, 0, size_E);

// A x1 = M b1 - dt K b2
M->Mult(b1, rhs);
K->AddMult(b2, rhs, -saved_gamma);
kspA->Mult(rhs, x1);
M->Mult(b1, RHS1);
K->AddMult(b2, RHS1, -saved_gamma);
kspA->Mult(RHS1, x1);

// x2 = b2 + dt x1
linalg::AXPBYPCZ(1.0, b2, saved_gamma, x1, 0.0, x2);

x.SetSubVector(idx1, x1);
x.SetSubVector(idx2, x2);
return 0;
}
};
Expand All @@ -270,24 +256,16 @@ TimeOperator::TimeOperator(const IoData &iodata, SpaceOperator &space_op,
int size_E = space_op.GetNDSpace().GetTrueVSize();
int size_B = space_op.GetRTSpace().GetTrueVSize();

// Set indices of Edot and E in the solution vector, sol = [Edot, E]
idx1.SetSize(size_E);
idx2.SetSize(size_E);
for (int i = 0; i < size_E; i++)
{
idx1[i] = i;
idx2[i] = i + size_E;
}

// Allocate space for solution vectors.
sol.SetSize(2 * size_E);
E.SetSize(size_E);
En.SetSize(size_E);
B.SetSize(size_B);
sol.UseDevice(true);
E.UseDevice(true);
En.UseDevice(true);
B.UseDevice(true);
sol.ReadWrite();
E.MakeRef(sol, size_E, size_E);

// Create ODE solver for 1st-order IVP.
mfem::TimeDependentOperator::Type type = mfem::TimeDependentOperator::IMPLICIT;
Expand Down Expand Up @@ -401,7 +379,6 @@ void TimeOperator::Init()
{
// Always use zero initial conditions.
sol = 0.0;
sol.GetSubVector(idx2, E);
B = 0.0;
if (use_mfem_integrator)
{
Expand All @@ -411,10 +388,9 @@ void TimeOperator::Init()

void TimeOperator::Step(double &t, double &dt)
{
sol.GetSubVector(idx2, En);
En = E;
double dt_input = dt;
ode->Step(sol, t, dt);
sol.GetSubVector(idx2, E);
// Ensure user-specified dt does not change.
dt = dt_input;

Expand Down
3 changes: 0 additions & 3 deletions palace/models/timeoperator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,9 +37,6 @@ class TimeOperator
double rel_tol, abs_tol;
bool use_mfem_integrator = false;

// Indices of the first and second halves of the solution vector, sol = [Edot, E]
mfem::Array<int> idx1, idx2;

// Discrete curl for B-field time integration (not owned).
const Operator *Curl;

Expand Down
Loading