Lahendaja koodi ettevalmistamine
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. Funktsioondim = 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 onINTEGER :: 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 lisaargumendigaIntegStuff = 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