Lahendaja koodi ettevalmistamine

cpII-2010-prax09-gmsh-20091210-00

Kõigil Elmeri võrrandilahendajatel on alljärgnev ühtne liides:

SUBROUTINE PoissonSolver( Model,Solver,dt,TransientSimulation )
USE SolverUtils
TYPE(Model) :: Model
TYPE(Solver_t), POINTER :: Solver
REAL(KIND=dp) :: dt
LOGICAL :: TransientSimulation
...
END SUBROUTINE PoissonSolver

Argument Model sisaldab viita kogu simulatsiooni kirjeldusele. Argument Solver sisaldab käesoleva võrrandilahendaja jaoks spetsiifilisi parameetreid. Argumendid dt ja TransientSimulation sisaldavad vastavalt jooksva ajasammu pikkust ja lippu, mis näitab, kas tegemist on staatilise või dünaamilise simulatsiooniga.

ElmerSolver otsib mudeli kirjelduse sisendfaili “Solver” sektsioonist võtmesõna “Procedure” alt. See peaks sisaldama viidet kompileeritud koodile:

  Procedure = "Poisson" "PoissonSolver"

Selle esimene parameeter võrdusmärgist paremal on kompileeritud koodi sisaldava faili nimi ja teine parameeter on alamfunktsiooni nimi, mida sellest koodist otsida.

“Solver” sektsioonis antakse veel väljamuutuja nimi (antud juhul Poisson) ja DOFs/sõlm (antud juhul 1).

Võrrandilahendaja põhiülesanne on lahendada üks või mitu väljamuutujat ElmerSolveri dünaamilise või staatilise iteratsioonitsükli sees. Käesoleval juhul kasutatakse lõplike elementide meetodit Poisson'i võrrandi diskreetimiseks ja siis lahendatakse võrrand ElmerSolveri utilliidiga SolveSystem .

Lahend saadakse järgmisel viisil:

  • Küsitakse mäluruum ElmerSolverilt ja kompilaatorilt muutujate jaoks. Maatriksstruktuur ja mäluruum lahendi ning RHS vektori jaoks on eraldatud juba enne võrrandilahendaja väljakutsumist.

    Maatriksi tüüp on Matrix_t ja seda võib kasutada järgmiste argumentidega:
    	TYPE(Matrix_t), POINTER :: StiffMatrix
    	StiffMatrix => Solver % Matrix
    	
    Tavaliselt pole vaja maatriksitüübi sisemise mälueralduse skeemi või sisemisi välju vaja teada, ent see annab antud viida edasi ElmerSolveri utiliitidele.

    Analoogselt saab kasutada jõuvektorit:
    	REAL(KIND=dp), POINTER :: ForceVector(:)
    	ForceVector => StiffMatrix % RHS
    	
    ja lahendivektorit:
    	TYPE(Variable_t), POINTER :: SOlution
    	Solution => Solver % Variable
    	
    Struktuur Variable_t sisaldab järgmisi välju:
    • DOFs: ühe sõlme vabadusastmete arv. See väärtus on informatsioonilise tähtsusega ja seda pole vajadust muuta.
    • Perm: täisarvude massiiv, mille väärtused on nullist erinevad sõlmede jaoks, mis kuuluvad antud võrrandi poolt arvutatavasse ruumalasse. Kirje Perm(i) sisaldab globaalmaatriksi reaindeksit (1 DOF jaoks) sõlmpunktile i. Võrrandilahendaja ei tohiks seda massiivi muuta.
    • Väärtused: ruum lahendivektori elementidele. Elemendid on järjestatid samal viisil kui maatriksi read, s.t. sõlme n väärtus Potential on salvestatud kui
      		val = Solution % Values( Solution % Perm(n) )
      		
  • Globaalse süsteemi täitmine nullidega, milleks kutsutakse välja funktsioon
    	CALL InitializeToZero( StiffMatrix, ForceVector )
    	
  • Käia läbi elemendid, mille jaoks võrrandit lahendatakse, saada elemendimaatriksid ja -vektorid ning lisada need globaalsele süsteemile:
    	DO i=1,Solver % NumberOfActiveElements
    	CurrentElement => Solver % Mesh % Elements( Solver % &
    	ActiveElements(i) )
    	...
    	CALL LocalMatrix( ... )
    	CALL UpdateGlobalEquations( ... )
    	END DO
    	CALL FinishAssembly( ... )
    	
    Siin on LocalMatrix kasutaja poolt loodud alamfunktsioon elemendimaatriksite ja -vektorite genereerimiseks. Näitekoodis kasutab LocalMatrix kolme funktsiooni ElmerSolveri utiliitidest. Funktsioon
    	  dim = CoordinateSystemDimension()
    	
    annab jooksva koordinaadisüsteemi dimensiooni (1, 2 või 3) sõltuvalt mudelifaili võtmesõnast “Coordinate System”. Funktsioon GaussPoint annab struktuuri, mis sisaldab integreerimispunkti lokaalkoordinaate ja kaalusid:
    	  TYPE(GaussIntegrationPoints_t) :: IntegStuff
    	IntegStuff = GaussPoints( Element )
    	
    GaussIntegrationPoints_t väljad on
    	INTEGER :: n
    	REAL(KIND=dp) :: u(:), v(:), w(:), s(:)
    	
    Täisarv n on valitud punktide arv. Massiivid u, v ja w on punktide lokaalkoordinaadid ja massiiv s sisaldab punktide kaalusid. GaussPoints-funktsiooni võib välja kutsuda ka lisaargumendiga
    	  IntegStuff = GaussPoints( Element, n )
    	
    kui antud elemendi jaoks integreerimispunktide vaikeväärtus pole sobiv.

    Integreerimistsükli sees kutsutakse välja funktsioon ElementInfo:
    	  TYPE(Element_t), POINTER :: Element
    	TYPE(Nodes_t) :: Nodes
    	REAL(KIND=dp) :: U,V,W,detJ, Basis(n), dBasisdx(n,3), &
    	ddBasisddx(n,3,3)
    	stat = ElementInfo( Element, Nodes, U, V, W, detJ, &
    	Basis, dBasisdx, ddBasisddx, .FALSE. )
    	
    Antud funktsioon annab elemendi jakobiaani determinandi (detJ), baasifunktsioonide väärtused (Basis(n)), baasifunktsioonide globaaltuletiste väärtused (dBasisdx(n,3)) ja baasifunktsioonide teist järku tuletiste väärtused (ddBasisddx(n,3,3)). Teist järku tuletised arvutatakse vaid siis, kui teise loogilise lipu väärtus on tõene. Elemendi sees arvutatakse punktis (U,V,W) kõik struktuuride Elements ja Nodes defineeritud väärtused.
  • Seatakse ääretingimused. Käesoleval juhul kasutatakse vaid Dirichlet' ääretingimusi, mida saab seada funktsiooniga SetDirichletBoundaries.
  • Lahendatakse süsteem funktsiooni SolveSystem väljakutsumisega.

Lahendaja Poisson.f90 täielik lähtetekst:

SUBROUTINE PoissonSolver( Model,Solver,dt,TransientSimulation )
!DEC$ATTRIBUTES DLLEXPORT :: PoissonSolver
!--------------------------------------------------------------
!**************************************************************
!
!  Solve the Poisson equation!
!
!  ARGUMENTS:
!
!  TYPE(Model_t) :: Model,
!     INPUT: All model information (mesh, materials, BCs, etc.)
!
!  TYPE(Solver_t) :: Solver
!     INPUT: Linear & nonlinear equation solver options
!
!  REAL(KIND=dp) :: dt,
!     INPUT: Timestep size for time dependent simulations
!
!  LOGICAL :: TransientSimulation
!     INPUT: Steady state or transient simulation
!
!**************************************************************
USE DefUtils
IMPLICIT NONE
!--------------------------------------------------------------
TYPE(Solver_t) :: Solver
TYPE(Model_t) :: Model
REAL(KIND=dp) :: dt
LOGICAL :: TransientSimulation
!--------------------------------------------------------------
! Local variables
!--------------------------------------------------------------
TYPE(Element_t),POINTER :: Element
LOGICAL :: AllocationsDone = .FALSE., Found
INTEGER :: n, t, istat
REAL(KIND=dp) :: Norm
TYPE(ValueList_t), POINTER :: BodyForce
REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), LOAD(:), FORCE(:)
SAVE STIFF, LOAD, FORCE, AllocationsDone
!--------------------------------------------------------------
!Allocate some permanent storage, this is done first time
!only:
!------------------------------------------------------------
IF ( .NOT. AllocationsDone ) THEN
N = Solver % Mesh % MaxElementNodes ! just big enough for
! elemental arrays
ALLOCATE( FORCE(N), LOAD(N), STIFF(N,N), STAT=istat )
IF ( istat /= 0 ) THEN
CALL Fatal( 'PoissonSolve', 'Memory allocation error.' )
END IF
AllocationsDone = .TRUE.
END IF
!Initialize the system and do the assembly:
!------------------------------------------
CALL DefaultInitialize()
DO t=1,Solver % NumberOfActiveElements
!--------------------------------------------------------------
Element => GetActiveElement(t)
n = GetElementNOFNodes()
!--------------------------------------------------------------
LOAD = 0.0d0
BodyForce => GetBodyForce()
IF ( ASSOCIATED(BodyForce) ) &
Load(1:n) = GetReal( BodyForce, 'Source', Found )
!Get element local matrix and rhs vector:
!----------------------------------------
CALL LocalMatrix(  STIFF, FORCE, LOAD, Element, n )
!Update global matrix and rhs vector from local matrix &
! vector:
!--------------------------------------------------------
CALL DefaultUpdateEquations( STIFF, FORCE )
!--------------------------------------------------------------
END DO
!--------------------------------------------------------------
CALL DefaultFinishAssembly()
!
!  Dirichlet BCs:
!  --------------
CALL DefaultDirichletBCs()
!
!  Solve the system and we are done:
!  ---------------------------------
Norm = DefaultSolve()
!--------------------------------------------------------------
CONTAINS
!--------------------------------------------------------------
SUBROUTINE LocalMatrix(  STIFF, FORCE, LOAD, Element, n )
!--------------------------------------------------------------
REAL(KIND=dp) :: STIFF(:,:), FORCE(:), LOAD(:)
INTEGER :: n
TYPE(Element_t), POINTER :: Element
!--------------------------------------------------------------
REAL(KIND=dp) :: Basis(n),dBasisdx(n,3),ddBasisddx(n,3,3)
REAL(KIND=dp) :: detJ, LoadAtIP
LOGICAL :: Stat
INTEGER :: t
TYPE(GaussIntegrationPoints_t) :: IP
TYPE(Nodes_t) :: Nodes
SAVE Nodes
!--------------------------------------------------------------
CALL GetElementNodes( Nodes )
STIFF = 0.0d0
FORCE = 0.0d0
!Numerical integration:
!----------------------
IP = GaussPoints( Element )
DO t=1,IP % n
!Basis function values & derivatives at the integration
! point:
!-------------------------------------------------------
stat = ElementInfo( Element, Nodes, IP % U(t), &
IP % V(t), IP % W(t),  detJ, Basis, dBasisdx, &
ddBasisddx, .FALSE. )
!The source term at the integration point:
!-----------------------------------------
LoadAtIP = SUM( Basis(1:n) * LOAD(1:n) )
!Finally, the elemental matrix & vector:
!---------------------------------------
STIFF(1:n,1:n) = STIFF(1:n,1:n) + IP % s(t) * DetJ * &
MATMUL( dBasisdx, TRANSPOSE( dBasisdx ) )
FORCE(1:n) = FORCE(1:n) + IP % s(t) * DetJ * LoadAtIP &
* Basis(1:n)
END DO
!--------------------------------------------------------------
END SUBROUTINE LocalMatrix
!--------------------------------------------------------------
END SUBROUTINE PoissonSolver
!--------------------------------------------------------------

Lahendaja kood kompileeritakse käsuga:

Temperature1D $ elmerf90 -o Poisson Poisson.f90

This document was translated from LATEX by HEVEA.

Licensed under the GNU Free Documentation License