/[escript]/trunk/doc/examples/geotutorial/steadystate_variablek.py
ViewVC logotype

Diff of /trunk/doc/examples/geotutorial/steadystate_variablek.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5287 by jfenwick, Wed Apr 9 05:41:57 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 24  __url__="https://launchpad.net/escript-f Line 24  __url__="https://launchpad.net/escript-f
24  # import tools  # import tools
25  from esys.escript import *  from esys.escript import *
26  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
27  from esys.dudley import Rectangle  try:
28        from esys.dudley import Rectangle
29        HAVE_DUDLEY = True
30    except ImportError:
31        HAVE_DUDLEY = False
32  from esys.weipa import saveVTK  from esys.weipa import saveVTK
33  # dimensions:  
34  L0=1.;L1=1.;  if not HAVE_DUDLEY:
35  # height of k interface:      print("Dudley module not available")
36  H=L1*0.75  else:
37  # bottom temperature:      # dimensions:
38  T_bot=100      L0=1.;L1=1.;
39  # location, size and value of heat source      # height of k interface:
40  xc=[0.3,0.4]; r=0.1; Qc=3000      H=L1*0.75
41  # two values for k      # bottom temperature:
42  k0=1; k1=10      T_bot=100
43  # create domain:      # location, size and value of heat source
44  mydomain=Rectangle(l0=L0,l1=L1,n0=20,n1=20)      xc=[0.3,0.4]; r=0.1; Qc=3000
45  x=mydomain.getX()      # two values for k
46  # set variable k      k0=1; k1=10
47  k=k0+(k1-k0)*wherePositive(x[1]-H)      # create domain:
48  # boundary temperature      mydomain=Rectangle(l0=L0,l1=L1,n0=20,n1=20)
49  T_D=T_bot/L1*(L1-x[1])      x=mydomain.getX()
50  # heat source      # set variable k
51  Q=Qc*whereNegative(length(x-xc)-r)      k=k0+(k1-k0)*wherePositive(x[1]-H)
52  # create PDE and set coefficients:      # boundary temperature
53  mypde=LinearPDE(mydomain)      T_D=T_bot/L1*(L1-x[1])
54  mypde.setSymmetryOn()      # heat source
55  # set PDE coefficients:      Q=Qc*whereNegative(length(x-xc)-r)
56  mypde.setValue(A=k*kronecker(mydomain),Y=Q, r=T_D, \      # create PDE and set coefficients:
57                  q=whereZero(x[1])+whereZero(x[1]-L1))      mypde=LinearPDE(mydomain)
58  # get temperature:      mypde.setSymmetryOn()
59  T=mypde.getSolution()      # set PDE coefficients:
60  # save as VTK for visualisation:      mypde.setValue(A=k*kronecker(mydomain),Y=Q, r=T_D, \
61  saveVTK("u.vtu",T=T)                      q=whereZero(x[1])+whereZero(x[1]-L1))
62        # get temperature:
63        T=mypde.getSolution()
64        # save as VTK for visualisation:
65        saveVTK("u.vtu",T=T)
66    

Legend:
Removed from v.5287  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26