1. Nathan Collier
  2. NumPorCodes

Commits

Nathan Collier  committed 647b516

final work

  • Participants
  • Parent commits c307aa4
  • Branches master

Comments (0)

Files changed (3)

File DSW/DSW.c

View file
  • Ignore whitespace
 #define __FUNCT__ "Forcing"
 PetscReal Forcing(PetscScalar *x,PetscScalar t)
 {
-  if(t>10) return 0.;
-  PetscScalar r0 = 3000, r = sqrt(x[0]*x[0]+x[1]*x[1]);
-  if(r>r0) return 0.;
-  PetscScalar tmp = 1-r*r/r0/r0;
-  return 0.001*tmp*tmp;
+  PetscScalar r0 = 1000, r = sqrt(x[0]*x[0]+x[1]*x[1]);
+  PetscScalar t0 = 5.*60.*0.5; t = fabs(t-t0);
+  if(r>r0 || t>t0) return 1e-12;
+  PetscScalar space = 1-r*r/r0/r0;
+  PetscScalar time = 1-t*t/t0/t0;
+  return 2e-4*time*time*space*space+1e-12;
 }
 
 #undef  __FUNCT__
 
   TS ts;
   ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
-  ierr = TSSetDuration(ts,1000000,1000.0);CHKERRQ(ierr);
+  ierr = TSSetDuration(ts,1000000,40.*60);CHKERRQ(ierr);
   ierr = TSSetTimeStep(ts,1.0e-2);CHKERRQ(ierr);
   ierr = TSSetType(ts,TSALPHA);CHKERRQ(ierr);
   ierr = TSAlphaSetRadius(ts,0.5);CHKERRQ(ierr);

File DSW/solution/anime.py

View file
  • Ignore whitespace
-from __future__ import print_function
 import numpy as np
 import pylab as plt
 from igakit.cad import grid
 from igakit.io import PetIGA
 from igakit.nurbs import NURBS
-import glob,sys
-
-class ProgressBar:
-    def __init__(self, iterations):
-        self.iterations = iterations
-        self.prog_bar = '[]'
-        self.fill_char = '*'
-        self.width = 50
-        self.__update_amount(0)
-
-    def animate(self, iter):
-        print('\r', self, end='')
-        sys.stdout.flush()
-        self.update_iteration(iter + 1)
-
-    def update_iteration(self, elapsed_iter):
-        self.__update_amount((elapsed_iter / float(self.iterations)) * 100.0)
-        self.prog_bar += '  %d of %s complete' % (elapsed_iter, self.iterations)
-
-    def __update_amount(self, new_amount):
-        percent_done = int(round((new_amount / 100.0) * 100.0))
-        all_full = self.width - 2
-        num_hashes = int(round((percent_done / 100.0) * all_full))
-        self.prog_bar = '[' + self.fill_char * num_hashes + ' ' * (all_full - num_hashes) + ']'
-        pct_place = (len(self.prog_bar) // 2) - len(str(percent_done))
-        pct_string = '%d%%' % percent_done
-        self.prog_bar = self.prog_bar[0:pct_place] + \
-            (pct_string + self.prog_bar[pct_place + len(pct_string):])
-
-    def __str__(self):
-        return str(self.prog_bar)
+import glob,sys,os
+from progress import ProgressBar
 
 def Limits(datafiles,z):
-    print("Finding water depth limits...")
+    print "Finding water depth limits..."
     hmin = []
     hmax = []
-    p = ProgressBar(len(datafiles)-1)
     for i in range(len(datafiles)):
         fname = datafiles[i]
-        p.animate(i)
         u = PetIGA().read_vec(fname,iga)
         h = u-z
         hmin.append(h.min())
         hmax.append(h.max())
-    print("Limits found: 0, %.3f" % max(hmax))
+    print "Limits found: 0, %.3f" % max(hmax)
     return min(min(hmin),0),max(hmax)
 
+def ComputeU(t0,ts,iga):
+    i = ts.searchsorted(t0)
+    u0 = PetIGA().read_vec("./data/dsw%d.dat" % (i  ),iga)
+    uf = PetIGA().read_vec("./data/dsw%d.dat" % (i+1),iga)
+    a = (t0-ts[i])/(ts[i+1]-ts[i])
+    return u0+a*(uf-u0)
+
+def AdvectSensors(t0,tf,ts,dts,xs,iga):
+    a = ts.searchsorted(t0)
+    b = ts.searchsorted(tf)
+    ind = np.arange(xs.shape[0])
+    for i in range(a,b+1):
+        if i == a:
+            u  = ComputeU(t0,ts,iga)
+            dt = ts[a+1]-t0
+        elif i == b:
+            u = PetIGA().read_vec("./data/dsw%d.dat" % i,iga)
+            dt = tf-ts[b]
+        else:
+            u = PetIGA().read_vec("./data/dsw%d.dat" % i,iga)
+            dt = dts[i]
+        uu = iga.evaluate(fields=u,u=xs[:,0],v=xs[:,1])
+        du = iga.gradient(fields=u,u=xs[:,0],v=xs[:,1])
+        zz = iga.evaluate(fields=z,u=xs[:,0],v=xs[:,1])
+        magdu = np.sqrt(du[...,0]*du[...,0]+du[...,1]*du[...,1])
+        ka = (((uu-zz).clip(1e-15))**alpha * magdu**(gamma-1.))/Cf
+        vs = du*ka[...,None]
+        print "MAX: %e" % ((magdu*ka).max())
+        for j in range(vs.shape[0]): 
+            if uu[j,j]-zz[j,j] < 0.0: vs[j,j,:] = 0.
+        xs += vs[ind,ind,:]*dt
+    return xs
+
+os.system("rm time*.png")
+Ns = 7
+xs = np.zeros((Ns,2)) # x,y position
+xs[0,:] = [ 670, 1730]
+xs[1,:] = [ 732, 1278]
+xs[2,:] = [ 800, 878]
+xs[3,:] = [ 790, 411]
+xs[4,:] = [ 744,-100]
+xs[5,:] = [ 619,-348]
+xs[6,:] = [  46,-355]
+
+# Parse logfile
+t = []
+dts = []
+log = file('../dsw.log').readlines()
+for line in log:
+    line = line.split()
+    if len(line) < 2: continue
+    if line[1] == "TS":
+        t.append(float(line[-1]))
+        dts.append(float(line[3]))
+t = np.asarray(t)
+dts = np.asarray(dts)
+
+# initialize DSW stuff
 iga = PetIGA().read("../iga.dat")
-z   = PetIGA().read_vec("dsw0.dat",iga)
+z   = PetIGA().read_vec("./data/dsw0.dat",iga)
 alpha = 5./3.
 gamma = 0.5
-Cf    = 0.001
-
-datafiles = glob.glob("dsw*.dat")
-hmin,hmax = Limits(datafiles,z)
-#hmin,hmax = 0.0,0.536748612783
+Cf    = 0.009
+datafiles = glob.glob("./data/dsw*.dat")
+hmin,hmax = 0,0.25 # = Limits(datafiles,z)
+z[0,0] = -2.5 # keeps blue out of the terrain colormap
+
+# render plots
+print "Rendering plots..."
+tfinal = 2000.
+dt = 60.
+nsteps = int(tfinal/dt)
+#p = ProgressBar(nsteps-1)
+for i in range(nsteps):
+    time = i*dt
+    #p.animate(i)
 
-print("Rendering plots...")
-p = ProgressBar(len(datafiles))
-for i in range(len(datafiles)):
-    fname = datafiles[i]
-    p.animate(i)
-    u = PetIGA().read_vec(fname,iga)
-    h = u-z
     fig = plt.figure(figsize=(6,6),tight_layout=True)
-    cb = plt.imshow(h,vmin=hmin,vmax=hmax)
-    plt.colorbar(cb)    
-    CS = plt.contour(z,10,colors='k')
-    plt.clabel(CS, fontsize=10, inline=1, fmt='%d')
-    fig.savefig(fname.split(".")[0]+".png")
-    plt.close(fig)
-
-"""
-plt.figure(figsize=(15,7.5))
 
-for i in range(ith,ith+1):
-    fname = "./dsw%d.dat" % i
-    h = u-z
-    print "Reading %s..." % fname
-    print "Depth min/max = %e / %e" % (h.min(),h.max())
-    
-    x = (iga.breaks()[0])[:-1]+10.
-    uu = iga.evaluate(fields=u,u=x,v=x)
-    du = iga.gradient(fields=u,u=x,v=x)
-    zz = iga.evaluate(fields=z,u=x,v=x)
-    dz = iga.gradient(fields=z,u=x,v=x)
-
-    magdu = np.sqrt(du[...,0]*du[...,0]+du[...,1]*du[...,1])
-    ka = ((uu-zz).clip(0))**alpha * magdu**(gamma-1.)
-    vel = du*ka[...,None]
-
-    plt.subplot(121)
-
-    plt.subplot(122)
-    vel = vel[::5,::5,:]
-    z = z[::-1,:]
-    vel = vel[::-1,...]
-    X,Y = plt.meshgrid(iga.breaks(0),iga.breaks(1))
-    CS = plt.contour(X,Y,z,20,linewidths=2)
-    #plt.clabel(CS, fontsize=10, inline=1, fmt='%.1f')
-    X,Y = plt.meshgrid(x[::5],x[::5])
-    plt.quiver(X,Y,vel[...,1],vel[...,0],alpha=0.75)
-
-
-    plt.show()
-
-
-"""
+    # plot water height
+    u = ComputeU(time,t,iga)
+    h = (u-z)*100.
+    cb = plt.imshow(h.clip(0),vmax=100*hmax,cmap='Blues')
+    plt.colorbar(cb)
+
+    # plot terrain
+    CS = plt.contour(z,40,cmap='terrain',alpha=0.9)
+
+    # update and plot sensors
+    xs=AdvectSensors(max(0,(i-1)*dt),time,t,dts,xs,iga)
+    for j in range(Ns):
+        xt = ((xs[j,:]/4000.)+0.5)
+        xt *= h.shape[0]
+        plt.plot(xt[0],xt[1],'or',markersize=3,mew=0)
+
+    fig.axes[0].set_xticklabels([])
+    fig.axes[0].set_yticklabels([])
+    fig.axes[0].set_xticks([])
+    fig.axes[0].set_yticks([])
+    plt.title('Water Depth (cm)')
+    fig.savefig('time%d.png' % i)
+    plt.close(fig)

File DSW/solution/post.py

View file
  • Ignore whitespace
     pass
 
 iga = PetIGA().read("../iga.dat")
-z   = PetIGA().read_vec("dsw0.dat",iga)
+z   = PetIGA().read_vec("./data/dsw0.dat",iga)
 alpha = 5./3.
 gamma = 0.5
 Cf    = 0.009
 
-plt.figure(figsize=(15,7.5))
+plt.figure(figsize=(10,10))
 
 for i in range(ith,ith+1):
-    fname = "./dsw%d.dat" % i
+    fname = "./data/dsw%d.dat" % i
     u = PetIGA().read_vec(fname,iga)
     h = u-z
     print "Reading %s..." % fname
     print "Depth min/max = %e / %e" % (h.min(),h.max())
     
-    x = (iga.breaks()[0])[:-1]+10.
+    x  = iga.breaks(0)[:-1]+10.
     uu = iga.evaluate(fields=u,u=x,v=x)
     du = iga.gradient(fields=u,u=x,v=x)
     zz = iga.evaluate(fields=z,u=x,v=x)
     dz = iga.gradient(fields=z,u=x,v=x)
 
     magdu = np.sqrt(du[...,0]*du[...,0]+du[...,1]*du[...,1])
-    ka = ((uu-zz).clip(0))**alpha * magdu**(gamma-1.)
+    ka = ((uu-zz).clip(1e-15))**alpha * magdu**(gamma-1.)/Cf
     vel = du*ka[...,None]
 
-    plt.subplot(121)
-    cb=plt.imshow(h)
+    z[0,0] = -2.5 # keeps blue out of the terrain colormap, maybe use browns
+    CS = plt.contour(z,40,cmap='terrain',alpha=0.9)
+    #cb = plt.imshow(h.clip(0),cmap='Blues')
+    cb = plt.imshow(vel[...,0])
+    plt.colorbar(cb)
+    plt.title('Water Depth')
+    """
+    plt.subplot(122)
+    cb=plt.imshow(np.log10(ka))
     plt.colorbar(cb)    
     CS = plt.contour(z,10,colors='k')
     plt.clabel(CS, fontsize=10, inline=1, fmt='%d')
-
-    plt.subplot(122)
-    vel = vel[::5,::5,:]
-    z = z[::-1,:]
-    vel = vel[::-1,...]
-    X,Y = plt.meshgrid(iga.breaks(0),iga.breaks(1))
-    CS = plt.contour(X,Y,z,20,linewidths=2)
-    #plt.clabel(CS, fontsize=10, inline=1, fmt='%.1f')
-    X,Y = plt.meshgrid(x[::5],x[::5])
-    plt.quiver(X,Y,vel[...,1],vel[...,0],alpha=0.75)
-
-
+    plt.title('log10(Diffusivity)')
+    """
     plt.show()