beuler solver timestep and Jacobian calculation #3107
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change | ||||
|---|---|---|---|---|---|---|
|
|
@@ -123,6 +123,15 @@ SNESSolver::SNESSolver(Options* opts) | |||||
| upper_its((*options)["upper_its"] | ||||||
| .doc("Iterations above which the next timestep is reduced") | ||||||
| .withDefault(static_cast<int>(maxits * 0.8))), | ||||||
| timestep_factor_on_failure((*options)["timestep_factor_on_failure"] | ||||||
| .doc("Multiply timestep on convergence failure") | ||||||
| .withDefault(0.5)), | ||||||
| timestep_factor_on_upper_its((*options)["timestep_factor_on_upper_its"] | ||||||
| .doc("Multiply timestep if iterations exceed upper_its") | ||||||
| .withDefault(0.9)), | ||||||
| timestep_factor_on_lower_its((*options)["timestep_factor_on_lower_its"] | ||||||
| .doc("Multiply timestep if iterations are below lower_its") | ||||||
| .withDefault(1.4)), | ||||||
| diagnose( | ||||||
| (*options)["diagnose"].doc("Print additional diagnostics").withDefault(false)), | ||||||
| diagnose_failures((*options)["diagnose_failures"] | ||||||
|
|
@@ -297,13 +306,21 @@ int SNESSolver::init() { | |||||
| SNESSetJacobian(snes, Jmf, Jmf, MatMFFDComputeJacobian, this); | ||||||
|
|
||||||
| } else { | ||||||
| // Calculate the Jacobian using finite differences. The finite | ||||||
| // difference Jacobian (Jfd) may be used for both operator and | ||||||
| // preconditioner or, if matrix_free_operator, in only the | ||||||
| // preconditioner. | ||||||
|
|
||||||
| // Create a vector to store interpolated output solution | ||||||
| // Used so that the timestep does not have to be adjusted, | ||||||
| // because that would require updating the preconditioner. | ||||||
| ierr = VecDuplicate(snes_x, &output_x); | ||||||
| CHKERRQ(ierr); | ||||||
|
|
||||||
| if (use_coloring) { | ||||||
| // Use matrix coloring. | ||||||
| // This greatly reduces the number of times the rhs() function | ||||||
| // needs to be evaluated when calculating the Jacobian. | ||||||
|
|
||||||
| // Use global mesh for now | ||||||
| Mesh* mesh = bout::globals::mesh; | ||||||
|
|
@@ -717,14 +734,18 @@ int SNESSolver::run() { | |||||
| CHKERRQ(ierr); | ||||||
| } | ||||||
|
|
||||||
| BoutReal target = simtime; | ||||||
| for (int s = 0; s < getNumberOutputSteps(); s++) { | ||||||
| target += getOutputTimestep(); | ||||||
|
|
||||||
| bool looping = true; | ||||||
| int snes_failures = 0; // Count SNES convergence failures | ||||||
| int steps_since_snes_failure = 0; | ||||||
| int saved_jacobian_lag = 0; | ||||||
| int loop_count = 0; | ||||||
| do { | ||||||
| if (simtime >= target) | ||||||
| break; // Could happen if step over multiple outputs | ||||||
| if (scale_vars) { | ||||||
| // Individual variable scaling | ||||||
| // Note: If variables are rescaled then the Jacobian columns | ||||||
|
|
@@ -802,9 +823,15 @@ int SNESSolver::run() { | |||||
| dt = timestep; | ||||||
| looping = true; | ||||||
| if (simtime + dt >= target) { | ||||||
| // Note: When the timestep is changed the preconditioner needs to be updated | ||||||
| // => Step over the output time and interpolate if not matrix free | ||||||
|
|
||||||
| if (matrix_free) { | ||||||
| // Ensure that the timestep goes to the next output time and then stops. | ||||||
| // This avoids the need to interpolate | ||||||
| dt = target - simtime; | ||||||
| } | ||||||
| looping = false; | ||||||
| } | ||||||
|
|
||||||
| if (predictor and (time1 > 0.0)) { | ||||||
|
|
@@ -851,9 +878,10 @@ int SNESSolver::run() { | |||||
| } | ||||||
|
|
||||||
| ++snes_failures; | ||||||
| steps_since_snes_failure = 0; | ||||||
|
|
||||||
| // Try a smaller timestep | ||||||
| timestep *= timestep_factor_on_failure; | ||||||
| // Restore state | ||||||
| VecCopy(x0, snes_x); | ||||||
|
|
||||||
|
|
@@ -872,6 +900,7 @@ int SNESSolver::run() { | |||||
| jacobian_pruned = false; // Reset flag. Will be set after pruning. | ||||||
| } | ||||||
| if (saved_jacobian_lag == 0) { | ||||||
| // This triggers a Jacobian recalculation | ||||||
| SNESGetLagJacobian(snes, &saved_jacobian_lag); | ||||||
| SNESSetLagJacobian(snes, 1); | ||||||
| } | ||||||
|
|
@@ -943,6 +972,7 @@ int SNESSolver::run() { | |||||
| } | ||||||
|
|
||||||
| simtime += dt; | ||||||
| ++steps_since_snes_failure; | ||||||
|
|
||||||
| if (diagnose) { | ||||||
| // Gather and print diagnostic information | ||||||
|
|
@@ -952,7 +982,6 @@ int SNESSolver::run() { | |||||
| simtime, timestep, nl_its, lin_its, static_cast<int>(reason)); | ||||||
| if (snes_failures > 0) { | ||||||
| output.write(", SNES failures: {}", snes_failures); | ||||||
| } | ||||||
| output.write("\n"); | ||||||
| } | ||||||
|
|
@@ -1003,24 +1032,74 @@ int SNESSolver::run() { | |||||
| #endif // PETSC_VERSION_GE(3,20,0) | ||||||
|
|
||||||
| if (looping) { | ||||||
| // Consider changing the timestep. | ||||||
| // Note: The preconditioner depends on the timestep, | ||||||
| // so if it is not recalculated the it will be less | ||||||
| // effective. | ||||||
| if ((nl_its <= lower_its) && (timestep < max_timestep) | ||||||
| && (steps_since_snes_failure > 2)) { | ||||||
| // Increase timestep slightly | ||||||
| timestep *= timestep_factor_on_lower_its; | ||||||
|
|
||||||
| if (timestep > max_timestep) { | ||||||
| timestep = max_timestep; | ||||||
| } | ||||||
|
|
||||||
| // Note: Setting the SNESJacobianFn to NULL retains | ||||||
| // previously set evaluation function. | ||||||
| // | ||||||
| // The SNES Jacobian is a combination of the RHS Jacobian | ||||||
| // and a factor involving the timestep. | ||||||
| // Depends on equation_form | ||||||
| // -> Probably call SNESSetJacobian(snes, Jfd, Jfd, NULL, fdcoloring); | ||||||
|
|
||||||
| if (static_cast<BoutReal>(lin_its) / nl_its > 4) { | ||||||
| // Recompute Jacobian (for now) | ||||||
| if (saved_jacobian_lag == 0) { | ||||||
| SNESGetLagJacobian(snes, &saved_jacobian_lag); | ||||||
| SNESSetLagJacobian(snes, 1); | ||||||
| } | ||||||
| } | ||||||
|
|
||||||
| } else if (nl_its >= upper_its) { | ||||||
| // Reduce timestep slightly | ||||||
| timestep *= timestep_factor_on_upper_its; | ||||||
|
|
||||||
| // Recompute Jacobian | ||||||
| if (saved_jacobian_lag == 0) { | ||||||
| SNESGetLagJacobian(snes, &saved_jacobian_lag); | ||||||
| SNESSetLagJacobian(snes, 1); | ||||||
| } | ||||||
| } | ||||||
| } | ||||||
| snes_failures = 0; | ||||||
| } while (looping); | ||||||
|
|
||||||
| if (!matrix_free) { | ||||||
| ASSERT2(simtime >= target); | ||||||
| ASSERT2(simtime - dt < target); | ||||||
| // Stepped over output timestep => Interpolate | ||||||
| // snes_x is the solution at t = simtime | ||||||
| // x0 is the solution at t = simtime - dt | ||||||
| // Calculate output_x at t = target | ||||||
| VecCopy(snes_x, output_x); | ||||||
|
|
||||||
| // Note: If simtime = target then alpha = 0 | ||||||
| // and output_x = snes_x | ||||||
| BoutReal alpha = (simtime - target) / dt; | ||||||
|
Contributor
There was a problem hiding this comment. warning: variable 'alpha' of type 'BoutReal' (aka 'double') can be declared 'const' [misc-const-correctness]
Suggested change
|
||||||
|
|
||||||
| // output_x <- alpha * x0 + (1 - alpha) * output_x | ||||||
| VecAXPBY(output_x, alpha, 1. - alpha, x0); | ||||||
|
|
||||||
| } else { | ||||||
| // Timestep was adjusted to hit target output time | ||||||
| output_x = snes_x; | ||||||
| } | ||||||
|
|
||||||
| // Put the result into variables | ||||||
| if (scale_vars) { | ||||||
| // scaled_x <- output_x * var_scaling_factors | ||||||
| int ierr = VecPointwiseMult(scaled_x, output_x, var_scaling_factors); | ||||||
| CHKERRQ(ierr); | ||||||
|
|
||||||
| const BoutReal* xdata = nullptr; | ||||||
|
|
@@ -1031,15 +1110,15 @@ int SNESSolver::run() { | |||||
| CHKERRQ(ierr); | ||||||
| } else { | ||||||
| const BoutReal* xdata = nullptr; | ||||||
| int ierr = VecGetArrayRead(output_x, &xdata); | ||||||
| CHKERRQ(ierr); | ||||||
| load_vars(const_cast<BoutReal*>(xdata)); | ||||||
| ierr = VecRestoreArrayRead(output_x, &xdata); | ||||||
| CHKERRQ(ierr); | ||||||
| } | ||||||
| run_rhs(target); // Run RHS to calculate auxilliary variables | ||||||
|
|
||||||
| if (call_monitors(target, s, getNumberOutputSteps()) != 0) { | ||||||
| break; // User signalled to quit | ||||||
| } | ||||||
| } | ||||||
|
|
||||||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
warning: statement should be inside braces [readability-braces-around-statements]