test/python/run_convection.py

##############################################################################
#
# Copyright (c) 2003-2018 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
#
##############################################################################
from __future__ import print_function, division

__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"

#
#
#    a very simple convection model in model frame:
#

import os
import sys
if sys.version_info >= (3,0):
    from io import StringIO
else:
    #specifically to avoid non-unicode default strings
    #being passed to a python3 style StringIO that expects unicode
    from StringIO import StringIO

import esys.escriptcore.utestselect as unittest
from esys.escriptcore.testing import *
from esys.escript import hasFeature
from esys.escript.modelframe import Link,Simulation
from esys.modellib.input import Sequencer,InterpolateOverBox,GaussianProfile,LinearCombination
from esys.modellib.flow import SteadyIncompressibleFlow
from esys.modellib.temperature import TemperatureAdvection
from esys.modellib.materials import SimpleEarthModel,GravityForce
from esys.modellib.visualization import WriteVTK

try:
    import esys.dudley
    HAVE_DUDLEY = True
except ImportError:
    HAVE_DUDLEY = False

try:
    import esys.finley
    from esys.modellib.geometry import RectangularDomain, ScalarConstrainerOverBox,VectorConstrainerOverBox
    HAVE_FINLEY = True
except ImportError:
    HAVE_FINLEY = False

# TODO: once Amesos2 can deal with block matrices uncomment
HAVE_DIRECT = hasFeature("PASO_DIRECT") #or hasFeature('trilinos')

try:
   WORKDIR=os.environ['MODELLIB_WORKDIR']
except KeyError as e:
   WORKDIR='.'

#if this func is inside the test case, the Link() calls blow up
def run(dom, stream):
    temp_val=InterpolateOverBox()
    temp_val.domain=Link(dom,"domain")
    temp_val.value_left_bottom_front=1.
    temp_val.value_right_bottom_front=1.
    temp_val.value_left_top_front=0.
    temp_val.value_right_top_front=0.
    temp_val.value_left_bottom_back=1.
    temp_val.value_right_bottom_back=1.
    temp_val.value_left_top_back=0.
    temp_val.value_right_top_back=0.

    temp_constraints=ScalarConstrainerOverBox()
    temp_constraints.domain=Link(dom)
    temp_constraints.top=1
    temp_constraints.bottom=1

    vel_constraints=VectorConstrainerOverBox()
    vel_constraints.domain=Link(dom)
    vel_constraints.left=[1,0,0]
    vel_constraints.right=[1,0,0]
    vel_constraints.top=[0,1,0]
    vel_constraints.bottom=[0,1,0]
    vel_constraints.front=[0,0,1]
    vel_constraints.back=[0,0,1]

    mat=SimpleEarthModel()
    mat.density0=1.
    mat.viscocity0=1.
    mat.rayleigh_number=10000.
    mat.alpha=0.001

    temp=TemperatureAdvection(debug=True)
    temp.domain=Link(dom)
    temp.density=Link(mat,"density0")
    temp.heat_capacity=Link(mat,"heat_capacity")
    temp.location_fixed_temperature=Link(temp_constraints,"location_of_constraint")
    temp.fixed_temperature=Link(temp_val,"out")
    temp.safety_factor=0.01
    mat.temperature=Link(temp,"temperature")
    grav=GravityForce()
    grav.domain=Link(dom,"domain")
    grav.direction=[0.,-1.,0.]
    grav.density=Link(mat,"density")
    grav.gravity=Link(mat,"gravity")


    vel=SteadyIncompressibleFlow(debug=True)
    vel.domain=Link(dom)
    vel.internal_force=Link(grav,"gravity_force")
    vel.viscosity=Link(mat,"viscosity")
    vel.location_prescribed_velocity=Link(vel_constraints,"location_of_constraint")
    vel.rel_tol=1.e-6
    temp.velocity=Link(vel,"velocity")

    sq=Sequencer()
    sq.t_end=0.001

    vis=WriteVTK()
    vis.t=Link(sq)
    vis.data0=Link(temp,"temperature")
    vis.data1=Link(vel,"velocity")
    vis.dt=0.0001
    vis.filename=os.path.join(WORKDIR,"temp.vtu")

    per=GaussianProfile()
    per.domain=Link(dom)
    per.x_c=[0.5,0.5,0.5]
    per.A=0.0001
    per.width=0.01
    per.r=0

    lc=LinearCombination()
    lc.f0=1.
    lc.v0=Link(per,"out")
    lc.f1=1.
    lc.v1=Link(temp_val,"out")
    temp.temperature=Link(lc,"out")

    s=Simulation([sq,vel_constraints, temp_constraints,Simulation([vel],debug=True),temp,vis],debug=True)
    s.writeXML(stream)
    s.run()

@unittest.skipUnless(HAVE_DIRECT, "Direct solver not available")
class Test_Convection(unittest.TestCase):
    def setUp(self):
        import sys
        self.old = sys.stdout
        sys.stdout = StringIO()

    def tearDown(self):
        import sys
        sys.stdout = self.old

    @unittest.skipUnless(HAVE_FINLEY, "Finley module not available")
    def test_order2(self):
        dom=RectangularDomain()
        dom.order=2
        run(dom, sys.stdout)

    @unittest.skipUnless(HAVE_DUDLEY and HAVE_FINLEY, "Dudley module not available")
    def test_order1(self):
        dom=RectangularDomain(esys.dudley)
        dom.order=1
        run(dom, sys.stdout)

if __name__ == '__main__':
    run_tests(__name__, exit_on_failure=True)
1	jfenwick	3981	##############################################################################
2	ksteube	1312	#
3	jfenwick	6651	# Copyright (c) 2003-2018 by The University of Queensland
4	jfenwick	3981	# http://www.uq.edu.au
5	ksteube	1312	#
6	ksteube	1809	# Primary Business: Queensland, Australia
7	jfenwick	6112	# Licensed under the Apache License, version 2.0
8			# http://www.apache.org/licenses/LICENSE-2.0
9	ksteube	1312	#
10	jfenwick	3981	# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11	jfenwick	4657	# Development 2012-2013 by School of Earth Sciences
12			# Development from 2014 by Centre for Geoscience Computing (GeoComp)
13	jfenwick	3981	#
14			##############################################################################
15	sshaw	5685	from __future__ import print_function, division
16	jgs	147
17	jfenwick	6651	__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
18	jfenwick	3981	http://www.uq.edu.au
19	ksteube	1809	Primary Business: Queensland, Australia"""
20	jfenwick	6112	__license__="""Licensed under the Apache License, version 2.0
21			http://www.apache.org/licenses/LICENSE-2.0"""
22	jfenwick	2344	__url__="https://launchpad.net/escript-finley"
23	ksteube	1809
24	jgs	147	#
25			#
26			# a very simple convection model in model frame:
27			#
28
29	sshaw	5074	import os
30			import sys
31			if sys.version_info >= (3,0):
32			from io import StringIO
33			else:
34			#specifically to avoid non-unicode default strings
35			#being passed to a python3 style StringIO that expects unicode
36			from StringIO import StringIO
37	sshaw	5288
38	caltinay	6325	import esys.escriptcore.utestselect as unittest
39			from esys.escriptcore.testing import *
40			from esys.escript import hasFeature
41			from esys.escript.modelframe import Link,Simulation
42			from esys.modellib.input import Sequencer,InterpolateOverBox,GaussianProfile,LinearCombination
43			from esys.modellib.flow import SteadyIncompressibleFlow
44			from esys.modellib.temperature import TemperatureAdvection
45			from esys.modellib.materials import SimpleEarthModel,GravityForce
46			from esys.modellib.visualization import WriteVTK
47
48	sshaw	5288	try:
49			import esys.dudley
50			HAVE_DUDLEY = True
51			except ImportError:
52			HAVE_DUDLEY = False
53
54			try:
55			import esys.finley
56			from esys.modellib.geometry import RectangularDomain, ScalarConstrainerOverBox,VectorConstrainerOverBox
57			HAVE_FINLEY = True
58			except ImportError:
59			HAVE_FINLEY = False
60
61	caltinay	6325	# TODO: once Amesos2 can deal with block matrices uncomment
62			HAVE_DIRECT = hasFeature("PASO_DIRECT") #or hasFeature('trilinos')
63	sshaw	5288
64			try:
65	gross	728	WORKDIR=os.environ['MODELLIB_WORKDIR']
66	sshaw	5017	except KeyError as e:
67	gross	728	WORKDIR='.'
68
69	sshaw	5017	#if this func is inside the test case, the Link() calls blow up
70			def run(dom, stream):
71	sshaw	4576	temp_val=InterpolateOverBox()
72			temp_val.domain=Link(dom,"domain")
73			temp_val.value_left_bottom_front=1.
74			temp_val.value_right_bottom_front=1.
75			temp_val.value_left_top_front=0.
76			temp_val.value_right_top_front=0.
77			temp_val.value_left_bottom_back=1.
78			temp_val.value_right_bottom_back=1.
79			temp_val.value_left_top_back=0.
80			temp_val.value_right_top_back=0.
81	sshaw	5288
82	sshaw	4576	temp_constraints=ScalarConstrainerOverBox()
83			temp_constraints.domain=Link(dom)
84			temp_constraints.top=1
85			temp_constraints.bottom=1
86	sshaw	5288
87	sshaw	4576	vel_constraints=VectorConstrainerOverBox()
88			vel_constraints.domain=Link(dom)
89			vel_constraints.left=[1,0,0]
90			vel_constraints.right=[1,0,0]
91			vel_constraints.top=[0,1,0]
92			vel_constraints.bottom=[0,1,0]
93			vel_constraints.front=[0,0,1]
94			vel_constraints.back=[0,0,1]
95	sshaw	5288
96	sshaw	4576	mat=SimpleEarthModel()
97			mat.density0=1.
98			mat.viscocity0=1.
99			mat.rayleigh_number=10000.
100			mat.alpha=0.001
101	sshaw	5288
102	sshaw	4576	temp=TemperatureAdvection(debug=True)
103			temp.domain=Link(dom)
104			temp.density=Link(mat,"density0")
105			temp.heat_capacity=Link(mat,"heat_capacity")
106			temp.location_fixed_temperature=Link(temp_constraints,"location_of_constraint")
107			temp.fixed_temperature=Link(temp_val,"out")
108			temp.safety_factor=0.01
109			mat.temperature=Link(temp,"temperature")
110			grav=GravityForce()
111			grav.domain=Link(dom,"domain")
112			grav.direction=[0.,-1.,0.]
113			grav.density=Link(mat,"density")
114			grav.gravity=Link(mat,"gravity")
115	sshaw	5288
116
117	sshaw	4576	vel=SteadyIncompressibleFlow(debug=True)
118			vel.domain=Link(dom)
119			vel.internal_force=Link(grav,"gravity_force")
120			vel.viscosity=Link(mat,"viscosity")
121			vel.location_prescribed_velocity=Link(vel_constraints,"location_of_constraint")
122			vel.rel_tol=1.e-6
123			temp.velocity=Link(vel,"velocity")
124	sshaw	5288
125	sshaw	4576	sq=Sequencer()
126			sq.t_end=0.001
127	sshaw	5288
128	sshaw	4576	vis=WriteVTK()
129			vis.t=Link(sq)
130			vis.data0=Link(temp,"temperature")
131			vis.data1=Link(vel,"velocity")
132			vis.dt=0.0001
133			vis.filename=os.path.join(WORKDIR,"temp.vtu")
134	sshaw	5288
135	sshaw	4576	per=GaussianProfile()
136			per.domain=Link(dom)
137			per.x_c=[0.5,0.5,0.5]
138			per.A=0.0001
139			per.width=0.01
140			per.r=0
141	sshaw	5288
142	sshaw	4576	lc=LinearCombination()
143			lc.f0=1.
144			lc.v0=Link(per,"out")
145			lc.f1=1.
146			lc.v1=Link(temp_val,"out")
147			temp.temperature=Link(lc,"out")
148	sshaw	5288
149	sshaw	4576	s=Simulation([sq,vel_constraints, temp_constraints,Simulation([vel],debug=True),temp,vis],debug=True)
150	sshaw	5017	s.writeXML(stream)
151	sshaw	4576	s.run()
152	jgs	147
153	caltinay	6321	@unittest.skipUnless(HAVE_DIRECT, "Direct solver not available")
154	sshaw	5017	class Test_Convection(unittest.TestCase):
155			def setUp(self):
156			import sys
157			self.old = sys.stdout
158	sshaw	5074	sys.stdout = StringIO()
159	sshaw	5288
160	sshaw	5017	def tearDown(self):
161			import sys
162			sys.stdout = self.old
163	sshaw	5288
164	caltinay	6321	@unittest.skipUnless(HAVE_FINLEY, "Finley module not available")
165	sshaw	5017	def test_order2(self):
166			dom=RectangularDomain()
167			dom.order=2
168			run(dom, sys.stdout)
169	sshaw	5288
170	caltinay	6321	@unittest.skipUnless(HAVE_DUDLEY and HAVE_FINLEY, "Dudley module not available")
171	sshaw	5017	def test_order1(self):
172			dom=RectangularDomain(esys.dudley)
173			dom.order=1
174			run(dom, sys.stdout)
175
176			if __name__ == '__main__':
177			run_tests(__name__, exit_on_failure=True)
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision