Changeset 8466 for trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

Timestamp:

Jul 11, 2012, 9:42:10 PM (13 years ago)

Author:

davies

Message:

Adding new 'okada_tsunami' routines, which are well tested
and can treat both simple and complex faults

File:

• trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c (modified) (6 diffs)

trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

@@ -493 +493 @@
                 //Bedslope approx 1:
                 //bedslope_work = g*hc*(ql[0])*length-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length;
+                bedslope_work = g*length*( hc*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
                 //
                 // Bedslope approx 2
@@ -505 +506 @@
 
                 // Bedslope approx 3
-                bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
+                //bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
                 //
                 xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
@@ -535 +536 @@
                     // Bedslope approx 1:
                     //bedslope_work = g*hc_n*(qr[0])*length-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
+                    bedslope_work = g*length*(hc_n*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
                     //
                     // Bedslope approx 2:
@@ -547 +549 @@
                     //
                     // Bedslope approx 3
-                    bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
+                    //bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
                     //
                     xmom_explicit_update[n] += normals[ki2]*bedslope_work;
@@ -1079 +1081 @@
       dqv[1] = a*dxv1 + b*dyv1;
       dqv[2] = a*dxv2 + b*dyv2;
-      
+    
+      // Idea: New limiter, computed using only neighbouring centroid values and local centroid value 
+      // Step1: Compute the value of stage at the centroid of triangle k ,
+      //       using only the neighbouring centroid stage values
+      //tmp = a*(x-x0) + b*(y-y0) + stage_centroid_values[k0];
+      // Step2: Now, compute a limiting factor, so that if gradients are multiplied by this,
+      //        then for a triangle with these limited gradients, based at the local centroid value
+      //        ,the re-constructed neighbour centroid values will not over shoot (if larger) or undershoot (if smaller)
+      //limiting_factor=1 - max_all_k0((tmp - stage_centroid_values[k])/(tmp - stage_centroid_values[k0]))
+      //limiting_factor=min(max(limiting_factor,0),1)
+      // Advantages: No limiting for linear problem. Worth a look anyway
+ 
       // Now we want to find min and max of the centroid and the 
       // vertices of the auxiliary triangle and compute jumps 
@@ -1088 +1101 @@
     
       // Limit the gradient
+      //if(hmin/hc<0.5){
       limit_gradient(dqv, qmin, qmax, beta_tmp); 
+      //}
       //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);