← Back to branch summary

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

« back to all changes in this revision

Viewing changes to tests/mphase_mms_p2p1_stress_form/mphase_mms_p2p1_stress_form.sage

  • 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/mphase_mms_p2p1_stress_form/mphase_mms_p2p1_stress_form.sage
 
1
y = var('y')
 
2
 
 
3
def function(phi_0, phi_x, phi_y, phi_xy, 
 
4
             f_sin_x, f_cos_x, f_sin_y, f_cos_y, f_sin_xy, f_cos_xy, 
 
5
             alpha_x, alpha_y, alpha_xy):
 
6
    
 
7
    f_0 = phi_0 
 
8
    f_x = phi_x*(f_sin_x*sin(alpha_x*x) + f_cos_x*cos(alpha_x*x)) 
 
9
    f_y = phi_y*(f_sin_y*sin(alpha_y*y) + f_cos_y*cos(alpha_y*y)) 
 
10
    f_xy = phi_xy*(f_sin_xy*sin(alpha_xy*x*y/pi) + f_cos_xy*cos(alpha_xy*x*y/pi)) 
 
11
    f = f_0 + f_x + f_y + f_xy
 
12
    return f
 
13
 
 
14
p = function(-1.0, 1.0, 1.0, 1.0,
 
15
             1.0, 0.0, 0.0, 1.0, 1.0, 0.0,
 
16
             1.0, 1.0, 1.0)
 
17
rho1 = 2.5
 
18
rho2 = 0.5
 
19
u1 = 0.25*cos(x)*cos(y) - x*cos(y)
 
20
v1 = sin(y)
 
21
u2 = sin(x)*cos(y)
 
22
v2 = sin(y)*sin(x) - cos(x)*sin(y)
 
23
             
 
24
vfrac1 = 0.8 #0.1 + y/20.0
 
25
vfrac2 = 1.0 - vfrac1 #1.0 - vfrac1
 
26
 
 
27
#print "DIVERGENCE = ", diff(vfrac1*u1,x) + diff(vfrac1*v1,y) + diff(vfrac2*u2,x) + diff(vfrac2*v2,y)
 
28
 
 
29
nu = 0.7
 
30
 
 
31
tau_xx1 = 2*nu*diff(u1,x) - (2.0/3.0)*nu*(diff(u1,x) + diff(v1,y))
 
32
tau_xy1 = nu*(diff(u1,y) + diff(v1,x))
 
33
tau_yy1 = 2*nu*diff(v1,y) - (2.0/3.0)*nu*(diff(u1,x) + diff(v1,y)) 
 
34
tau_yx1 = nu*(diff(u1,y) + diff(v1,x))
 
35
 
 
36
tau_xx2 = 2*nu*diff(u2,x) - (2.0/3.0)*nu*(diff(u2,x) + diff(v2,y))
 
37
tau_xy2 = nu*(diff(u2,y) + diff(v2,x))
 
38
tau_yy2 = 2*nu*diff(v2,y) - (2.0/3.0)*nu*(diff(u2,x) + diff(v2,y)) 
 
39
tau_yx2 = nu*(diff(u2,y) + diff(v2,x))
 
40
 
 
41
Su1 = vfrac1*rho1*u1*diff(u1,x) + vfrac1*rho1*v1*diff(u1,y) - diff(vfrac1*tau_xx1, x) - diff(vfrac1*tau_xy1, y) + vfrac1*diff(p,x)  
 
42
Sv1 = vfrac1*rho1*u1*diff(v1,x) + vfrac1*rho1*v1*diff(v1,y) - diff(vfrac1*tau_yx1, x) - diff(vfrac1*tau_yy1, y) + vfrac1*diff(p,y)  
 
43
Su2 = vfrac2*rho2*u2*diff(u2,x) + vfrac2*rho2*v2*diff(u2,y) - diff(vfrac2*tau_xx2, x) - diff(vfrac2*tau_xy2, y) + vfrac2*diff(p,x)  
 
44
Sv2 = vfrac2*rho2*u2*diff(v2,x) + vfrac2*rho2*v2*diff(v2,y) - diff(vfrac2*tau_yx2, x) - diff(vfrac2*tau_yy2, y) + vfrac2*diff(p,y)  
 
45
  
 
46
print 'from math import sin, cos, tanh, pi'
 
47
print ''
 
48
print 'def u1(X):'
 
49
print '    return', str(u1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
50
print ''
 
51
print 'def v1(X):'
 
52
print '    return', str(v1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
53
print ''  
 
54
print 'def u2(X):'
 
55
print '    return', str(u2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
56
print ''
 
57
print 'def v2(X):'
 
58
print '    return', str(v2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
59
print ''  
 
60
print 'def vfrac1(X):'
 
61
print '    return', str(vfrac1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
62
print ''
 
63
print 'def vfrac2(X):'
 
64
print '    return', str(vfrac2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
65
print ''  
 
66
print 'def p(X):'
 
67
print '    return', str(p).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
68
print '' 
 
69
print 'def rho1(X):'
 
70
print '    return', str(rho1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
71
print ''
 
72
print 'def rho2(X):'
 
73
print '    return', str(rho2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
74
print '' 
 
75
print 'def forcing_u1(X):'
 
76
print '    return', str(Su1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
77
print ''
 
78
print 'def forcing_v1(X):'
 
79
print '    return', str(Sv1).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
80
print ''
 
81
print 'def forcing_u2(X):'
 
82
print '    return', str(Su2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
83
print ''
 
84
print 'def forcing_v2(X):'
 
85
print '    return', str(Sv2).replace('e^', 'exp').replace('^', '**').replace('000000000000', '').replace('x', 'X[0]').replace('y', 'X[1]')
 
86
print ''
 
87
print 'def velocity1(X):'
 
88
print '   return [u1(X), v1(X)]'
 
89
print ''
 
90
print 'def velocity2(X):'
 
91
print '   return [u2(X), v2(X)]'
 
92
print ''
 
93
print 'def forcing_velocity1(X):'
 
94
print '   return [forcing_u1(X), forcing_v1(X)]'
 
95
print ''
 
96
print 'def forcing_velocity2(X):'
 
97
print '   return [forcing_u2(X), forcing_v2(X)]'
 
98
print ''