modellib/py_src/probe.py


##############################################################################
#
# Copyright (c) 2003-2020 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)
# Development from 2019 by School of Earth and Environmental Sciences
#
##############################################################################

from __future__ import print_function, division

__copyright__="""Copyright (c) 2003-2020 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"

from esys.escript.modelframe import Model,ParameterSet
from esys.escript import Data
from esys.escript.util import *


class EvaluateExpression(ParameterSet):
       """
       Return the evaluation of an expression at current time t and 
       locations in the domain

       :warning: this class use python's eval function!!!!! 
                 Please use input.InterpolateOverBox is possible!!!!
       :note: Instance variable expression - expression or list of expressions defining 
                 expression value (in)
       :note: Instance variable out - current value of the expression (callable)
       """

       def __init__(self,**kwargs):
           """
           Set up parameters
           """
           super(EvaluateExpression, self).__init__(**kwargs)
           self.declareParameter(domain=None, \
                                 t=0., \
                                 expression="x[0]")

       def out(self):
         x=self.domain.getX()
         t=self.t
         if isinstance(self.expression,str):
           out=eval(self.expression)
         else:
           out=Data(0,(len(self.expression),),x.getFunctionSpace())
           for i in range(len(self.expression)): out[i]=eval(self.expression[i])
         return out

class Probe(Model):
       """
       Tests values against a expression which may depend on time and spatial 
       coordinates.
       
       It prints out the relative error in each time step and the maximum
       relative error over all time steps at the end.

       :warning: this class uses python's eval function!!!!!
       :ivar value: values to be tested (in)
       :ivar expression: expressions defining expression values to test against. If None only value is reported.  (in)
       :ivar line_tag: tag to be used when printing error.  (in)
       :ivar t: current time  (in)
       :ivar max_error: maximum error  (out)
       :ivar t_max: time of maximum error  (out)

       """

       def __init__(self,**kwargs):
           """
           Set up parameters
           """
           super(Probe, self).__init__(**kwargs)
           self.declareParameter(expression=None, \
                                 value=0., \
                                 t=0., \
                                 line_tag="PROBE")

       def doInitialization(self):
           """
           Initializes values
           """
           self.t_max=None
           self.max_error=0.

       def doStepPostprocessing(self,dt):
            t=self.t
            print("%s : time %e"%(self.line_tag,t))
            v=self.value
            x=v.getFunctionSpace().getX()
            if v.getRank()==0:
              if self.expression==None:
                 print("%s : true (min,max)= (%e,%e)"%(self.line_tag,inf(v),sub(v)))
              else:
                ref=eval(self.expression)
                err=Lsup(v-ref)/Lsup(ref)
                print("%s : true (min,max)= (%e,%e)"%(self.line_tag,inf(ref),sup(ref)))
                print("%s : (min,max,rel error)= (%e,%e,%e)"%(self.line_tag,inf(v),sup(v),err))
            else:
              err=0
              for i in range(v.getRank()):
                 vi=v[i]
                 if self.expression==None:
                    print("%s :   component %d: true (min,max)= (%e,%e)"%(self.line_tag,i,inf(vi)))
                 else:
                    refi=eval(self.expression[i])
                    erri=Lsup(vi-refi)/Lsup(refi)
                    print("%s :   component %d true (min,max)= (%e,%e)"%(self.line_tag,i,inf(refi),sup(refi)))
                    print("%s :   component %d (min,max,rel error)= (%e,%e,%e,%e,%e)"%(self.line_tag,i,inf(vi),sup(vi),erri))
                    err=max(err,erri)
              if not self.expression==None: print("%s :   maximum error %e",err)
            
            if not self.expression==None: 
               if err>self.max_error:
                   self.t_max=t
                   self.max_error=err

       def doFinalization(self):
          """
          Print out the maximum error.
          """
          if not self.t_max==None: print("%s : == maximum error %e at time %e == "%(self.line_tag,self.max_error,self.t_max))

# vim: expandtab shiftwidth=4:
1	ksteube	1809
2	jfenwick	3981	##############################################################################
3	ksteube	1312	#
4	uqaeller	6939	# Copyright (c) 2003-2020 by The University of Queensland
5	jfenwick	3981	# http://www.uq.edu.au
6	ksteube	1312	#
7	ksteube	1809	# 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	ksteube	1312	#
11	jfenwick	3981	# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12	jfenwick	4657	# Development 2012-2013 by School of Earth Sciences
13			# Development from 2014 by Centre for Geoscience Computing (GeoComp)
14	uqaeller	6939	# Development from 2019 by School of Earth and Environmental Sciences
15	jfenwick	3981	#
16			##############################################################################
17	jgs	144
18	sshaw	5706	from __future__ import print_function, division
19
20	uqaeller	6939	__copyright__="""Copyright (c) 2003-2020 by The University of Queensland
21	jfenwick	3981	http://www.uq.edu.au
22	ksteube	1809	Primary Business: Queensland, Australia"""
23	jfenwick	6112	__license__="""Licensed under the Apache License, version 2.0
24			http://www.apache.org/licenses/LICENSE-2.0"""
25	jfenwick	2344	__url__="https://launchpad.net/escript-finley"
26	elspeth	628
27	jgs	149	from esys.escript.modelframe import Model,ParameterSet
28	caltinay	4004	from esys.escript import Data
29	jgs	149	from esys.escript.util import *
30	jgs	144
31
32	jgs	148	class EvaluateExpression(ParameterSet):
33	jgs	149	"""
34			Return the evaluation of an expression at current time t and
35			locations in the domain
36	jgs	144
37	jfenwick	2625	:warning: this class use python's eval function!!!!!
38	jgs	149	Please use input.InterpolateOverBox is possible!!!!
39	jfenwick	6647	:note: Instance variable expression - expression or list of expressions defining
40	jfenwick	2061	expression value (in)
41	jfenwick	6647	:note: Instance variable out - current value of the expression (callable)
42	jgs	149	"""
43	jgs	147
44	gross	918	def __init__(self,**kwargs):
45	jgs	149	"""
46			Set up parameters
47			"""
48	gross	918	super(EvaluateExpression, self).__init__(**kwargs)
49	jgs	147	self.declareParameter(domain=None, \
50			t=0., \
51			expression="x[0]")
52
53			def out(self):
54			x=self.domain.getX()
55			t=self.t
56			if isinstance(self.expression,str):
57			out=eval(self.expression)
58			else:
59			out=Data(0,(len(self.expression),),x.getFunctionSpace())
60			for i in range(len(self.expression)): out[i]=eval(self.expression[i])
61			return out
62
63	jgs	144	class Probe(Model):
64	jgs	149	"""
65			Tests values against a expression which may depend on time and spatial
66			coordinates.
67
68			It prints out the relative error in each time step and the maximum
69			relative error over all time steps at the end.
70	jgs	144
71	jfenwick	2625	:warning: this class uses python's eval function!!!!!
72			:ivar value: values to be tested (in)
73			:ivar expression: expressions defining expression values to test against. If None only value is reported. (in)
74			:ivar line_tag: tag to be used when printing error. (in)
75			:ivar t: current time (in)
76			:ivar max_error: maximum error (out)
77			:ivar t_max: time of maximum error (out)
78	jgs	144
79	jgs	149	"""
80	jgs	144
81	gross	918	def __init__(self,**kwargs):
82	jgs	149	"""
83			Set up parameters
84			"""
85	gross	918	super(Probe, self).__init__(**kwargs)
86	jgs	147	self.declareParameter(expression=None, \
87	jgs	144	value=0., \
88	jgs	147	t=0., \
89	jgs	144	line_tag="PROBE")
90
91	jgs	147	def doInitialization(self):
92	jgs	149	"""
93			Initializes values
94			"""
95	jgs	144	self.t_max=None
96	jgs	147	self.max_error=0.
97	jgs	144
98	jgs	147	def doStepPostprocessing(self,dt):
99			t=self.t
100	jfenwick	3892	print("%s : time %e"%(self.line_tag,t))
101	jgs	144	v=self.value
102			x=v.getFunctionSpace().getX()
103			if v.getRank()==0:
104	jgs	147	if self.expression==None:
105	jfenwick	3892	print("%s : true (min,max)= (%e,%e)"%(self.line_tag,inf(v),sub(v)))
106	jgs	144	else:
107	jgs	147	ref=eval(self.expression)
108	jgs	144	err=Lsup(v-ref)/Lsup(ref)
109	jfenwick	3892	print("%s : true (min,max)= (%e,%e)"%(self.line_tag,inf(ref),sup(ref)))
110			print("%s : (min,max,rel error)= (%e,%e,%e)"%(self.line_tag,inf(v),sup(v),err))
111	jgs	144	else:
112			err=0
113			for i in range(v.getRank()):
114			vi=v[i]
115	jgs	147	if self.expression==None:
116	jfenwick	3892	print("%s : component %d: true (min,max)= (%e,%e)"%(self.line_tag,i,inf(vi)))
117	jgs	144	else:
118	jgs	147	refi=eval(self.expression[i])
119	jgs	144	erri=Lsup(vi-refi)/Lsup(refi)
120	jfenwick	3892	print("%s : component %d true (min,max)= (%e,%e)"%(self.line_tag,i,inf(refi),sup(refi)))
121			print("%s : component %d (min,max,rel error)= (%e,%e,%e,%e,%e)"%(self.line_tag,i,inf(vi),sup(vi),erri))
122	jgs	144	err=max(err,erri)
123	jfenwick	3892	if not self.expression==None: print("%s : maximum error %e",err)
124	jgs	144
125	jgs	147	if not self.expression==None:
126			if err>self.max_error:
127			self.t_max=t
128			self.max_error=err
129	jgs	144
130	jgs	147	def doFinalization(self):
131	jgs	149	"""
132	sshaw	4576	Print out the maximum error.
133			"""
134	jfenwick	3892	if not self.t_max==None: print("%s : == maximum error %e at time %e == "%(self.line_tag,self.max_error,self.t_max))
135	jgs	149
136			# vim: expandtab shiftwidth=4:
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision