Changeset 747


Ignore:
Timestamp:
10/05/18 13:42:49 (6 years ago)
Author:
dubos
Message:

devel/unstructured : more fixes to mixed precision

Location:
codes/icosagcm/devel
Files:
1 added
15 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/Python/dynamico/meshes.py

    r689 r747  
    1111import parallel 
    1212from util import list_stencil, Base_class, inverse_list 
     13from precision import zeros 
    1314 
    1415radian=180/math.pi 
    15  
    16 def zeros(dims): return np.zeros([n for n in dims if n>1], dtype=unst.np_num) 
    17  
     16   
    1817#------------------- Hybrid mass-based coordinate ------------- 
    1918 
    … …  
    170169                    (indexu(x,y+1),.25), 
    171170                    (indexu(x-1,y+1),.25))) 
    172         Aiv=np.zeros(nx*ny)+dx*dy # Ai=Av=dx*dy 
     171        Aiv=np.zeros(nx*ny) 
     172        Aiv[:]=dx*dy # Ai=Av=dx*dy 
    173173        unst.init_mesh(llm,nqdyn,2*nx*ny,nx*ny,nx*ny,4,4,4, 
    174174                  primal_nb,primal_edge,primal_ne, 
    … …  
    181181    def field_z(self,n=1): return zeros((n,self.ny,self.nx,self.llm)) 
    182182    def field_w(self,n=1): return zeros((n,self.ny,self.nx,self.llm+1)) 
    183     def field_u(self,n=1):  
    184         if n==1: 
    185             return np.zeros((self.ny,2*self.nx,self.llm)) 
    186         else: 
    187             return np.zeros((n,self.ny,2*self.nx,self.llm)) 
     183    def field_u(self,n=1): return zeros((self.ny,2*self.nx,self.llm)) 
    188184    def field_ps(self,n=1): return zeros((n,self.ny,self.nx)) 
    189     def ucomp(self,u): return u[:,range(0,2*self.nx,2),:] 
    190     def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u 
     185    def ucomp(self,u):  
     186        return u[range(0,2*self.nx,2),:] if self.ny==1 else u[:,range(0,2*self.nx,2),:] 
     187    def set_ucomp(self,uv,u): 
     188        if self.ny==1 : uv[range(0,2*self.nx,2),:]=u 
     189        else : uv[:,range(0,2*self.nx,2),:]=u 
    191190    def vcomp(self,u): return u[:,range(1,2*self.nx,2),:] 
    192191 
    … …  
    288287        return self.de*(ulon*np.cos(self.angle_e)+ulat*np.sin(self.angle_e)) 
    289288    def ucov3D(self, ulon, ulat): 
    290         ucov = np.zeros((self.edge_num,ulon.shape[1])) 
     289        ucov = zeros((self.edge_num,ulon.shape[1])) 
    291290        for edge in range(self.edge_num): 
    292291            angle=self.angle_e[edge] 

  • codes/icosagcm/devel/Python/dynamico/time_step.py

    r615 r747  
    8585    def cjl(self,b): # converts Butcher tableau b(j,l) to update coefficients c(j,l) 
    8686        return np.array([ b[j+1,:]-b[j,:] for j in range(self.nstage)]) 
    87     def __init__(self,bwd_fast_slow, dt,nstage, bsjl,bfjl): 
     87    def __init__(self,precision,bwd_fast_slow, dt,nstage, bsjl,bfjl): 
     88        if bwd_fast_slow is None: 
     89            print 'SIRK with Fortran-side time stepping' 
     90            precision=np.float64 
     91        else: 
     92            print 'SIRK with Python-side time stepping at precision ', precision 
     93 
    8894        self.bwd_fast_slow = bwd_fast_slow 
    8995        self.dt, self.nstage = dt, nstage 
    90         self.csjl = dt*self.cjl(bsjl) 
    91         self.cfjl = dt*self.cjl(bfjl) 
    92         self.tauj = np.array([dt*bfjl[j,j] for j in range(self.nstage)]) 
     96        csjl = dt*self.cjl(bsjl) 
     97        cfjl = dt*self.cjl(bfjl) 
     98        self.csjl, self.cfjl = [x.astype(precision) for x in csjl, cfjl] 
     99        self.tauj = np.array([dt*bfjl[j,j] for j in range(self.nstage)], dtype=precision) 
     100        print 'Types of csjl, cfjl, tauj :', self.csjl.dtype, self.cfjl.dtype, self.tauj.dtype 
    93101    def next(self,zj): 
    94102        csjl, cfjl, tauj = self.csjl, self.cfjl, self.tauj 
    … …  
    100108            for l in range(0,j+1): # here zj is initially y(j) and finally z(j+1) 
    101109                zj = [z + csjl[j,l]*s + cfjl[j,l]*f for z,s,f in zip(zj, slowj[l], fastj[l]) ] 
     110#            print 'SIRK : zj has type ', [ z.dtype for z in zj] 
    102111        return zj         
    103112    def advance(self, flow, N): 
    … …  
    113122    return SIRK(bwd_fast_slow,dt,1,bsjl,bfjl) 
    114123         
    115 def ARK2(bwd_fast_slow, dt, a32 = (3.+2.*np.sqrt(2.))/6.): 
     124def ARK2(bwd_fast_slow, dt, precision=None, a32 = (3.+2.*np.sqrt(2.))/6.): 
    116125    """ ARK2 scheme by Giraldo et al. (2013), noted ARK2(2,3,2) in Weller et al. (2013) """  
    117126    gamma = 1.-1./np.sqrt(2.) 
    … …  
    120129    bsjl = np.array([[0,0,0],[2.*gamma,0,0],[1.-a32,a32,0],wj]) 
    121130    bfjl = np.array([[0,0,0],[gamma,gamma,0],[delta,delta,gamma],wj]) 
    122     return SIRK(bwd_fast_slow,dt,3,bsjl,bfjl) 
     131    return SIRK(precision,bwd_fast_slow,dt,3,bsjl,bfjl) 
    123132         
    124133def UJ3(bwd_fast_slow, dt): 
    … …  
    135144######### Explicit RK schemes written as ARK with two identical Butcher tableaus ######## 
    136145 
    137 def RK4(bwd_fast_slow, dt): 
     146def RK4(bwd_fast_slow, dt, precision=None): 
    138147    """ 4-stage 4th-order RK """  
    139148    bjl = np.array([[0,0,0,0], 
    … …  
    142151                    [0.,0.,1.,0], 
    143152                    [1/6.,1/3.,1/3.,1/6.]]) 
    144     return SIRK(bwd_fast_slow,dt,4,bjl,bjl) 
     153    return SIRK(precision, bwd_fast_slow,dt,4,bjl,bjl) 
    145154 
    146155######### Simplified explicit RK schemes (2nd-order) ######## 

  • codes/icosagcm/devel/Python/dynamico/wrap.py

    r663 r747  
    1111# defines converter functions depending on type of Python argument 
    1212def noop(x): return x 
    13 py2c = { np.ndarray:loc, np.float64:c_double, float:c_double, int:c_int,  
     13py2c = { np.ndarray:loc, np.float64:c_double, np.float32:c_float, int:c_int,  
    1414         c_void_p:noop, str:noop, c_int:noop } 
    1515 

  • codes/icosagcm/devel/Python/src/unstructured.pyx

    r746 r747  
    239239        # h*s = h => uniform buoyancy s=1 => shallow-water 
    240240        dh, du_slow, du_fast, hs, buf = self.fmass(), self.fu(), self.fu(), h.copy(), self.geopot 
     241        assert type(tau) is np_num, 'tau must be of type unstructured.np_num.' 
    241242        ker.dynamico_caldyn_unstructured(tau, self.ms, h, hs, u, self.geopot, buf, 
    242243                  self.s, self.ps, self.pk, self.hflux, self.qv, 
    … …  
    296297              left,right,down,up,trisk_deg,trisk, 
    297298              Ai, Av, fv, le_de, Riv2, wee): 
     299 
     300    print 'Types of Ai, Av, fv, le_de, Riv2, wee : ', Ai.dtype,Av.dtype,fv.dtype,le_de.dtype,Riv2.dtype,wee.dtype 
     301    for var,varname in zip((Ai,Av,fv,le_de,Riv2,wee), ('Ai','Av','fv','le_de','Riv2','wee')): 
     302        assert var.dtype == np.float64, '%s must be double precision'%varname 
     303 
    298304    setvars( ('llm','nqdyn','edge_num','primal_num','dual_num', 
    299305              'max_trisk_deg','max_primal_deg','max_dual_deg'), 

  • codes/icosagcm/devel/Python/test/py/DCMIP2008c5.py

    r687 r747  
    136136plt.savefig('fig_DCMIP2008c5/Ai.png'); plt.close() 
    137137 
    138  
    139 if False: # time stepping in Python 
     138if False: # time stepping in Python => set precision to np_num 
    140139    caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,thermo,params,params.g) 
    141     scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7) 
     140    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, precision=unst.np_num, a32=0.7) 
    142141    def next_flow(m,S,u):  
    143142        m,S,u = scheme.advance((m,S,u),nt) 
    … …  
    179178z, vol = mesh.field_mass(), mesh.field_mass() 
    180179m,S,u,Phi,W = flow0  
     180 
     181print 'types of m,S,u :', m.dtype, S.dtype, u.dtype 
     182 
    181183caldyn_step.geopot[:,0]=Phi[:,0] 
    182184plots(0) 

  • codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py

    r656 r747  
    44from dynamico import DCMIP 
    55from dynamico import meshes 
     6from dynamico import precision as prec 
    67from dynamico.meshes import Cartesian_mesh as Mesh 
    78import math as math 
    … …  
    5152    params.rho_bot = 1e6/gas_bot.v 
    5253     
    53     return thermo, mesh, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas 
     54    return thermo, mesh, params, prec.asnum([mass_ik,Sik,ujk,Phi_il,Wil]), gas 
    5455 
    5556#Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.8 
    … …  
    6364# compute hybrid coefs from initial distribution of mass 
    6465mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0]) 
     66print 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk] 
    6567unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk) 
    6668 

  • codes/icosagcm/devel/Python/test/py/RSW_2D.py

    r632 r747  
    1818 
    1919x1,x2,yy = mesh.xx[:,:,0]-1., mesh.xx[:,:,0]+1., mesh.yy[:,:,0] 
    20 h0 = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+ 
     20h0, u0 = mesh.field_mass(), mesh.field_u()  # flow initally at rest 
     21h0[:] = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+ 
    2122            np.exp(-2.*(x2*x2+yy*yy))) 
    22 u0 = mesh.field_u()  # flow initally at rest 
    2323flow0=(h0,u0) 
    2424 
    … …  
    3232scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt) 
    3333 
     34print 'types : yy, u0, h0, u0', yy.dtype, u0.dtype, h0.dtype, u0.dtype 
     35 
    3436flow=flow0 
    3537for i in range(10): 
    3638    h,u=flow 
     39    print 'types : h,u', h.dtype, u.dtype 
    3740    caldyn.bwd_fast_slow(flow, 0.) 
    3841    plt.figure(); plt.pcolor(mesh.x,mesh.y,caldyn.qv) 
    … …  
    4851        unst.elapsed=0. 
    4952        flow = scheme.advance(flow,N) 
     53        h,u=flow 
     54        print 'types : h,u', h.dtype, u.dtype 
    5055        # flops 
    5156        # mass flux : 1add+1mul per edge          =>  4 

  • codes/icosagcm/devel/Python/test/py/RSW_2D_mesh.py

    r719 r747  
    4545dx,dy=Lx/nx,Ly/ny 
    4646 
    47 filename, llm, nqdyn, g, f, radius = 'cart128.nc', 1, 1, 1., 1., 1. 
     47filename, llm, nqdyn, g, f, radius = 'cart128.nc', 1, 1, 1., 1., None 
    4848unst.setvar('g',g) 
    4949 
    … …  
    6464 
    6565x1,x2,yy = xx-1., xx+1., yy 
     66u0 = mesh.field_u()  # flow initally at rest 
    6667h0 = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+ 
    6768            np.exp(-2.*(x2*x2+yy*yy))) 
    68 u0 = mesh.field_u()  # flow initally at rest 
    69 flow0=(h0,u0) 
     69flow0=meshes.asnum([h0,u0]) 
    7070 
    7171#print('h0 : ', h0[1234]-1.) 
    … …  
    7979dt=T/N 
    8080print N,dt,Lx/nx 
    81 scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt) 
     81#scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, precision=unst.np_num) 
     82scheme = time_step.RK4(caldyn.bwd_fast_slow, dt, precision=unst.np_num) 
    8283 
    8384flow=flow0 
    … …  
    9091for i in range(10): 
    9192    h,u=flow 
    92     flow, fast, slow = caldyn.bwd_fast_slow(flow, 0.) 
     93    flow, fast, slow = caldyn.bwd_fast_slow(flow, unst.zero_num) 
    9394    junk, du_fast = fast 
    9495    dh, du_slow = slow 
    … …  
    102103    plt.figure(); plt.pcolor(x,y,reshape(caldyn.qv)) 
    103104    plt.colorbar(); plt.title('potential vorticity') 
    104     plt.savefig('fig_RSW_2D_mesh/%02d.png'%i) 
     105    plt.savefig('fig_RSW_2D_mesh/%03d.png'%i) 
     106    plt.close() 
    105107#    plt.figure(); plt.pcolor(mesh.x,mesh.y,h) 
    106108#    plt.colorbar(); plt.title('h') 
    … …  
    112114        unst.elapsed=0. 
    113115        flow = scheme.advance(flow,N) 
     116#        flow = scheme.advance(flow,5) 
    114117        # flops 
    115118        # mass flux : 1add+1mul per edge          =>  4 

  • codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py

    r680 r747  
    11print 'Loading DYNAMICO modules ...' 
    22from dynamico import unstructured as unst 
     3from dynamico.precision import asnum 
    34from dynamico.meshes import MPAS_Format, Unstructured_Mesh 
    45from dynamico import time_step 
    … …  
    3435dt = T/nt 
    3536#scheme = time_step.RKn_simple(1,caldyn.bwd_fast_slow, dt) 
    36 scheme = time_step.RK4(caldyn.bwd_fast_slow, dt) 
     37scheme = time_step.RK4(caldyn.bwd_fast_slow, dt, precision=unst.np_num) 
    3738     
    3839print dx, dt, dt*c0/dx 
    … …  
    6162fu_perp = mesh.ucov2D(0.*ulon,2*Omega*np.sin(mesh.lat_e)*ulon) 
    6263 
    63 flow = (Phi0,mesh.ucov2D(ulon,0.*ulon)) 
     64flow = asnum([Phi0, mesh.ucov2D(ulon,0.*ulon)]) 
     65print 'type of Phi0,u0 : ', flow[0].dtype, flow[1].dtype 
     66 
    6467for i in range(N): 
    6568    h,u=flow 

  • codes/icosagcm/devel/Python/test/py/bubble.py

    r631 r747  
    1 import math as math 
    2 import numpy as np 
    3 import matplotlib.pyplot as plt 
    4 import matplotlib.animation as manimation 
    5  
    61from dynamico.meshes import Cartesian_mesh as Mesh 
    72import dynamico.dyn as dyn 
    … …  
    105import dynamico.advect as advect 
    116from dynamico import unstructured as unst 
     7 
     8import math as math 
     9import numpy as np 
     10import matplotlib.pyplot as plt 
     11import matplotlib.animation as manimation 
     12 
    1213 
    1314def reldiff(a,b): 

  • codes/icosagcm/devel/Python/test/py/slice_GW_NH.py

    r666 r747  
    2929 
    3030    u0=mesh.field_u()  
    31     u0[:,range(0,2*nx,2),:] = u*mesh.dx # Doppler shift by u 
     31    u0[range(0,2*nx,2),:] = u*mesh.dx # Doppler shift by u 
    3232    flow0=(m0ik,S0ik,u0,Phi0_il,0*Phi0_il) 
    3333 
    … …  
    8181        plt.figure(figsize=(12,3)) 
    8282        ucomp = mesh.ucomp(u) 
    83         plt.pcolor(xx,zz,ucomp[0,:,:]/dx) 
     83        plt.pcolor(xx,zz,ucomp[:,:]/dx) 
    8484        plt_axes() 
    8585        plt.title("u(t = %f s)" % ((i+1)*dt*nt)) 

  • codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py

    r656 r747  
    2323 
    2424    u0=mesh.field_u() 
    25     u0[:,range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s 
     25    u0[range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s 
    2626    flow0=(m0ik,S0ik,u0) 
    2727 
    … …  
    6060        plt.figure(figsize=(12,3)) 
    6161        ucomp = mesh.ucomp(u) 
    62         plt.pcolor(xx,zz,ucomp[0,:,:]/dx) 
     62        plt.pcolor(xx,zz,ucomp[:,:]/dx) 
    6363        plt.ylim(0,10.) 
    6464        plt.colorbar() 

  • codes/icosagcm/devel/make_icosa

    r744 r747  
    227227  ./build --job $job 
    228228fi 
     229 
     230# force recompilation of cython modules 
     231touch Python/src/*.pyx 

  • codes/icosagcm/devel/make_python

    r721 r747  
    3232        update "$x" "$DYNAMICO_ROOT/src/$2/$x" 
    3333    done 
     34 
    3435} 
    3536 
    3637function cmd_kernels() # this function is invoked by : './make_python kernels' 
    3738{ 
     39    cmd_clean 
    3840    cd $KERNELS 
    3941    ./codegen hexagonal hex_master unstructured 

  • codes/icosagcm/devel/src/unstructured/data_unstructured.F90

    r744 r747  
    205205    Riv2(:,:)=Riv2_(:,:) 
    206206    wee(:,:) = wee_(:,:) 
    207     PRINT *, MAXVAL(ABS(Ai)) 
    208     PRINT *, MAXVAL(ABS(Av)) 
    209     PRINT *, MAXVAL(ABS(fv)) 
    210     PRINT *, MAXVAL(ABS(le_de)) 
    211     PRINT *, MAXVAL(ABS(Riv2)) 
    212     PRINT *, MAXVAL(ABS(wee)) 
    213     PRINT *, MINVAL(right), MAXVAL(right) 
    214     PRINT *, MINVAL(right), MAXVAL(left) 
     207    PRINT *, 'Max Ai : ',    MAXVAL(ABS(Ai)) 
     208    PRINT *, 'Max Av : ',    MAXVAL(ABS(Av)) 
     209    PRINT *, 'Max fv : ',    MAXVAL(ABS(fv)) 
     210    PRINT *, 'Max le_de : ', MAXVAL(ABS(le_de)) 
     211    PRINT *, 'Max Riv2 : ',  MAXVAL(ABS(Riv2)) 
     212    PRINT *, 'Max wee : ',   MAXVAL(ABS(wee)) 
     213    PRINT *, MINVAL(right),  MAXVAL(right) 
     214    PRINT *, MINVAL(right),  MAXVAL(left) 
    215215    PRINT *,' ... Done.' 
    216216    IF(nb_threads==0) nb_threads=OMP_GET_MAX_THREADS() 
Note: See TracChangeset for help on using the changeset viewer.