require 'narray' require 'numru/gphys' require 'numru/derivative' ############################################################ =begin =module NumRu::GPhys::Derivative in derivative.rb ==Index * (()) * (()) * First derivative (2nd Order difference use three point.) * (()) * First derivative (using center difference method) =module NumRu::GPhys::Derivative Module functions of Derivative Operater for GPhys. ---threepoint_O2nd_deriv(gp, dim_or_dimname, bc=NumRu::Derivative::LINEAR_EXT)) Derivate (({z})) respect to (({dim})) th dimension with 2nd Order difference. return an NArray which result of the difference (()) divided difference (({x})) (in other wards, (s**2*z_{i+1} + (t**2 - s**2)*f_{i} - t**2*f_{i-1}) / (s*t*(s + t)): now s represents (x_{i} - x_{i-1}) ,t represents (x_{i+1} - x_{i}) and _{i} represents the suffix of {i} th element in the (()) th dimension of array. ). ARGUMENTS * z (NArray): a NArray which you want to derivative. * x (NArray): a NArray represents the dimension which derivative respect to. z.rank must be 1. * dim (Numeric): a Numeric represents the dimention which derivative respect to. you can give number count backward ((())<0), but (()) must be > 0. * bc (Numeric) : a Numeric which represent boundary condition. now only LINEAR_EXT(=1) supported. LINEAR_EXT load (()) which extend boundary with lenear value. RETURN VALUE * O2nd_deriv_data (NArray): (s**2*z_{i+1} + (t**2 - s**2)*f_{i} - t**2*f_{i-1}) / (s*t*(s + t)) ---cderiv(gp, dim_or_dimname, bc=NumRu::Derivative::LINEAR_EXT) Derivate (({gp})) respect to (({dim})) th or (({dimname})) dimension with center difference. return a GPhys which result of the difference (()) divided difference axis. ( in other wards, (z_{i+1} - z_{i-1}) / (x_{i+1} - x_{i-1}): now x is axis array which you wants to respects to, _{i} represents the suffix of {i} th element in the (()) th dimension of array. ). ARGUMENTS * gp (GPhys): a GPhys which you want to derivative. * dim_or_dimname (Numeric or String): a Numeric or String represents the dimension which derivate respect to. you can give number count backward ((())<0), but (()) must be > 0. * bc (Numeric) : a Numeric which represent boundary condition. now only NumRu::Derivative::LINEAR_EXT(=1) supported. LINEAR_EXT load (()) which extend boundary with linear value. RETURN VALUE * a GPhys =end ############################################################ module NumRu class GPhys module Derivative module_function def threepoint_O2nd_deriv(gp,dim,bc=NumRu::Derivative::LINEAR_EXT) # <> if (dim.is_a?(Numeric) && dim < 0) dim += gp.rank elsif dim.is_a?(String) dim = gp.axnames.index(dim) end # <> v_data = gp.data; v_x = gp.coord(dim) # get varray n_data = v_data.val; n_x = v_x.val # get narray n_dgpdx = NumRu::Derivative::threepoint_O2nd_deriv(n_data,n_x,dim,bc) # derivate n_data name = "d#{gp.name}_d#{v_x.name}" # ex. "dT_dx" v_dgpdx = VArray.new(n_dgpdx, gp.data, name) # make varray g_dgpdx = GPhys.new( gp.grid_copy, v_dgpdx ) # make gphys # <> u_data = v_data.units; u_x = v_x.units # get units g_dgpdx.units = u_data/u_x # set units if v_data.get_att('long_name') && v_x.get_att('long_name') long_name = "d_#{v_data.get_att('long_name')}_d_#{v_x.get_att('long_name')}" else long_name = name end g_dgpdx.data.set_att("long_name",long_name) # set long_name return g_dgpdx end def cderiv(gp,dim,bc=NumRu::Derivative::LINEAR_EXT) # <> if (dim.is_a?(Numeric) && dim < 0) dim += gp.rank elsif dim.is_a?(String) dim = gp.axnames.index(dim) end # <> v_data = gp.data; v_x = gp.coord(dim) # get varray n_data = v_data.val; n_x = v_x.val # get narray n_dgpdx = NumRu::Derivative::cderiv(n_data,n_x,dim,bc) # derivate n_data name = "d#{gp.name}_d#{v_x.name}" # ex. "dT_dx" v_dgpdx = VArray.new(n_dgpdx, gp.data, name) # make varray g_dgpdx = GPhys.new( gp.grid_copy, v_dgpdx ) # make gphys # <> u_data = v_data.units; u_x = v_x.units # get units g_dgpdx.units = u_data/u_x # set units if v_data.get_att('long_name') && v_x.get_att('long_name') long_name = "d_#{v_data.get_att('long_name')}_d_#{v_x.get_att('long_name')}" else long_name = name end g_dgpdx.data.set_att("long_name",long_name) # set long_name return g_dgpdx end end end end ###################################################### ## < test > if $0 == __FILE__ include NumRu include NMath #### define method for tests. ## # genarate narray def gen_x(n) 2*PI*NArray.sfloat(n).indgen!/(n-1) end def gen_y(n) PI*(NArray.sfloat(n).indgen!/(n-1)-0.5) end def gen_z1(n) NArray.sfloat(n).indgen!/(n-1) end def gen_z2(n) exp(-NArray.sfloat(n).indgen!/(n-1)) end def make_gp1D(f,x) ax = Axis.new.set_pos( VArray.new( x , {"long_name"=>"longitude", "units"=>"rad"}, "lon" )) data = VArray.new( f, {"long_name"=>"temperature", "units"=>"K"}, "t" ) return GPhys.new( Grid.new(ax), data) end def make_gp3D(f,x,y,z) vax = VArray.new( x, {"long_name"=>"longitude", "units"=>"rad"}, "lon" ) vay = VArray.new( y, {"long_name"=>"latitude", "units"=>"rad"}, "lat" ) vaz = VArray.new( z, {"long_name"=>"altitude", "units"=>"m"}, "z" ) axx = Axis.new.set_pos(vax) axy = Axis.new.set_pos(vay) axz = Axis.new.set_pos(vaz) data = VArray.new( f, {"long_name"=>"temperature", "units"=>"K"}, "t" ) return GPhys.new( Grid.new(axx, axy, axz), data) end def show_attr(bef_deriv, aft_deriv) fm = "%-15s%-15s%-10s%s" printf(fm, "", "", "", "\n") printf(fm, "name", bef_deriv.data.name, aft_deriv.data.name, "\n") aft_deriv.data.att_names.each{|nm| printf(fm, nm, bef_deriv.data.get_att(nm).to_s, aft_deriv.data.get_att(nm).to_s, "\n") } end ## # test def test1(x1, dim_or_dimname) f1 = sin(x1) gp1 = make_gp1D(f1, x1) deriv = GPhys::Derivative::cderiv(gp1, dim_or_dimname) dfdx1 = deriv.data.val dfdx2 = cos(x1) p(dfdx1) if $VERBOSE diff = (dfdx1 - dfdx2)[1..-2].abs err = diff.mean print "dfdx - kaiseki_kai (except boundary): " print err, "\t", diff.max,"\n" err end def test2(x2, dim_or_dimname) f2 = sin(x2) gp2 = make_gp1D(f2, x2) deriv = GPhys::Derivative::threepoint_O2nd_deriv(gp2, dim_or_dimname) dfdx1 = deriv.data.val dfdx2 = cos(x2) p(dfdx1) if $VERBOSE diff = (dfdx1 - dfdx2)[1..-2].abs err = diff.mean print "dfdx - kaiseki_kai (except boundary): " print err, "\t", diff.max,"\n" err end def test3 x3 = gen_x(11) f3 = sin(x3) gp3 = make_gp1D(f3, x3) deriv = GPhys::Derivative::cderiv(gp3, "lon") dfdx1 = deriv.data.val dfdx2 = cos(x3) p(dfdx1) if $VERBOSE show_attr(gp3, deriv) end def test4(dim) nx = 10; ny = 5; nz = 5 x = gen_x(nx) y = gen_y(ny) z1 = gen_z1(nz) f = sin(x).reshape(nx,1,1) * sin(y).reshape(1,ny,1) * sin(z1).reshape(1,1,nz) gp = make_gp3D(f, x, y, z1) print "**** equally spaced grid (2nd order) ****\n" deriv = GPhys::Derivative::threepoint_O2nd_deriv(gp, dim) dfdx1 = deriv.data.val dfdx2 = sin(x).reshape(nx,1,1) * cos(y).reshape(1,ny,1) * sin(z1).reshape(1,1,nz) p(dfdx1) if $VERBOSE diff = (dfdx1 - dfdx2)[1..-2].abs err = diff.mean print "dfdx - kaiseki_kai (except boundary): " print err, "\t", diff.max,"\n" print "**** check attribute ****\n" show_attr(gp, deriv) end ## main routine of test --------------------------------------------- [0, "lon"].each do |dim| print "************** dim.type === #{dim.type} **************\n" print "**** equally spaced grid ****\n" er1 = test1( gen_x(11), dim ) er2 = test1( gen_x(21), dim ) print "error change from nx=11->21: ", er2/er1,"\n" print "**** equally spaced grid (2nd order) ****\n" er1 = test2( gen_x(11), dim ) er2 = test2( gen_x(21), dim ) print "error change from nx=11->21: ", er2/er1,"\n" print "**** non-uniform grid ****\n" p 'x(11):',gen_z2(11),'x(21):',gen_z2(21) er1 = test1( gen_z2(11), dim ) er2 = test1( gen_z2(21), dim ) print "error change from nx=11->21: ", er2/er1,"\n" print "**** non-uniform grid (2nd order) ****\n" er1 = test2( gen_z2(11), dim ) er2 = test2( gen_z2(21), dim ) print "error change from nx=11->21: ", er2/er1,"\n" end print "************** multi-D *************\n" print "******** dimname == 'lat' ********\n" test4("lat") print "******** dim == -1(z) ********\n" test4(-1) end