← Back to branch summary

~fluidity-core/fluidity/sea-ice-branch

« back to all changes in this revision

Viewing changes to tests/mms_sediment/mms_sediment.py

  • Committer: Simon Mouradian
  • Date: 2012-10-19 10:35:59 UTC
  • mfrom: (3520.32.371 fluidity)
  • Revision ID: simon.mouradian06@imperial.ac.uk-20121019103559-y36qa47phc69q8sc
merge from trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/mms_sediment/mms_sediment.py
 
1
# This file was *autogenerated* from the file tests/mms_sediment/mms_sediment.sage.
 
2
from sage.all_cmdline import *   # import sage library
 
3
_sage_const_2 = Integer(2); _sage_const_2p39 = RealNumber('2.39'); _sage_const_0p65 = RealNumber('0.65'); _sage_const_0p66 = RealNumber('0.66'); _sage_const_10p0 = RealNumber('10.0'); _sage_const_0p02 = RealNumber('0.02'); _sage_const_0p01 = RealNumber('0.01'); _sage_const_0p07 = RealNumber('0.07'); _sage_const_1p0 = RealNumber('1.0'); _sage_const_1p1 = RealNumber('1.1'); _sage_const_1p3 = RealNumber('1.3'); _sage_const_1p4 = RealNumber('1.4'); _sage_const_1p7 = RealNumber('1.7'); _sage_const_0p1 = RealNumber('0.1'); _sage_const_0p0 = RealNumber('0.0'); _sage_const_0p3 = RealNumber('0.3'); _sage_const_0p2 = RealNumber('0.2'); _sage_const_0p5 = RealNumber('0.5'); _sage_const_0p4 = RealNumber('0.4'); _sage_const_0p6 = RealNumber('0.6'); _sage_const_0p9 = RealNumber('0.9'); _sage_const_3p0 = RealNumber('3.0'); _sage_const_2p3 = RealNumber('2.3'); _sage_const_2p1 = RealNumber('2.1'); _sage_const_2p0 = RealNumber('2.0'); _sage_const_5p0 = RealNumber('5.0'); _sage_const_0p75 = RealNumber('0.75'); _sage_const_0p288 = RealNumber('0.288'); _sage_const_0p33 = RealNumber('0.33'); _sage_const_1p625 = RealNumber('1.625'); _sage_const_1p3en7 = RealNumber('1.3e-7')
 
4
y = var('y')
 
5
 
 
6
def function(phi_0, phi_x, phi_y, phi_xy, 
 
7
             f_sin_x, f_cos_x, f_sin_y, f_cos_y, f_sin_xy, f_cos_xy, 
 
8
             alpha_x, alpha_y, alpha_xy):
 
9
    
 
10
    f_0 = phi_0 
 
11
    f_x = phi_x*(f_sin_x*sin(alpha_x*x) + f_cos_x*cos(alpha_x*x)) 
 
12
    f_y = phi_y*(f_sin_y*sin(alpha_y*y) + f_cos_y*cos(alpha_y*y)) 
 
13
    f_xy = phi_xy*(f_sin_xy*sin(alpha_xy*x*y/pi) + f_cos_xy*cos(alpha_xy*x*y/pi)) 
 
14
    f = f_0 + f_x + f_y + f_xy
 
15
    return f
 
16
 
 
17
u = function(_sage_const_0p5 , _sage_const_1p0 , _sage_const_0p3 , _sage_const_0p0 , 
 
18
             _sage_const_0p0 , _sage_const_1p0 , _sage_const_0p5 , _sage_const_0p5 , _sage_const_1p0 , _sage_const_0p0 ,
 
19
             _sage_const_0p9 , _sage_const_0p75 , _sage_const_1p0 )
 
20
p = function(-_sage_const_1p0 , _sage_const_1p0 , _sage_const_1p0 , _sage_const_1p0 ,
 
21
              _sage_const_1p0 , _sage_const_0p0 , _sage_const_0p0 , _sage_const_1p0 , _sage_const_1p0 , _sage_const_0p0 ,
 
22
              _sage_const_1p0 , _sage_const_1p0 , _sage_const_1p0 )
 
23
s1 = function(_sage_const_0p1 , -_sage_const_0p01 , -_sage_const_0p02 , _sage_const_0p07 , 
 
24
              _sage_const_1p0 , _sage_const_0p0 , _sage_const_0p0 , _sage_const_1p0 , _sage_const_1p0 , _sage_const_0p0 ,
 
25
              _sage_const_1p7 , _sage_const_2p1 , _sage_const_1p3 )
 
26
s2 = function(_sage_const_0p1 , _sage_const_0p01 , -_sage_const_0p01 , -_sage_const_0p1 , 
 
27
              _sage_const_1p0 , _sage_const_0p0 , _sage_const_0p0 , _sage_const_1p0 , _sage_const_1p0 , _sage_const_0p0 ,
 
28
              _sage_const_1p4 , _sage_const_3p0 , _sage_const_0p6 )
 
29
s1_d = function(_sage_const_10p0 , _sage_const_0p2 , _sage_const_0p4 , -_sage_const_0p2 , 
 
30
                _sage_const_1p0 , _sage_const_0p0 , _sage_const_0p0 , _sage_const_1p0 , _sage_const_1p0 , _sage_const_0p0 ,
 
31
                _sage_const_1p7 , _sage_const_1p1 , _sage_const_0p3 )
 
32
s2_d = function(_sage_const_5p0 , -_sage_const_0p1 , -_sage_const_0p2 , _sage_const_0p1 , 
 
33
                _sage_const_1p0 , _sage_const_0p0 , _sage_const_0p0 , _sage_const_1p0 , _sage_const_1p0 , _sage_const_0p0 ,
 
34
                _sage_const_1p4 , _sage_const_1p3 , _sage_const_2p3 )
 
35
v = integral(-diff(u,x),y)  # divergence free  
 
36
 
 
37
nu = _sage_const_1p0 /(_sage_const_1p0  - (s1 + s2)/_sage_const_0p65 )**_sage_const_1p625  # concentration dependent viscosity
 
38
 
 
39
s1_R = _sage_const_0p33 
 
40
s2_R = _sage_const_0p66 
 
41
rho = s1_R*s1 + s2_R*s2 
 
42
 
 
43
tau_xx = _sage_const_2 *nu*diff(u,x)            
 
44
tau_xy = nu*(diff(u,y) + diff(v,x))
 
45
tau_yy = _sage_const_2 *nu*diff(v,y)            
 
46
tau_yx = nu*(diff(u,y) + diff(v,x))  
 
47
 
 
48
s1_us = _sage_const_0p33 *(_sage_const_1p0  - (s1 + s2))**_sage_const_2p39  # hindered settling
 
49
s2_us = _sage_const_0p66 *(_sage_const_1p0  - (s1 + s2))**_sage_const_2p39  # hindered settling
 
50
 
 
51
Su = u*diff(u,x) + v*diff(u,y) - diff(tau_xx, x) - diff(tau_xy, y) - rho + diff(p,x)  
 
52
Sv = u*diff(v,x) + v*diff(v,y) - diff(tau_yx, x) - diff(tau_yy, y) - rho + diff(p,y)  
 
53
 
 
54
Ss1 = (u + s1_us)*diff(s1,x) + (v + s1_us)*diff(s1,y) - _sage_const_1p0 *(diff(s1, x, x) + diff(s1, y, y)) 
 
55
Ss2 = (u + s2_us)*diff(s2,x) + (v + s2_us)*diff(s2,y) - _sage_const_1p0 *(diff(s2, x, x) + diff(s2, y, y)) 
 
56
 
 
57
# EROSION RATE
 
58
 
 
59
s1_D = _sage_const_1p0 
 
60
s2_D = _sage_const_2p0 
 
61
 
 
62
R_p1 = ((s1_R*s1_D**_sage_const_3p0 )**_sage_const_0p5 )/nu
 
63
R_p2 = ((s2_R*s2_D**_sage_const_3p0 )**_sage_const_0p5 )/nu
 
64
 
 
65
mu = (s1_d*s1_D + s2_d*s2_D)/(s1_d + s2_d)
 
66
sigma = ((s1_d*(s1_D - mu)**_sage_const_2p0  + s2_d*(s2_D - mu)**_sage_const_2p0 )/(s1_d + s2_d))**_sage_const_0p5 
 
67
lambda_m = _sage_const_1p0  - _sage_const_0p288 *sigma
 
68
 
 
69
tau_b = u**_sage_const_2p0  + v**_sage_const_2p0     # drag shear stress // |cD*|u|*u| // cD = 1.0
 
70
 
 
71
d_50 = s1_D
 
72
 
 
73
Z_1 = (lambda_m*(tau_b**_sage_const_0p5 )*(R_p1**_sage_const_0p6 )*(s1_D/d_50)**_sage_const_0p2 )/s1_us
 
74
Z_2 = (lambda_m*(tau_b**_sage_const_0p5 )*(R_p2**_sage_const_0p6 )*(s2_D/d_50)**_sage_const_0p2 )/s2_us
 
75
 
 
76
F_1 = s1_d / (s1_d + s2_d)
 
77
F_2 = s2_d / (s1_d + s2_d)
 
78
 
 
79
A = _sage_const_1p3en7 
 
80
 
 
81
E_1 = s1_us*F_1*(A*Z_1**_sage_const_5p0 )/(_sage_const_1p0  + A*(Z_1**_sage_const_5p0 )/_sage_const_0p3 )
 
82
E_2 = s2_us*F_2*(A*Z_2**_sage_const_5p0 )/(_sage_const_1p0  + A*(Z_2**_sage_const_5p0 )/_sage_const_0p3 )
 
83
 
 
84
# DEPOSIT RATE
 
85
 
 
86
D_1_x = (u + s1_us)*s1
 
87
D_1_y = (v + s1_us)*s1
 
88
D_2_x = (u + s2_us)*s2
 
89
D_2_y = (v + s2_us)*s2
 
90
  
 
91
print 'from math import sin, cos, tanh, pi, e, sqrt'
 
92
print ''
 
93
print 'def u(X):'
 
94
print '    return', str(u).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
95
print ''
 
96
print 'def v(X):'
 
97
print '    return', str(v).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
98
print ''  
 
99
print 'def p(X):'
 
100
print '    return', str(p).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
101
print ''  
 
102
print 'def s1(X):'
 
103
print '    return', str(s1).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
104
print ''  
 
105
print 'def s2(X):'
 
106
print '    return', str(s2).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
107
print ''
 
108
print 'def s1_d(X):'
 
109
print '    return', str(s1_d).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
110
print ''  
 
111
print 'def s2_d(X):'
 
112
print '    return', str(s2_d).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
113
print ''
 
114
print 'def rho(X):'
 
115
print '    return', str(rho).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
116
print ''
 
117
print 'def nu(X):'
 
118
print '    return', str(nu).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
119
print ''
 
120
print 'def mu(X):'
 
121
print '    return', str(mu).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
122
print ''
 
123
print 'def sigma(X):'
 
124
print '    return', str(sigma).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
125
print ''
 
126
print 'def tau_b(X):'
 
127
print '    return', str(tau_b).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
128
print ''
 
129
print 'def E_1(X):'
 
130
print '    if X[1] < 1e-6:'
 
131
print '        return', str(E_1).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
132
print '    else:'
 
133
print '        return 0.0'
 
134
print ''
 
135
print 'def E_2(X):'
 
136
print '    if X[1] < 1e-6:'
 
137
print '        return', str(E_2).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
138
print '    else:'
 
139
print '        return 0.0'
 
140
print ''
 
141
print 'def D_1(X):'
 
142
print '    if X[1] < 1e-6:'
 
143
print '        return', str(-D_1_y).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
144
print '    else:'
 
145
print '        return 0.0'
 
146
print ''
 
147
print 'def D_2(X):'
 
148
print '    if X[1] < 1e-6:'
 
149
print '        return', str(-D_2_y).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
150
print '    else:'
 
151
print '        return 0.0'
 
152
print ''
 
153
print 'def forcing_u(X):'
 
154
print '    return', str(Su).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
155
print ''
 
156
print 'def forcing_v(X):'
 
157
print '    return', str(Sv).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
158
print ''
 
159
print 'def forcing_s1(X):'
 
160
print '    return', str(Ss1).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
161
print ''
 
162
print 'def forcing_s2(X):'
 
163
print '    return', str(Ss2).replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
164
print ''
 
165
print 'def velocity(X):'
 
166
print '    return [u(X), v(X)]'
 
167
print ''
 
168
print 'def forcing_velocity(X):'
 
169
print '    return [forcing_u(X), forcing_v(X)]'