Module: NumRu::GAnalysis::PDE

Defined in:
lib/numru/ganalysis/pde.rb,
ext/numru/ganalysis/pde_ext.c

Class Method Summary collapse

Class Method Details

.SOR_2D_2ndorder(z, a, b, c, d, e, f, dx, dy, ome, eps: nil, maxi: nil, ignore_non_eliptic: false) ⇒ Object

Solve 2-dim 2nd order PDE by the succesive over-relaxation method

ARGUMENTS.

  • z [in-out, 2D dfloat NArray] : When in, initial values. When out, the result. Boundary values are not changed (Dirichlet problem) (developper’s memo: Neuman cond can be supported as an option)

  • a, b, c, d, e, f [in, 2D dfloat NArray with the same shape as z] : Specifies the PDF as

    a z_xx + 2*b z_xy + c z_yy + d z_x + e z_y = f
         ^^^

    Please note that 2*b is the coefficient of z_xy. x is the first (0th) NArray dimension.

  • dx, dy [Float] : grid spacings

  • ome [Float] : the over-relaxation parameter 1<=ome<2. For a large-scale poisson equation, should be ~ 2 (e.g. 1.9).

  • eps [optional, Float] : conversion threshold (e.g. 1e-5). Compared with |increment|/|F| (| | is the L2 norm). When |f| == 0 (homogenous), compared simply with |increment|.

  • maxi [optional, Integer]: max iteration times(default, 5*.max). Since SOR (or the Gauss-Seidel method) acts like a stepwise diffusion, it must be greater than the along-x and along-y grid sizes.

  • ignore_non_eliptic : if true (non-nil,non-false), do not raise error when the coeeficients are not entirely eliptic.

RETURN VALUE

  • |last increment|/|F| of |last increment| (to be compared with eps)



36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
 
# File 'lib/numru/ganalysis/pde.rb', line 36

def SOR_2D_2ndorder(z, a, b, c, d, e, f, dx, dy, ome, eps: nil, maxi: nil,
                    ignore_non_eliptic: false)
  lf = Math.sqrt( (f**2).sum )  # L2 norm of f
  lf = 1.0 if lf==0   # when |f|==0 (homogeneous), no-division by |f|.
  res = nil
  if !maxi
    maxi = [ z.shape[0], z.shape[1] ].max * 5  # default max itr times
  end
  convd = false
  for i in 0...maxi
    res = SOR_2D_2ndorder_1step(z, a, b, c, d, e, f, dx, dy, ome,
                                ignore_non_eliptic)
    #p res/lf, z
    if (eps && res/lf < eps)
      convd = true
      break
    end
  end
  if (eps && !convd)
    raise("Max iteration times (#{maxi}) reached before convergence" +
          " (#{res/lf} > #{eps}). Change ome (#{ome}) or specify maxi.")
  end
  #p i,res/lf #, z
  res/lf
end

.SOR_2D_2ndorder_1step ⇒ Object

Eq: A z_xx + 2 B z_xy + C z_yy + D z_x + E z_y = F

^^^

Assumption: dx and dy are constants (uniform). Boundary values are not changed. –> change afterwords if not Dirichlet.

Basic discretization: by central differences:

A/dx^2 (z[i-1,j]-2z[i,j]+z[i+1,j])
+ B/(2dxdy) ((z[i+1,j+1]-z[i-1,j+1]-z[i+1,j-1]+z[i-1,j-1])
+ C/dy^2 (z[i,j-1]-2z[i,j]+z[i,j+1])
+ D/(2dx) (z[i+1,j]-z[i-1,j]) + E/(2dy) (z[i,j+1]-z[i,j-1]) - F = 0

  • Z [2D Float(=DFLOAT) NArray] (inout) – to be overwritten

  • A, B, … [2D Float(=DFLOAT) NArray] coeffitients

  • OME [Float] constant, which should 1<= OME < 2 (Gauss-Seidel if OME=1)

  • Return value: res = |H dz| (L2 norm), where H = 2(diag(A)/dx^2+diag©/dy^2). Convergence can be judged by res/|F| < eps



37
38
39
 
# File 'ext/numru/ganalysis/pde_ext.c', line 37

static VALUE
SOR_2D_2ndorder_1step(obj, Z, A, B, C, D, E, F, DX, DY, OME, Ignore_non_eliptic)
VALUE obj;