Package com.nr.pde

Source Code of com.nr.pde.Mgfas

package com.nr.pde;

import static com.nr.NRUtil.*;
import static java.lang.Math.*;

import org.netlib.util.doubleW;

/*"
Full multigrid algorithm for FAS solution of nonlinear elliptic equation, here equation (20.6.44)
on a square domain of side 1, so that h = 1/(n-1)
*/
public class Mgfas {
  int n,ng;
  double[][] uj,uj1;
  //NRvector<NRmatrix<double> *> rho;
  double[][][] rho;   // Vector of pointers to ρ on each level
 
  public double[][] u;

  /**
   *
   * @param u
   * @param maxcyc
   */
  /*
   On input u[0..n-1][0..n-1] contains the right-hand side ρ, while on output it returns the
   solution. The dimension n must be of the form 2^j+1 for some integer j. (j is actually the number
   of grid levels used in the solution, called ng below). maxcyc is the maximum number of V-cycles to be used at each level
   */
  public Mgfas(final double[][] u, final int maxcyc) {
    this.u = u;
    n = u.length;
    ng = 0;
    int nn=n;
    while ((nn >>>= 1)!=0) ng++;
    if ((n-1) != (1 << ng))
      throw new IllegalArgumentException("n-1 must be a power of 2 in mgfas.");
    nn=n;
    int ngrid=ng-1;
    rho = new double[ng][][];   // Allocate storage for r, h, s on grid ng-1
    //rho.resize(ng);                  // and fill it with the input r, h, s
    rho[ngrid]=new double[nn][nn];   // Similarly allocate storage and fill r, h, s by
    copyAssign(rho[ngrid],u);           // restriction on all coarse grids
    while (nn > 3) {
      nn=nn/2+1;
      rho[--ngrid]=new double[nn][nn];
      rstrct(rho[ngrid],rho[ngrid+1]);
    }
    nn=3;
    uj=new double[nn][nn];
    slvsm2(uj,rho[0]);    // Initial solution on coarsest grid
    for(int j=1;j<ng;j++) {  // nested iteration loop
      nn=2*nn-1;
      uj1=uj;
      uj=new double[nn][nn];
      double[][] temp = new double[nn][nn];
      interp(uj,uj1)// interpolate from grid j-1 to next finer grid j
      for (int jcycle=0;jcycle<maxcyc;jcycle++) {        // V-cycle loop
        doubleW trerr= new doubleW(1.0);   // R.h.s is dummy
        mg(j,uj,temp,rho,trerr);
        lop(temp,uj);   // Form residual ||d_h||
        matsub(temp,rho[j],temp);
        double res=anorm2(temp);
        if (res < trerr.val) break;    // No more V-cycles needed if residual small enough
      }
    }
    this.u = uj;   // return solution in u
  }
 
 
  /*
   Matrix addition: Adds a[0..n-1][0..n-1] to b[0..n-1][0..n-1] and returns reslt in c[0..n-1][0..n-1]
   */
  public static void matadd(final double[][] a, final double[][] b, final double[][] c) {
    int n=a.length;
    for (int j=0;j<n;j++)
      for (int i=0;i<n;i++)
        c[i][j]=a[i][j]+b[i][j];
  }
 
  /*
   Matrix subtraction: Subtracts b[0..n-1][0..n-1] from a[0..n-1][0..n-1] and returns result in c[0..n-1][0..n-1]
   */
  public static void matsub(final double[][] a, final double[][] b, final double[][] c) {
    int n=a.length;
    for (int j=0;j<n;j++)
      for (int i=0;i<n;i++)
        c[i][j]=a[i][j]-b[i][j];
  }
 
/*
   Solution of equation (20.6.44) on the coarsest grid, where h = 1/2. The right-hand side is
   input in rhs[0..2][0..2] and the solution is returned in u[0..2][0..2]
   */
  public static void slvsm2(final double[][] u, final double[][] rhs){
    double h=0.5;
    for (int i=0;i<3;i++)
      for (int j=0;j<3;j++)
        u[i][j]=0.0;
    double fact=2.0/(h*h);
    double disc=sqrt(fact*fact+rhs[1][1]);
    u[1][1]= -rhs[1][1]/(fact+disc);
  }
 
  /*
   Red-black Gauss-Seidel relaxation for equation (20.6.44). The current value of the solution
   u[0..n-1][0..n-1] is updated, using the right-hand side function rhs[0..n-1][0..n-1]
   */
  public static void relax2(final double[][] u, final double[][] rhs) {
    int n=u.length;
    int jsw=1;
    double h=1.0/(n-1);
    double h2i=1.0/(h*h);
    double foh2 = -4.0*h2i;
    for (int ipass=0;ipass<2;ipass++,jsw=3-jsw) {
      int isw=jsw;
      for (int j=1;j<n-1;j++,isw=3-isw) {
        for (int i=isw;i<n-1;i+=2) {
          double res=h2i*(u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]-
            4.0*u[i][j])+u[i][j]*u[i][j]-rhs[i][j];
          u[i][j] -= res/(foh2+2.0*u[i][j]);
        }
      }
    }
  }
 
  /*
   Half-weighting restriction. If nc is the coarse-grid dimension, the fine-grid solution is input
   in uf[0..2*nc-2][0..2*nc-2]. The coarse-grid solution obtained by restriction is returned in uc[0..nc-1][0..nc-1]
   */
  public static void rstrct(final double[][] uc, final double[][] uf) {
    int nc=uc.length;
    int ncc=2*nc-2;
    for (int jf=2,jc=1;jc<nc-1;jc++,jf+=2) {   // Interior points
      for (int iif=2,ic=1;ic<nc-1;ic++,iif+=2) {
        uc[ic][jc]=0.5*uf[iif][jf]+0.125*(uf[iif+1][jf]+uf[iif-1][jf]
          +uf[iif][jf+1]+uf[iif][jf-1]);
      }
    }
    for (int jc=0,ic=0;ic<nc;ic++,jc+=2) {   // Boundary points
      uc[ic][0]=uf[jc][0];
      uc[ic][nc-1]=uf[jc][ncc];
    }
    for (int jc=0,ic=0;ic<nc;ic++,jc+=2) {
      uc[0][ic]=uf[0][jc];
      uc[nc-1][ic]=uf[ncc][jc];
    }
  }
 
  /*
   Given u[0..n-1][0..n-1], returns \phi_h(u_h) for eqn. (20.6.44) in ot[0..n-1][0..n-1]
   */
  public static void lop(final double[][] out, final double[][] u) {
    int n=u.length;
    double h=1.0/(n-1);
    double h2i=1.0/(h*h);
    for (int j=1;j<n-1;j++)    // Interior points
      for (int i=1;i<n-1;i++)
        out[i][j]=h2i*(u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]-
          4.0*u[i][j])+u[i][j]*u[i][j];
    for (int i=0;i<n;i++)    // Boundary points
      out[i][0]=out[i][n-1]=out[0][i]=out[n-1][i]=0.0;
  }
 
  /*
   Coarse-to-fine prolongation by bilinear interpolation. If nf is the fine-grid dimension, the coarse-grid solution is
   input as uc[0..nc-1][0..nc-1], where nc = nf/2+1. The fine-grid solution is returned in f[0..nf-1][0..nf-1]
   */
  public static void interp(final double[][] uf, final double[][] uc){
    int nf=uf.length;
    int nc=nf/2+1;
    for (int jc=0;jc<nc;jc++)   // Do elements that are copies
      for (int ic=0;ic<nc;ic++) uf[2*ic][2*jc]=uc[ic][jc];
    for (int jf=0;jf<nf;jf+=2)   // Do even-numbered columns, interpolating verically
      for (int iif=1;iif<nf-1;iif+=2)
        uf[iif][jf]=0.5*(uf[iif+1][jf]+uf[iif-1][jf]);
    for (int jf=1;jf<nf-1;jf+=2)   // Do odd-numbered columns, interpolating horizontally
      for (int iif=0;iif<nf;iif++)
        uf[iif][jf]=0.5*(uf[iif][jf+1]+uf[iif][jf-1]);
  }
 
  // Returns the Euclidean norm of the matrix a[0..n-1][0..n-1]
  public static double anorm2(final double[][] a) {
    double sum=0.0;
    int n=a.length;
    for (int j=0;j<n;j++)
      for (int i=0;i<n;i++)
        sum += a[i][j]*a[i][j];
    return sqrt(sum)/n;
  }
 
  /*
   Recursive multigrid iteration. On input, j is the current level and u is the current value of the solution.
   For the first call on a given level, the right-hand side is zero, and the argument rhs is dummy.
   This is signaled by inputting trerr positive. Subsequent recursive calls supply
   a nonzero rhs as in equation (20.6.33). Tis is signaled by inputting trerr negative. rho is the vector
   of pointers to ρ on each level. On output u contains the improved solution at the current level. When the
   first call on a given level is made, the relative truncation error is returned in trerr
   */
  public static void mg(final int j, final double[][] u, final double[][] rhs,
    final double[][][]rho, final doubleW trerr) {
    final int NPRE=1,NPOST=1;
    // Number of relaxation sweeps before and after the coarse-grid correction is computed
    final double ALPHA=0.33;   // Relates the estimated truncation error to the norm of the residual
    doubleW dum= new doubleW(-1.0);
    int nf=u.length;
    int nc=(nf+1)/2;
    double[][] temp = new double[nf][nf];
    if (j == 0) {   // Bottom of V:  Solve on coarsest grid.
      matadd(rhs,rho[j],temp);
      slvsm2(u,temp);
    } else {   // On downward stoke of the V
      double[][] v= new double[nc][nc],ut=new double[nc][nc],tau=new double[nc][nc],tempc=new double[nc][nc];
      for (int jpre=0;jpre<NPRE;jpre++) {   // Pre-smoothing
        if (trerr.val < 0.0) { 
          matadd(rhs,rho[j],temp);
          relax2(u,temp);
        }
        else
          relax2(u,rho[j]);
      }
      rstrct(ut,u);
      copyAssign(v, ut);
      lop(tau,ut);
      lop(temp,u);
      if (trerr.val < 0.0)
        matsub(temp,rhs,temp);
      rstrct(tempc,temp);
      matsub(tau,tempc,tau);
      if (trerr.val > 0.0)
        trerr.val=ALPHA*anorm2(tau);
      mg(j-1,v,tau,rho,dum);
      matsub(v,ut,tempc);
      interp(temp,tempc);
      matadd(u,temp,u);
      for (int jpost=0;jpost<NPOST;jpost++) {
        if (trerr.val < 0.0) {
          matadd(rhs,rho[j],temp);
          relax2(u,temp);
        }
        else
          relax2(u,rho[j]);
      }
    }
  }
}
TOP

Related Classes of com.nr.pde.Mgfas

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.