axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / PDEPACK.spad

)abbrev package PDEPACK AnnaPartialDifferentialEquationPackage
++ Author: Brian Dupee
++ Date Created: June 1996
++ Date Last Updated: December 1997
++ Description: 
++ AnnaPartialDifferentialEquationPackage is an uncompleted
++ package for the interface to NAG PDE routines.  It has been realised that
++ a new approach to solving PDEs will need to be created.

AnnaPartialDifferentialEquationPackage() : SIG == CODE where

 LEDF  ==> List Expression DoubleFloat
 EDF  ==> Expression DoubleFloat
 LDF  ==> List DoubleFloat
 MDF  ==> Matrix DoubleFloat
 DF  ==> DoubleFloat
 LEF  ==> List Expression Float
 EF  ==> Expression Float
 MEF  ==> Matrix Expression Float
 LF  ==> List Float
 F  ==> Float
 LS  ==> List Symbol
 ST  ==> String
 LST  ==> List String
 INT  ==> Integer
 NNI  ==> NonNegativeInteger
 RT  ==> RoutinesTable
 PDEC  ==> Record(start:DF, finish:DF, grid:NNI, boundaryType:INT, 
                   dStart:MDF, dFinish:MDF)
 PDEB  ==> Record(pde:LEDF, constraints:List PDEC,
                   f:List LEDF, st:ST, tol:DF)
 IFL  ==> List(Record(ifail:INT,instruction:ST))
 Entry  ==> Record(chapter:ST, type:ST, domainName: ST, 
                    defaultMin:F, measure:F, failList:IFL, explList:LST)
 Measure  ==> Record(measure:F,name:ST, explanations:LST)

 SIG ==> with

  solve : (NumericalPDEProblem) -> Result
    ++ solve(PDEProblem) is a top level ANNA function to solve numerically 
    ++ a system of partial differential equations.  
    ++
    ++ The method used to perform the numerical
    ++ process will be one of the routines contained in the NAG numerical
    ++ Library.  The function predicts the likely most effective routine
    ++ by checking various attributes of the system of PDE's and calculating
    ++ a measure of compatibility of each routine to these attributes.
    ++
    ++ It then calls the resulting `best' routine.
    ++
    ++ ** At the moment, only Second Order Elliptic Partial Differential
    ++ Equations are solved **

  solve : (NumericalPDEProblem,RT) -> Result
    ++ solve(PDEProblem,routines) is a top level ANNA function to solve numerically a system
    ++ of partial differential equations.  
    ++
    ++ The method used to perform the numerical
    ++ process will be one of the routines contained in the NAG numerical
    ++ Library.  The function predicts the likely most effective routine
    ++ by checking various attributes of the system of PDE's and calculating
    ++ a measure of compatibility of each routine to these attributes.
    ++
    ++ It then calls the resulting `best' routine.
    ++
    ++ ** At the moment, only Second Order Elliptic Partial Differential
    ++ Equations are solved **

  solve : (F,F,F,F,NNI,NNI,LEF,List LEF,ST,DF) -> Result
    ++ solve(xmin,ymin,xmax,ymax,ngx,ngy,pde,bounds,st,tol) is a top level 
    ++ ANNA function to solve numerically a system of partial differential 
    ++ equations.  This is defined as a list of coefficients (\axiom{pde}),
    ++ a grid (\axiom{xmin}, \axiom{ymin}, \axiom{xmax}, \axiom{ymax}, 
    ++ \axiom{ngx}, \axiom{ngy}), the boundary values (\axiom{bounds}) and a
    ++ tolerance requirement (\axiom{tol}).  There is also a parameter 
    ++ (\axiom{st}) which should contain the value "elliptic" if the PDE is
    ++ known to be elliptic, or "unknown" if it is uncertain.  This causes the
    ++ routine to check whether the PDE is elliptic.
    ++
    ++ The method used to perform the numerical
    ++ process will be one of the routines contained in the NAG numerical
    ++ Library.  The function predicts the likely most effective routine
    ++ by checking various attributes of the system of PDE's and calculating
    ++ a measure of compatibility of each routine to these attributes.
    ++
    ++ It then calls the resulting `best' routine.
    ++
    ++ ** At the moment, only Second Order Elliptic Partial Differential
    ++ Equations are solved **

  solve : (F,F,F,F,NNI,NNI,LEF,List LEF,ST) -> Result
    ++ solve(xmin,ymin,xmax,ymax,ngx,ngy,pde,bounds,st) is a top level 
    ++ ANNA function to solve numerically a system of partial differential 
    ++ equations.  This is defined as a list of coefficients (\axiom{pde}),
    ++ a grid (\axiom{xmin}, \axiom{ymin}, \axiom{xmax}, \axiom{ymax}, 
    ++ \axiom{ngx}, \axiom{ngy}) and the boundary values (\axiom{bounds}).  
    ++ A default value for tolerance is used.  There is also a parameter 
    ++ (\axiom{st}) which should contain the value "elliptic" if the PDE is
    ++ known to be elliptic, or "unknown" if it is uncertain.  This causes the
    ++ routine to check whether the PDE is elliptic.
    ++
    ++ The method used to perform the numerical
    ++ process will be one of the routines contained in the NAG numerical
    ++ Library.  The function predicts the likely most effective routine
    ++ by checking various attributes of the system of PDE's and calculating
    ++ a measure of compatibility of each routine to these attributes.
    ++
    ++ It then calls the resulting `best' routine.
    ++
    ++ ** At the moment, only Second Order Elliptic Partial Differential
    ++ Equations are solved **

  measure : (NumericalPDEProblem) -> Measure
    ++ measure(prob) is a top level ANNA function for identifying the most
    ++ appropriate numerical routine from those in the routines table
    ++ provided for solving the numerical PDE
    ++ problem defined by \axiom{prob}.
    ++
    ++ It calls each \axiom{domain} of \axiom{category}
    ++ \axiomType{PartialDifferentialEquationsSolverCategory} in turn to 
    ++ calculate all measures and returns the best the name of 
    ++ the most appropriate domain and any other relevant information.
    ++ It predicts the likely most effective NAG numerical
    ++ Library routine to solve the input set of PDEs
    ++ by checking various attributes of the system of PDEs and calculating
    ++ a measure of compatibility of each routine to these attributes.

  measure : (NumericalPDEProblem,RT) -> Measure
    ++ measure(prob,R) is a top level ANNA function for identifying the most
    ++ appropriate numerical routine from those in the routines table
    ++ provided for solving the numerical PDE
    ++ problem defined by \axiom{prob}.
    ++
    ++ It calls each \axiom{domain} listed in \axiom{R} of \axiom{category}
    ++ \axiomType{PartialDifferentialEquationsSolverCategory} in turn to 
    ++ calculate all measures and returns the best the name of 
    ++ the most appropriate domain and any other relevant information.
    ++ It predicts the likely most effective NAG numerical
    ++ Library routine to solve the input set of PDEs
    ++ by checking various attributes of the system of PDEs and calculating
    ++ a measure of compatibility of each routine to these attributes.

 CODE ==> add

  import PDEB, d03AgentsPackage, ExpertSystemToolsPackage, NumericalPDEProblem

  zeroMeasure:Measure -> Result
  measureSpecific:(ST,RT,PDEB) -> Record(measure:F,explanations:ST)
  solveSpecific:(PDEB,ST) -> Result
  changeName:(Result,ST) -> Result
  recoverAfterFail:(PDEB,RT,Measure,Integer,Result) -> _
    Record(a:Result,b:Measure)

  zeroMeasure(m:Measure):Result ==
    a := coerce(0$F)$AnyFunctions1(F)
    text:= coerce("No available routine appears appropriate")$AnyFunctions1(ST)
    r := construct([[result@Symbol,a],[method@Symbol,text]])$Result
    concat(measure2Result m,r)$ExpertSystemToolsPackage

  measureSpecific(name:ST,R:RT,p:PDEB):Record(measure:F,explanations:ST) ==
    name = "d03eefAnnaType" => measure(R,p)$d03eefAnnaType
    --name = "d03fafAnnaType" => measure(R,p)$d03fafAnnaType
    error("measureSpecific","invalid type name: " name)$ErrorFunctions

  measure(P:NumericalPDEProblem,R:RT):Measure ==
    p:PDEB := retract(P)$NumericalPDEProblem
    sofar := 0$F
    best := "none" :: ST
    routs := copy R
    routs := selectPDERoutines(routs)$RT
    empty?(routs)$RT => 
      error("measure", "no routines found")$ErrorFunctions
    rout := inspect(routs)$RT
    e := retract(rout.entry)$AnyFunctions1(Entry)
    meth := empty()$LST
    for i in 1..# routs repeat
      rout := extract!(routs)$RT
      e := retract(rout.entry)$AnyFunctions1(Entry)
      n := e.domainName
      if e.defaultMin > sofar then
        m := measureSpecific(n,R,p)
        if m.measure > sofar then
          sofar := m.measure
          best := n
        str:LST := [string(rout.key)$Symbol "measure: " 
                    outputMeasure(m.measure)$ExpertSystemToolsPackage " - " 
                     m.explanations]
      else 
        str := [string(rout.key)$Symbol " is no better than other routines"]
      meth := append(meth,str)$LST
    [sofar,best,meth]

  measure(P:NumericalPDEProblem):Measure == measure(P,routines()$RT)

  solveSpecific(p:PDEB,n:ST):Result ==
    n = "d03eefAnnaType" => PDESolve(p)$d03eefAnnaType
    --n = "d03fafAnnaType" => PDESolve(p)$d03fafAnnaType
    error("solveSpecific","invalid type name: " n)$ErrorFunctions

  changeName(ans:Result,name:ST):Result ==
    sy:Symbol := coerce(name "Answer")$Symbol
    anyAns:Any := coerce(ans)$AnyFunctions1(Result)
    construct([[sy,anyAns]])$Result

  recoverAfterFail(p:PDEB,routs:RT,m:Measure,iint:Integer,r:Result):
                                            Record(a:Result,b:Measure) ==
    while positive?(iint) repeat
      routineName := m.name
      s := recoverAfterFail(routs,routineName(1..6),iint)$RT
      s case "failed" => iint := 0
      (s = "no action")@Boolean => iint := 0
      fl := coerce(s)$AnyFunctions1(ST)
      flrec:Record(key:Symbol,entry:Any):=[failure@Symbol,fl]
      m2 := measure(p::NumericalPDEProblem,routs)
      zero?(m2.measure) => iint := 0
      r2:Result := solveSpecific(p,m2.name)
      m := m2
      insert!(flrec,r2)$Result
      r := concat(r2,changeName(r,routineName))$ExpertSystemToolsPackage
      iany := search(ifail@Symbol,r2)$Result
      iany case "failed" => iint := 0
      iint := retract(iany)$AnyFunctions1(Integer)
    [r,m]

  solve(P:NumericalPDEProblem,t:RT):Result ==
    routs := copy(t)$RT
    m := measure(P,routs)
    p:PDEB := retract(P)$NumericalPDEProblem
    zero?(m.measure) => zeroMeasure m
    r := solveSpecific(p,n := m.name)
    iany := search(ifail@Symbol,r)$Result
    iint := 0$Integer
    if (iany case Any) then
      iint := retract(iany)$AnyFunctions1(Integer)
    if positive?(iint) then
      tu:Record(a:Result,b:Measure) := recoverAfterFail(p,routs,m,iint,r)
      r := tu.a
      m := tu.b
    expl := getExplanations(routs,n(1..6))$RoutinesTable
    expla := coerce(expl)$AnyFunctions1(LST)
    explaa:Record(key:Symbol,entry:Any) := ["explanations"::Symbol,expla]
    r := concat(construct([explaa]),r)
    concat(measure2Result m,r)$ExpertSystemToolsPackage

  solve(P:NumericalPDEProblem):Result == solve(P,routines()$RT)

  solve(xmi:F,xma:F,ymi:F,yma:F,nx:NNI,ny:NNI,pe:LEF,bo:List
                LEF,s:ST,to:DF):Result ==
    cx:PDEC := [f2df xmi, f2df xma, nx, 1, empty()$MDF, empty()$MDF]
    cy:PDEC := [f2df ymi, f2df yma, ny, 1, empty()$MDF, empty()$MDF]
    p:PDEB := [[ef2edf e for e in pe],[cx,cy],
                [[ef2edf u for u in w] for w in bo],s,to]
    solve(p::NumericalPDEProblem,routines()$RT)

  solve(xmi:F,xma:F,ymi:F,yma:F,nx:NNI,ny:NNI,pe:LEF,bo:List
                LEF,s:ST):Result ==
    solve(xmi,xma,ymi,yma,nx,ny,pe,bo,s,0.0001::DF)