Changeset 221


Ignore:
Timestamp:
Aug 25, 2004, 5:26:37 PM (21 years ago)
Author:
ole
Message:

Painful testing of second order extrapolation.
One bug found: The protection against infinitesimal heighs at centroid. This was not done in the previous version
and was consequently removed.

Location:
inundation/ga/storm_surge/pyvolution
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/config.py

    r205 r221  
    5959#use_extensions = False   #Do not use C-extensions
    6060
    61 use_psyco = True  #Use psyco optimisations
    62 #use_psyco = False  #Do not use psyco optimisations
     61#use_psyco = True  #Use psyco optimisations
     62use_psyco = False  #Do not use psyco optimisations
    6363
    6464
  • inundation/ga/storm_surge/pyvolution/quantity.py

    r215 r221  
    7878            v2 = self.vertex_values[i, 2]
    7979           
    80             self.centroid_values[i] = (v0 + v1 + v2)/3
     80            self.centroid_values[i] = (v0 + v1 + v2)/3.0
    8181           
    8282            self.edge_values[i, 0] = 0.5*(v1 + v2)
     
    139139                #Special case where we have only one neighbouring volume.
    140140
    141                 k0 = k  #Self
    142141                #Find index of the one neighbour
    143                 for k1 in self.domain.neighbours[k,:]:
    144                     if k1 >= 0:
     142                for k0 in self.domain.neighbours[k,:]:
     143                    if k0 >= 0:
    145144                        break
    146                 assert k1 != k0
    147                 assert k1 >= 0
     145
     146                assert k0 != k
     147                assert k0 >= 0
     148
     149                k1 = k  #Self
    148150
    149151                #Get data
     
    180182                x2, y2 = self.domain.centroids[k2] #V2 centroid
    181183
     184                #if k < 4:
     185                #    print 'k, k0, k1, k2 = %d, %d, %d, %d' %(k, k0, k1, k2)
     186                #    #print 'k = %d: (%f, %f), (%f, %f), (%f, %f)'\
     187                #    #      %(k, x0, y0, x1, y1, x2, y2)
     188                #    print 'k = %d: (%f, %f, %f)' %(k, q0, q1, q2)   
     189               
    182190                #Gradient
    183191                a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
     
    266274       
    267275        a, b = self.compute_gradients()
     276        ##for k in range(4):
     277        ##    print 'Gradients %d: %16.12f, %16.12f' %(k, a[k], b[k])
    268278
    269279        V = self.domain.get_vertex_coordinates()
  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r215 r221  
    272272
    273273def distribute_to_vertices_and_edges(domain):
    274    
    275     protect_against_infinitesimal_heights_centroid(domain)
     274
     275    #print 'Calling distrib within SW'
    276276    if domain.order == 1:
     277        protect_against_infinitesimal_heights_centroid(domain)       
    277278        extrapolate_first_order(domain)
    278279    elif domain.order == 2:
     280        protect_against_negative_heights_centroid(domain)       
    279281        extrapolate_second_order(domain)       
    280282    else:
     
    291293    minimum allowed height
    292294    """
     295    #FIXME: Used only in first order
    293296   
    294297    #Water levels at centroids
     
    316319               
    317320
     321def protect_against_negative_heights_centroid(domain):
     322    """Adjust height and momentum at centroid if height is less than zero
     323    """
     324
     325    #FIXME: USed only in second order
     326    #Water levels at centroids
     327    wc = domain.quantities['level'].centroid_values
     328
     329    #Bed elevations at centroids
     330    zc = domain.quantities['elevation'].centroid_values
     331
     332    #Water depths at centroids   
     333    hc = wc - zc
     334
     335    #Momentums at centroids
     336    xmomc = domain.quantities['xmomentum'].centroid_values
     337    ymomc = domain.quantities['ymomentum'].centroid_values       
     338
     339    for k in range(domain.number_of_elements):
     340        #Protect against infinitesimal heights and high velocities
     341        if hc[k] < 0.0:
     342            #Control level and height
     343            wc[k] = zc[k]; hc[k] = 0.0
     344
     345            #Control momentum
     346            xmomc[k] = ymomc[k] = 0.0
     347               
    318348
    319349
  • inundation/ga/storm_surge/pyvolution/test_shallow_water.py

    r218 r221  
    960960        #assert allclose(domain.max_timestep, 0.0396825396825)
    961961        #print domain.quantities['level'].centroid_values
     962
     963
     964    def test_flatbed_first_order(self):
     965        from mesh_factory import rectangular
     966        from shallow_water import Domain,\
     967             Reflective_boundary, Dirichlet_boundary
     968
     969        from Numeric import array
     970
     971        #Create basic mesh
     972        N = 8
     973        points, vertices, boundary = rectangular(N, N)
     974
     975        #Create shallow water domain
     976        domain = Domain(points, vertices, boundary)
     977        domain.smooth = False
     978        domain.visualise = False
     979        domain.default_order=1
     980
     981        # Boundary conditions
     982        Br = Reflective_boundary(domain)
     983        Bd = Dirichlet_boundary([0.2,0.,0.])
     984       
     985        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
     986        domain.check_integrity()
     987
     988 
     989        #Evolution
     990        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
     991            pass
     992            #domain.write_time()
     993
     994        assert allclose(domain.min_timestep, 0.0140413643926)
     995        assert allclose(domain.max_timestep, 0.0140947355753)
     996
     997        for i in range(3):
     998            assert allclose(domain.quantities['level'].edge_values[:4,i],
     999                            [0.10730244,0.12337617,0.11200126,0.12605666])
     1000
     1001            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
     1002                            [0.07610894,0.06901572,0.07284461,0.06819712])
     1003
     1004            assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
     1005                            [-0.0060238, -0.00157404, -0.00309633, -0.0001637])                         
     1006       
     1007    def test_flatbed_second_order(self):
     1008        from mesh_factory import rectangular
     1009        from shallow_water import Domain,\
     1010             Reflective_boundary, Dirichlet_boundary
     1011
     1012        from Numeric import array
     1013
     1014        #Create basic mesh
     1015        N = 8
     1016        points, vertices, boundary = rectangular(N, N)
     1017
     1018        #Create shallow water domain
     1019        domain = Domain(points, vertices, boundary)
     1020        domain.smooth = False
     1021        domain.visualise = False
     1022        domain.default_order=2
     1023
     1024        # Boundary conditions
     1025        Br = Reflective_boundary(domain)
     1026        Bd = Dirichlet_boundary([0.2,0.,0.])
     1027
     1028        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
     1029        domain.check_integrity()
     1030 
     1031        #Evolution
     1032        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
     1033            pass
     1034
     1035
     1036        assert allclose(domain.min_timestep, 0.0210448446782)
     1037        assert allclose(domain.max_timestep, 0.0210448446782)
     1038
     1039        #print domain.quantities['level'].vertex_values[:4,0]
     1040        #print domain.quantities['xmomentum'].vertex_values[:4,0]
     1041        #print domain.quantities['ymomentum'].vertex_values[:4,0]
     1042       
     1043        assert allclose(domain.quantities['level'].vertex_values[:4,0],
     1044                        [0.00101913,0.05352143,0.00104852,0.05354394])
     1045        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
     1046                        [0.00090581,0.03685719,0.00088303,0.03687313])
     1047        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
     1048                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
     1049
     1050
     1051       
     1052    def test_flatbed_second_order_distribute(self):
     1053        #Use real data from pyvolution 2
     1054        #painfully setup and extracted.
     1055        from mesh_factory import rectangular
     1056        from shallow_water import Domain,\
     1057             Reflective_boundary, Dirichlet_boundary
     1058
     1059        from Numeric import array
     1060
     1061        #Create basic mesh
     1062        N = 8
     1063        points, vertices, boundary = rectangular(N, N)
     1064
     1065        #Create shallow water domain
     1066        domain = Domain(points, vertices, boundary)
     1067        domain.smooth = False
     1068        domain.visualise = False
     1069        domain.default_order=domain.order=2
     1070
     1071        # Boundary conditions
     1072        Br = Reflective_boundary(domain)
     1073        Bd = Dirichlet_boundary([0.2,0.,0.])
     1074
     1075        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
     1076        domain.check_integrity()
     1077
     1078
     1079
     1080        for V in [False, True]:
     1081            if V:
     1082                #Set centroids as if system had been evolved           
     1083                L = zeros(2*N*N, Float)
     1084                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
     1085                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
     1086                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
     1087                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
     1088                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
     1089                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
     1090                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
     1091                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
     1092                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
     1093                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
     1094                          0.00000000e+000, 5.57305948e-005]
     1095       
     1096                X = zeros(2*N*N, Float)
     1097                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
     1098                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
     1099                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
     1100                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
     1101                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
     1102                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
     1103                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
     1104                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
     1105                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
     1106                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
     1107                          0.00000000e+000, 4.57662812e-005]
     1108
     1109                Y = zeros(2*N*N, Float)
     1110                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
     1111                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
     1112                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
     1113                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
     1114                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
     1115                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
     1116                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
     1117                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
     1118                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
     1119                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
     1120                        0.00000000e+000, -2.57635178e-005]
     1121
     1122       
     1123                domain.set_quantity('level', L, 'centroids')
     1124                domain.set_quantity('xmomentum', X, 'centroids')
     1125                domain.set_quantity('ymomentum', Y, 'centroids')
     1126       
     1127                domain.check_integrity()
     1128            else:   
     1129                #Evolution
     1130                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
     1131                    pass
     1132                assert allclose(domain.min_timestep, 0.0210448446782)
     1133                assert allclose(domain.max_timestep, 0.0210448446782)
     1134       
     1135
     1136            #Centroids were correct but not vertices.
     1137            #Hence the check of distribute below.
     1138            assert allclose(domain.quantities['level'].centroid_values[:4],
     1139                            [0.00721206,0.05352143,0.01009108,0.05354394])
     1140
     1141            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
     1142                            [0.00648352,0.03685719,0.00850733,0.03687313])
     1143
     1144            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
     1145                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
     1146
     1147            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
     1148            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
     1149
     1150            assert allclose(domain.quantities['xmomentum'].centroid_values[17],
     1151                            0.00028606084)
     1152
     1153            import copy
     1154            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
     1155            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
     1156       
     1157            domain.distribute_to_vertices_and_edges()
     1158
     1159            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
     1160
     1161            assert allclose(domain.quantities['xmomentum'].centroid_values[17],
     1162                            0.00028606084)
     1163
     1164
     1165            assert allclose(domain.quantities['level'].vertex_values[:4,0],
     1166                            [0.00101913,0.05352143,0.00104852,0.05354394])
     1167
     1168            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
     1169                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
     1170
     1171            #This was the culprit. First triangles vertex 0 had an
     1172            #x-momentum of 0.0064835 instead of 0.00090581 and
     1173            #third triangle had 0.00850733 instead of 0.00088303
     1174            #print domain.quantities['xmomentum'].vertex_values[:4,0]
     1175       
     1176            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
     1177                            [0.00090581,0.03685719,0.00088303,0.03687313])       
     1178
     1179
    9621180
    9631181       
Note: See TracChangeset for help on using the changeset viewer.