1. Nathan Collier
  2. NumPorCodes

Commits

Nathan Collier  committed c307aa4

added DSW work

  • Participants
  • Parent commits 7aa6b25
  • Branches master

Comments (0)

Files changed (5)

File DSW/DSW.c

View file
+#include "petiga.h"
+
+PetscScalar TOL_GRADU = 1e-12, TOL_FO = 0.95;
+
+typedef struct { 
+  IGA iga;
+  PetscScalar m,gamma,alpha,Cf,kappa_max,h2;
+} AppCtx;
+
+#undef __FUNCT__
+#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;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Diffusivity"
+void Diffusivity(PetscScalar u,PetscScalar *u_x,PetscScalar z,PetscScalar *kappa,PetscScalar *kappa_u,AppCtx *user)
+{
+  PetscReal du = u_x[0]*u_x[0]+u_x[1]*u_x[1];
+  PetscReal alpha = user->alpha, gamma = user->gamma, Cf = user->Cf;
+  PetscBool ok = (u >= z && sqrt(du) > TOL_GRADU);
+  if (kappa) {
+    *kappa = 0;
+    if(ok) *kappa = pow(u-z,alpha)*pow(du,0.5*(gamma-1.))/Cf;
+  }
+  if (kappa_u) {
+    kappa_u[0] = 0;
+    kappa_u[1] = 0;
+    if(ok){
+      kappa_u[0] = pow(u-z,alpha)*(gamma-1.)*pow(du,0.5*(gamma-1.)-1.)/Cf;
+      kappa_u[1] = alpha*pow(u-z,alpha-1.)*pow(du,0.5*(gamma-1.))/Cf;
+    }
+  }
+  return;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Residual"
+PetscErrorCode Residual(IGAPoint p,PetscReal dt,
+                        PetscReal shift,const PetscScalar *V,
+                        PetscReal t,const PetscScalar *U,
+                        PetscScalar *R,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+
+  PetscInt nen;
+  IGAPointGetSizes(p,0,&nen,0);
+
+  PetscScalar u_t,u,u_x[2],z;
+  IGAPointFormValue(p,V,&u_t);
+  IGAPointFormValue(p,U,&u);
+  IGAPointFormGrad (p,U,&u_x[0]);
+  IGAPointFormValue(p,p->property,&z);
+
+  PetscScalar kappa,f;
+  Diffusivity(u,u_x,z,&kappa,NULL,user);
+  f=Forcing(p->point,t);
+
+  user->kappa_max = PetscMax(kappa,user->kappa_max);
+
+  const PetscReal *N0,(*N1)[2];
+  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+
+  PetscInt a;
+  for (a=0; a<nen; a++) {
+    PetscReal Na    = N0[a];
+    PetscReal Na_x  = N1[a][0];
+    PetscReal Na_y  = N1[a][1];
+    R[a] = Na*u_t + kappa*(Na_x*u_x[0]+Na_y*u_x[1]) - Na*f;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Jacobian"
+PetscErrorCode Jacobian(IGAPoint p,PetscReal dt,
+			PetscReal shift,const PetscScalar *V,
+			PetscReal t,const PetscScalar *U,
+			PetscScalar *J,void *ctx)
+{
+ AppCtx *user = (AppCtx *)ctx;
+
+  PetscInt nen;
+  IGAPointGetSizes(p,0,&nen,0);
+
+  PetscScalar u_t,u,u_x[2],z;
+  IGAPointFormValue(p,V,&u_t);
+  IGAPointFormValue(p,U,&u);
+  IGAPointFormGrad (p,U,&u_x[0]);
+  IGAPointFormValue(p,p->property,&z);
+
+  PetscScalar kappa,kappa_u[2];
+  Diffusivity(u,u_x,z,&kappa,kappa_u,user);
+
+  const PetscReal *N0,(*N1)[2];
+  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+
+  PetscInt a,b;
+  for (a=0; a<nen; a++) {
+    PetscReal Na    = N0[a];
+    PetscReal Na_x  = N1[a][0];
+    PetscReal Na_y  = N1[a][1];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb    = N0[b];
+      PetscReal Nb_x  = N1[b][0];
+      PetscReal Nb_y  = N1[b][1];
+      J[b*nen+a]  = shift*Na*Nb;
+      J[b*nen+a] += kappa*(Na_x*Nb_x+Na_y*Nb_y);
+      J[b*nen+a] += kappa_u[0]*(Na_x*u_x[0]+Na_y*u_x[1]);
+      J[b*nen+a] += kappa_u[1]*Nb;
+    }
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "OutputMonitor"
+PetscErrorCode OutputMonitor(TS ts,PetscInt it_number,PetscReal c_time,Vec U,void *mctx)
+{
+  PetscFunctionBegin;
+  PetscErrorCode ierr;
+  AppCtx *user = (AppCtx *)mctx;
+  char           filename[256];
+  sprintf(filename,"./solution/dsw%d.dat",it_number);
+  ierr = IGAWriteVec(user->iga,U,filename);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "FormInitialCondition"
+PetscErrorCode FormInitialCondition(IGA iga,Vec U)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  Vec          Ul;
+  PetscScalar *u;
+  PetscInt     n;
+  ierr = IGAGetLocalVec(iga,&Ul);CHKERRQ(ierr);
+  ierr = VecGetArray(Ul,&u);CHKERRQ(ierr);
+  ierr = VecGetSize(Ul,&n);CHKERRQ(ierr);
+  ierr = PetscMemcpy(u,iga->propertyA,n*sizeof(PetscScalar));CHKERRQ(ierr);
+  ierr = VecRestoreArray(Ul,&u);CHKERRQ(ierr);
+  ierr = IGALocalToGlobal(iga,Ul,U,INSERT_VALUES);CHKERRQ(ierr);
+  ierr = VecDestroy(&Ul);CHKERRQ(ierr);
+  PetscFunctionReturn(0); 
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "KISSAdaptivity"
+PetscErrorCode KISSAdaptivity(TS ts,PetscReal t,Vec X,Vec Xdot, PetscReal *nextdt,PetscBool *ok,void *ctx)
+{
+  PetscErrorCode       ierr;
+  PetscReal            dt;
+  SNES                 snes;
+  SNESConvergedReason  snesreason;
+  PetscInt             sit;
+  ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
+  ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
+  ierr = SNESGetConvergedReason(snes,&snesreason);CHKERRQ(ierr);
+  ierr = SNESGetIterationNumber(snes,&sit);CHKERRQ(ierr);
+
+  // default
+  *ok = PETSC_TRUE;
+  *nextdt = dt;
+  
+  // Check the worst-case Fourier number
+  AppCtx     *user = (AppCtx *)ctx;
+  MPI_Comm    comm;
+  PetscScalar kappa_max,Fo;
+  ierr = IGAGetComm(user->iga,&comm);CHKERRQ(ierr);
+  ierr = MPI_Allreduce(&(user->kappa_max),&kappa_max,1,MPIU_REAL,MPIU_MAX,comm);CHKERRQ(ierr);
+  user->kappa_max = 0;
+  Fo = kappa_max*dt/user->h2;
+
+  if (Fo > TOL_FO) { // Fo number too large
+    *ok = PETSC_FALSE;
+    *nextdt = TOL_FO*user->h2/kappa_max;
+    *nextdt = PetscMin(0.1*dt,*nextdt);
+    PetscFunctionReturn(0);
+  }
+  if (snesreason < 0) { // SNES failed
+    *ok = PETSC_FALSE;
+    *nextdt = 0.1*dt;
+    PetscFunctionReturn(0);
+  }
+  if (sit < 3) { // too easy, increase dt
+    *ok = PETSC_TRUE;
+    *nextdt = 1.1*dt;
+    PetscFunctionReturn(0);
+  }
+  if (sit > 5) { // too hard, decrease dt
+    *ok = PETSC_TRUE;
+    *nextdt = 0.5*dt;
+    PetscFunctionReturn(0);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  AppCtx user;
+  user.alpha     = 5.0/3.0;
+  user.gamma     = 0.5;
+  user.Cf        = 0.009;
+  user.m         = 1.0+user.alpha/user.gamma;
+  user.kappa_max = 0;
+  PetscScalar L = 4e3, dx = 10;
+  PetscInt p = 1,k = 0,N = (PetscInt)(L/dx);
+  PetscBool output = PETSC_TRUE;
+  user.h2 = dx*dx;
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,2);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+
+  IGAAxis axis0;
+  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
+  ierr = IGAAxisSetDegree(axis0,p);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,N,-0.5*L,0.5*L,k);CHKERRQ(ierr);
+  IGAAxis axis;
+  ierr = IGAGetAxis(iga,1,&axis);CHKERRQ(ierr);
+  ierr = IGAAxisCopy(axis0,axis);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+  ierr = IGAWrite(iga,"iga.dat");CHKERRQ(ierr);
+  user.iga = iga;
+  
+  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetUserIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
+
+  TS ts;
+  ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
+  ierr = TSSetDuration(ts,1000000,1000.0);CHKERRQ(ierr);
+  ierr = TSSetTimeStep(ts,1.0e-2);CHKERRQ(ierr);
+  ierr = TSSetType(ts,TSALPHA);CHKERRQ(ierr);
+  ierr = TSAlphaSetRadius(ts,0.5);CHKERRQ(ierr);
+  ierr = TSAlphaSetAdapt(ts,KISSAdaptivity,&user);CHKERRQ(ierr); 
+  if (output)  {ierr = TSMonitorSet(ts,OutputMonitor,&user,NULL);CHKERRQ(ierr);}
+  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+
+  Vec U;
+  ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
+  ierr = FormInitialCondition(iga,U);CHKERRQ(ierr);
+  ierr = TSSolve(ts,U);CHKERRQ(ierr);
+
+  ierr = VecDestroy(&U);CHKERRQ(ierr);
+  ierr = TSDestroy(&ts);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

File DSW/makefile

View file
+TARGETS = DSW
+
+ALL: ${TARGETS}
+clean::
+	-@${RM} ${TARGETS}
+
+include ${PETIGA_DIR}/conf/petigavariables
+include ${PETIGA_DIR}/conf/petigarules
+
+DSW: DSW.o
+	${CLINKER} -o $@ $^ ${PETIGA_LIB}
+	${RM} -f $<

File DSW/solution/anime.py

View file
+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)
+
+def Limits(datafiles,z):
+    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))
+    return min(min(hmin),0),max(hmax)
+
+iga = PetIGA().read("../iga.dat")
+z   = PetIGA().read_vec("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
+
+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()
+
+
+"""

File DSW/solution/post.py

View file
+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
+
+ith = 1
+try:
+    ith = int(sys.argv[1])
+except:
+    pass
+
+iga = PetIGA().read("../iga.dat")
+z   = PetIGA().read_vec("dsw0.dat",iga)
+alpha = 5./3.
+gamma = 0.5
+Cf    = 0.009
+
+plt.figure(figsize=(15,7.5))
+
+for i in range(ith,ith+1):
+    fname = "./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.
+    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)
+    cb=plt.imshow(h)
+    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.show()
+
+

File DSW/terrain/extract.py

View file
+import numpy as np
+import pylab as plt
+from igakit.cad import grid
+from igakit.io import PetIGA
+from igakit.nurbs import NURBS
+
+def LaplacianSmoothing(A,it):
+    for i in range(it):
+        B = np.zeros(A.shape)
+        B[:-1,:] += A[1:,:]
+        B[1:,:]  += A[:-1,:]
+        B[:,:-1] += A[:,1:]
+        B[:,1:]  += A[:,:-1]
+        B[1:-1,1:-1] /= 4.
+        B[1:-1,0] /= 3.
+        B[1:-1,-1] /= 3.
+        B[0,1:-1] /= 3.
+        B[-1,1:-1] /= 3.
+        B[0,0] /= 2.
+        B[-1,0] /= 2.
+        B[0,-1] /= 2.
+        B[-1,-1] /= 2.
+        A = np.copy(B)
+    return A
+
+data = np.fromfile('./n37w119/floatn37w119_1.flt',dtype='f4').reshape((3612,-1))
+clip = np.copy(data[1200:2000,1600:2400])[::2,::2]
+
+plt.subplot(121)
+cb=plt.imshow(clip)
+plt.colorbar(cb)
+clip = LaplacianSmoothing(clip,100)
+x0 = clip.min()
+dx = clip.max()-clip.min()
+clip = (clip-x0)/dx *10.
+plt.subplot(122)
+cb=plt.imshow(clip)
+plt.colorbar(cb)
+plt.show()
+
+z = grid((clip.shape[0]-1,clip.shape[1]-1),degree=1,continuity=0,limits=[-2e3,2e3])
+nrbs = NURBS(z.knots,z.control,fields=clip)
+PetIGA().write('../elevation.dat',nrbs,control=False,fields=True,nsd=2)