/[escript]/trunk/dudley/test/python/RT2D.py
ViewVC logotype

Annotation of /trunk/dudley/test/python/RT2D.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (21 months ago) by jfenwick
File MIME type: text/x-python
File size: 3450 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 sshaw 5706
2     ##############################################################################
3     #
4 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
5 sshaw 5706 # http://www.uq.edu.au
6     #
7     # Primary Business: Queensland, Australia
8 jfenwick 6112 # Licensed under the Apache License, version 2.0
9     # http://www.apache.org/licenses/LICENSE-2.0
10 sshaw 5706 #
11     # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14     #
15     ##############################################################################
16    
17     from __future__ import print_function, division
18    
19 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 sshaw 5706 http://www.uq.edu.au
21     Primary Business: Queensland, Australia"""
22 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
23     http://www.apache.org/licenses/LICENSE-2.0"""
24 sshaw 5706 __url__="https://launchpad.net/escript-finley"
25    
26 jfenwick 3774 from esys.escript import *
27     from esys.escript.models import StokesProblemCartesian
28     from esys.dudley import Rectangle
29     from esys.weipa import saveVTK
30    
31    
32     #physical properties
33 sshaw 5705 rho1 = 1000 #fluid density on bottom
34     rho2 = 1010 #fluid density on top
35     eta1 = 100.0 #fluid viscosity on bottom
36     eta2 = 100.0 #fluid viscosity on top
37 jfenwick 3774 g=10.0
38    
39     #solver settings
40     dt = 0.001
41     t_step = 0
42     t_step_end = 2000
43     TOL = 1.0e-5
44     max_iter=400
45     verbose=True
46     useUzawa=True
47    
48     #define mesh
49     l0=0.9142
50     l1=1.0
51     n0=200
52     n1=200
53    
54     mesh=Rectangle(l0=l0, l1=l1, order=2, n0=n0, n1=n1)
55     #get mesh dimensions
56     numDim = mesh.getDim()
57     #get element size
58     h = Lsup(mesh.getSize())
59    
60     #level set parameters
61     tolerance = 1.0e-6
62     reinit_max = 30
63     reinit_each = 3
64     alpha = 1
65     smooth = alpha*h
66    
67     #boundary conditions
68     x = mesh.getX()
69     #left + bottom + right + top
70     b_c = whereZero(x[0])*[1.0,0.0] + whereZero(x[1])*[1.0,1.0] + whereZero(x[0]-l0)*[1.0,0.0] + whereZero(x[1]-l1)*[1.0,1.0]
71    
72     velocity = Vector(0.0, ContinuousFunction(mesh))
73     pressure = Scalar(0.0, ContinuousFunction(mesh))
74     Y = Vector(0.0,Function(mesh))
75    
76     #define initial interface between fluids
77     xx = mesh.getX()[0]
78     yy = mesh.getX()[1]
79     func = Scalar(0.0, ContinuousFunction(mesh))
80     h_interface = Scalar(0.0, ContinuousFunction(mesh))
81     h_interface = h_interface + (0.02*cos(math.pi*xx/l0) + 0.2)
82     func = yy - h_interface
83     func_new = func.interpolate(ReducedSolution(mesh))
84    
85     #Stokes Cartesian
86     solution=StokesProblemCartesian(mesh,debug=True)
87     solution.setTolerance(TOL)
88    
89     #level set
90     levelset = LevelSet(mesh, func_new, reinit_max, reinit_each, tolerance, smooth)
91    
92     while t_step <= t_step_end:
93     #update density and viscosity
94     rho = levelset.update_parameter(rho1, rho2)
95     eta = levelset.update_parameter(eta1, eta2)
96    
97     #get velocity and pressure of fluid
98     Y[1] = -rho*g
99     solution.initialize(fixed_u_mask=b_c,eta=eta,f=Y)
100     velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,useUzawa=useUzawa)
101    
102     #update the interface
103     func = levelset.update_phi(velocity, dt, t_step)
104    
105     print("##########################################################")
106     print("time step:", t_step, " completed with dt:", dt)
107     print("Velocity: min =", inf(velocity), "max =", Lsup(velocity))
108     print("##########################################################")
109    
110     #save interface, velocity and pressure
111     saveVTK("phi2D.%2.4i.vtu"%t_step,interface=func,velocity=velocity,pressure=pressure)
112     #courant condition
113     dt = 0.4*Lsup(mesh.getSize())/Lsup(velocity)
114     t_step += 1

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/dudley/test/python/RT2D.py:5567-5588 /branches/lapack2681/finley/test/python/RT2D.py:2682-2741 /branches/pasowrap/dudley/test/python/RT2D.py:3661-3674 /branches/py3_attempt2/dudley/test/python/RT2D.py:3871-3891 /branches/restext/finley/test/python/RT2D.py:2610-2624 /branches/ripleygmg_from_3668/dudley/test/python/RT2D.py:3669-3791 /branches/symbolic_from_3470/dudley/test/python/RT2D.py:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/test/python/RT2D.py:3517-3974 /branches/trilinos_from_5897/dudley/test/python/RT2D.py:5898-6118 /release/4.0/dudley/test/python/RT2D.py:5380-5406 /trunk/ripley/test/python/dudley/test/python/RT2D.py:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26