From bromley@physics.utah.edu Tue Nov 2 09:41:04 2004 Date: Fri, 22 Oct 2004 18:07:10 -0600 (MDT) From: Benjamin Bromley To: rprice@phys.utb.edu, lburko@a5198.bates.edu, friedman@thales.phys.uwm.edu, jsread@uwm.edu, uryu@sissa.it, bromley@physics.utah.edu, burko@physics.utah.edu, cbeetle@physics.fau.edu, lau@phys.utb.edu, mc@physics.utah.edu, michael@phys.utb.edu, owen@tapir.caltech.edu Subject: PSW Report Hi, Here is progress report N%P+1. I hope all is going well for everyone. best wishes, Ben PSW Progress Report B. C. Bromley 2004-October-22 The 3-d scalar model: an eigen-approach executive summary. We have a 3-d finite difference code to solve the problem for the linear scalar field produced by two point sources rotating at fixed separation, a. Our coordinates are simple spherical polars, chosen so that the polar axis is aligned with the sources to mimic the behavior of adapted (Roswell) coordinates at large distance r>>a. In a previous report (2004-09-24), we noted that solutions obtained with this code tend to be noisy, perhaps as a result of grid geometry. The noise evidently appears to be high-frequency, carried by high-order spherical harmonic modes. To cure this problem, we found eigenvectors of the {\em finite difference form} of the angular part of the Laplace operator. These eigens are the grid-based analogs of the usual spherical harmonics, Y_lm. We modified the finite-difference solver by projecting some of the finite-difference equations (involving distances r>>a, near the outer boundary) onto these eigens. We were then able to remove the high-frequency noise by explicity setting high-order "eigenmodes" to zero. Preliminary results are very promising, although some kinks still need to be ironed out. I. Background. Toward solving the binary black hole problem in the periodic standing wave approximation, we have worked on much simpler toy models, including the 3-D linear, scalar field of a rotating pair of sources. Numerical solution of this toy is straightfoward when we use a finite-difference approach in spherical coordinates aligned with the rotation axis, as in Andrade et al. 2004 (Paper I). However, the finite-difference method has limitations, mostly in that it requires the inversion of a fairly large matrix. In Paper I, we used parallel processing and 16 GB of shared memory to find the field on a modest grid, with n_r x n_theta x n_phi equal to roughly 100 x 20 x 20. That was our limit, given the hardware available to us at the time, and the software installed on the machine. In last week's report (2004-10-15), Richard suggested that we continue to use a finite difference method, but with more sophisticated matrix-inversion software and bigger computers. We've tried a variety of the off-the-sheld iterative sparse solvers, including Krylov methods (e.g. the biconjugate gradient), but none seem to converge for our linear, rotating scalar model, although the reason is not well-understood (at least not by this reporter). However, there are direct solvers available which have yet to be tried (e.g., scalapack), and thanks to the LSU group's offer of computing resources, we will be able to reassess our limitations with finite difference solvers in spherical coordinates. Meanwhile, we are also pursuing alternatives, including a finite difference method in adapted (Roswell) coordinates and spectral methods. The UWM group has implemented a spherical harmonic code that works very well for linear and nonlinear cases with modest values of lambda (|lambda| < 1; pls see Paper I and the recent UWM report). Here, we report some progress on the use of adapted coordinates. Although we concentrate on finite differences with spherical polar coordinates aligned with the sources, these coordinates are identical to Roswells in the limit large distance r from the sources. We believe the major difficulties we have encountered with finite-diff+Roswells lie in this large-distance regime. In our last report on this topic we discussed what the problem was: spurious solutions to the FDE which amount to high-frequency "noise" in our computed solution as compared with the exact solution. Our best guess as to the origin of this noise is a poor match between the symmetry of the solution and the symmetry of our finite difference grid. In this report, we argue that we have found a way to control this noise, at least the high-frequency part. II. Method. Particulars of our finite-difference approach are as follows: * our model is (in Paper I terminology) 3-D scalar, Omega=0.3, lambda=0. * coordinates are "source aligned" spherical polars (z-axis || source axis). * inner boundary is Dirichlet and is located outside of the source's orbit. The sources do not lie in the computational domain. Our goal is to suppress the noise which appears when we solve the finite-difference equations in these source-aligned coordinates. To do this we break up the finite-difference form of the Laplace operator L, a matrix, into radial (L_r) and angular (L_ang) parts. We find eigenvectors of the L_ang part, call them q_i, with eigenvalues w_i (pls see footnote for details). These are similar to the spherical harmonics Y_lm, and a quick inspection of the w_i and the form of the q_i (functions of theta and phi grid locations) can tell us how to relate the q_i to the Y_lm. Of course, in the limit of small grid spacing, the q_i become the Y_lm. Next we rewrite the finite-difference equations in terms of coefficients of the field projected on to the basis set q_i. In other words, we are making a simple linear transformation of the finite-difference equations into equations explicitly involving the individual eigenvectors. Some of these equations can be modified to so that the amplitude of mode q_i in the solution is set to a particular value. Since the q_i have a one-to-one correspondence with the Y_lm's (the eigens become the spherical Y_lm's in the limit of small grid spacing), we can be guided by our intuition based on the true spherical harmonics to determine how to modify the original finite difference equations. In our scalar 3-D toy model, we know that high-order Y_lm's are not important at large r, so we set the mode amplitudes of the corresponding "high-order" q_i to be zero as well. This has the benefit of allowing us to stomp on the high-frequency "noise". III. Results and discussion. The procedure described above seems to work well. The first test of the eigen-code is to see if it gives the same result as the original "pure" finite-difference method. Indeed, if we keep all eigenvectors in our transformation of the finite difference equations into the q_i eigenbasis, we recover the original result. And, when we "zero out" high-order q_i (corresponding to hig-order Y_lm) we get a beautifully smooth result, significantly better than the pure FDE on comparison with the exact solution. The bad news is that we are picking up an unwanted standing wave solution. The good news is that the original finite difference code also sees to have this problem. Thus the new method does exactly what it is supposed to do: It gets rid of noise, and in this case, enough noise to see what the problem is. We think the remaining difficulties lie in a coding error in the way boundaries are handled or the way the Omega terms in the field equation are treated (see Paper I). Nonetheless, the fact that we can replace a fraction of finite difference equations with equations involving our new eigenmodes means that we have successfully merged a pure finite-difference code with a "spectral" method. If the outcome of this procedure is favorable, it would just produce a very warm feeling that a pure spectral approach is the way to fly. **** Footnote: Actually, the q_i are not quite eigens of L_ang, but solve the generalized eigenvalue problem P L_ang q_i = w_i P q_i where P is a "measure matrix", the equivalent of "sin(theta)" in the orthogonality relation integral_{phi=0..2pi,theta=0..pi} sin theta dtheta dphi Y_lm Y_l'm' = Kroneker-delta_l,l';m,m' Richard wrote up a very nice set of notes in latex on this. %%%%%%%%%%%%% ADDENDUM FROM OCT 30 %%%%%%%%%%%%%%%% From bromley@physics.utah.edu Tue Nov 2 09:41:30 2004 Date: Sat, 30 Oct 2004 17:11:12 -0600 (MDT) From: Benjamin Bromley To: rprice@phys.utb.edu Subject: m ============================================================================ PSW supplemental report Finite-difference methods and source-aligned coordinates. executive summary. We tried a lot of stuff. Much of it was useless because of very sloppy coding. But now it works^*. ^* Some details will be discussed elsewhere. ============================================================================ The details. The code I've been working on is spherical polar, source aligned, Omega=0.3, a=1, lambda=0. It has three basic switches governing the method of solution of the field, psi, with outgoing boundary conditions: * method FD. Set up and solve the straightforward finite diff equations. Use source aligned coordinates throughout. * method EV. In the inner part of the computational domain, use the FD equations. At large radii (e.g., in the wave zone), transform the FD equations by projecting them onto a basis consisting of the eigenvectors of the angular part of the Laplace operator. In the continuum limit, this means doing a projection onto Y_lm. Throw away (that is, set to zero) those eigenmodes with high eigenvalues, since in the outer region, only numerical noise can couple to these modes in any significant way. We chose to keep only quadrupole terms. Also, implement the outgoing-wave boundary conditions using only monopole and quadrupole terms and conditions derived from the continuum mathematics. * method BF. This "back-n-forth" approach uses FD in inner and outer radial zones of the computational domain, and the eigenvector projection in an intermediate radial zone. The FD method is subject to "resonant cavity" behavior. We think numerical noise originates from a mismatch between the grid geometry and the outgoing wave symmetry can get amplified, sometimes producing a ratty solution, depending on the grid details. The poor solutions typically have spurious high-frequency components. The EV method seems to work well in general. No "resonant cavity" behavior was found. In fact, the high-frequency noise can not exist at large radii (r>0.5*r_max, with r_max typically 30a-60a), since only monopole and quadrupole terms are allowed there. Any numerical noise in the EV method appears as a small amount of ingoing radiation that gets amplified a little in the transition zone. The BF method behaves somewhere between the FD and EV approaches. This might be a sign that the noise originates mostly, but not entirely from the outer boundary. There is some anecdotal support for this idea: We generally see improvement in the computed solution compared to the exact when the EV approach is used only to specify the outer boundary, while FD is used elsewhere in the computational domain. This method could be interesting to PSW for a couple of reasons. One is might be related to the spectral approach. I need to think this through, but the eigenmethod might really speed up the ugly integrals we were doing to get projection coefficients needed for the nonlinear case (don't remember how the weight function plays into this?). Another is...well, there must be another reason. Anyway, there is something fun about play with Y_lm's that have non-integral l's and m's. Any thought's on this? we might improve our outer BC's if the "right" m's were used. I can guess an m by comparing with a projection onto the exact solution, maybe.