A generalized weighted residual method for RFP plasma simulation. Jan Scheffel Fusion Plasma Physics Alfvén Laboratory, KTH Stockholm, Sweden. OUTLINE. • What is the GWRM? • ODE example • SIR - a globally convergent root solver • Accuracy • Efficiency • Discussion
Fusion Plasma Physics
Alfvén Laboratory, KTH
Time differencing numerical initial value schemes (even implicit)
require extremelymany time steps for problems of physical
interest, where there are several separated time scales.
Causality is already embedded in the governing PDE’s –
- there is no need to mimic causality by time stepping.
Spectral methods (solution expanded in basis functions) are
successful in the spatial domain – why not employ them also in
the time domain?
By expanding in time + physical space + physical parameters, the
computational result will be semi-analytical. (Analytic basis
functions with numerical coefficients).
Ideal for scaling studies, for example.
What is the simulation
Generalized Weighted Spectral Method (GWRM) ?
Fully spectral weighted residual method for semi-analytical solution of initial value partial differential equations.
The method is acausal and thus avoids time step limitations.
The spectral coefficients are determined by iterative solution of a nonlinear system of algebraic eqs, for which a globally convergent semi-implicit root solver (SIR) has been developed.
Consider a system of parabolic or hyperbolic initial-value PDE’s, symbolically written as
D is a nonlinear matrix operator, f is a forcing term.
D and f contains both physical variables and physical free parameters (denoted p).
Initial u(0,x;p) + (Dirichlet, Neumann or Robin) boundary u(t,xB;p) conditions.
Integrate in time:
Solution u(t,x;p) is approximated by finite, first kind Chebyshev polynomial series.
Definition: Chebyshev polynomial Tn(x) = cos(n arccosx).
For simplicity – here single equation, one spatial dimension x, one physical parameter p.
The Weighted Residual of the GWRM is given by
The TP-WRM coefficients are now obtained from the nonlinear system of algebraic equations
The initial state is expanded as
Solution to be determined:
The coefficients aqrs
by iterations, using
a root solver.
for 1≤ q ≤ K + 1
Boundary conditions enter here
Light a match – a model of flame propagation:
y – flame radius
green - exact solution
d = 0.05; # Chebyshev modes K = 20, # time domains Nt = 1, error = 0.01
The GWRM - well adapted for iterative methods for two reasons:
1) Basic Chebyshev coefficient equations are of the standard iterative form
2) Initial estimate of solution vector can be chosen sufficiently close to the solution
by reducing the solution time interval
Instead of using direct iteration,
the Semi-Implicit Root solver (SIR)
finds the roots to the equations
or, in matrix form
has the same solutions as the original system, but contains free parameters in the form
of the components of the matrix A. The parameters can be chosen to control
the gradients of the hypersurfaces . Adjusting these parameters, global, quasi-monotonous
and superlinear convergence is attained. In SIR,
are finite and is controlled; it produces limited step lengths, quasi-monotonous convergence;
and approaches zero after some initial iterations.
Newton’s method is a special case of the present method, when all
• Rapid second order convergence is generally approached after some iteration steps.
• Relationship to Newton’s method - approximately similar numerical work;
inversion of a Jacobian matrix at each iteration step.
Newton’s method simulation
Newton’s method with linesearch simulation
f ≠ 0
Solution compared to Lax-Wendroff (explicit time differencing):
GWRM parameters - (S,M,N) = (2,7,6), 13 iterations.
L-W marginally stable parameters -
TP-WRM solution, using
two spatial subdomains
TP-WRM solution error, as
compared to exact solution
Result: GWRM 50 % faster than L-W
for same accuracy.
GWRM Burger equation solution, simulation
including viscosity dependence u = u(t,x;v)
Initial conditionj(x) = x(1 - x) and boundary condition u(t,0;v) = u(t,1;v) = 0.
Solution shown versus x and v at time t = 2.5. Here K = 8, L = 10, and M = 2.
Wave equation, forced (hyperbolic)
fast time scale)
(slow + rapid time scale)
(K,L) = (6,8)
∆x = 1/30, 900 time steps
∆x = 1/30, 100 time steps
GWRM work so far:
• The time- and parameter-generalized weighted residual method,
J. Scheffel, 2008. (GWRM method outlined)
• Semi-analytical solution of initial-value problems,
D. Lundin, 2006. (Resistive MHD stability of RFP and z-pinch)
• Application of the time- and parameter generalized weighted reidual
method to systems of nonlinear equations,
D. Jackson, 2007. (Navier-Stokes equations, Rayleigh-Taylor instability)
• Further development and implementation of the GWRM
A. Mirza, ongoing Ph D studies (Application of GWRM to nonlinear resistive MHD)
• Solution of systems of nonlinear equations, a semi-implicit approach,
J. Scheffel, 2006. (SIR outlined)
• Studies of a semi-implicit root solver,
C. Håkansson, M Sc Thesis. (Efficient SIR compared to other methods)
• The GWRM is shown to be accurate for spatially smooth solutions - convergence including sharp gradients should be further studied.
• Efficiency is central; SIR involves Jacobian matrix inversion of Chebyshev coefficient eqs - N eqs takes O[N3] operations.
• Methods to improve efficiency have been developed - temporal subdomains and spatial subdomain techniques using overlapping domains.
• Further benchmarking of efficiency for MHD relevant test problems should be carried out as well as corresponding comparisons with implicit methods.
• SIR efficiency and linking to GWRM is presently being optimized.
for solution of initial value partial differential equations, has been outlined.
by Chebyshev series, semi-analytical solutions can be obtained as
(“Semi-analytical”: expansion in basis functions with numerical coefficients.)
Computed solutions thus contain r- t- and parameter dependence explicitly.
The method is global and avoids time step limitations.
system of algebraic equations, for which an efficient semi-implicit root solver
(SIR) has been developed.
To improve efficiency, a spatial subdomain approach has been developed.
Problems in fluid mechanics and MHD will be addressed.
plasma pressure in stochastic magnetic field geometries, in particular
operational limits in reversed-field pinches.