Skip to content

Dambreak simulation using OpenFOAM, Pythonflu, Mayavi2

April 17, 2011

In this tutorial am going to demonstrate dambreak simulation using Openfoam, pythonflu, mayavi2 in Ubuntu 10.10.

Install following required package

1. OpenFOAM (see the link: https://dhastha.wordpress.com/2011/03/28/openfoam-introduction/)

2. Mayavi2 (link: https://dhastha.wordpress.com/2010/12/18/mayavi2-introduction/)

3. Pythonflu (link: https://dhastha.wordpress.com/2011/04/17/pythonflu-intro/)

copy  /opt/openfoam171/tutorials  into specified folder.

Open your terminal and type

cd tutorials/multiphase/interFoam/laminar/damBreak

then type following commands

blockMesh

setFields

Copy the following python code and run it as “python dambreak.py”

#!/usr/bin/env python

#----------------------------------------------------------------------------
def _createFields( runTime, mesh ):
    from Foam.OpenFOAM import ext_Info, nl
    from Foam.OpenFOAM import IOdictionary, IOobject, word, fileName
    from Foam.finiteVolume import volScalarField
        
    ext_Info() << "Reading field p\n" << nl
    p_rgh = volScalarField( IOobject( word( "p_rgh" ),
                                  fileName( runTime.timeName() ),
                                  mesh,
                                  IOobject.MUST_READ,
                                  IOobject.AUTO_WRITE ),
                        mesh )
    
    ext_Info() << "Reading field alpha1\n" << nl
    alpha1 = volScalarField( IOobject( word( "alpha1" ),
                                  fileName( runTime.timeName() ),
                                  mesh,
                                  IOobject.MUST_READ,
                                  IOobject.AUTO_WRITE ),
                                mesh )
    
    ext_Info() << "Reading field U\n" << nl
    from Foam.finiteVolume import volVectorField
    U = volVectorField( IOobject( word( "U" ),
                                  fileName( runTime.timeName() ),
                                  mesh,
                                  IOobject.MUST_READ,
                                  IOobject.AUTO_WRITE ),
                        mesh )
                        
    from Foam.finiteVolume.cfdTools.incompressible import createPhi
    phi = createPhi( runTime, mesh, U )
    
    
    ext_Info() << "Reading transportProperties\n" << nl
    from Foam.transportModels import twoPhaseMixture
    twoPhaseProperties = twoPhaseMixture(U, phi)
    
    rho1 = twoPhaseProperties.rho1()
    rho2 = twoPhaseProperties.rho2()
    
    # Need to store rho for ddt(rho, U)
    rho = volScalarField( IOobject( word( "rho" ),
                                    fileName( runTime.timeName() ),
                                    mesh,
                                    IOobject.READ_IF_PRESENT ),
                          alpha1 * rho1 + ( 1.0 - alpha1 ) * rho2,
                          alpha1.ext_boundaryField().types() )
    rho.oldTime()
    
    # Mass flux
    # Initialisation does not matter because rhoPhi is reset after the
    # alpha1 solution before it is used in the U equation.
    from Foam.finiteVolume import surfaceScalarField
    rhoPhi = surfaceScalarField( IOobject( word( "rho*phi" ),
                                           fileName( runTime.timeName() ),
                                           mesh,
                                           IOobject.NO_READ,
                                           IOobject.NO_WRITE ),
                                 rho1 * phi )
    
    # Construct interface from alpha1 distribution
    from Foam.transportModels import interfaceProperties
    interface = interfaceProperties( alpha1, U, twoPhaseProperties )


    # Construct incompressible turbulence model
    from Foam import incompressible
    turbulence = incompressible.turbulenceModel.New( U, phi, twoPhaseProperties ) 
    
    from Foam.finiteVolume.cfdTools.general.include import readGravitationalAcceleration
    g = readGravitationalAcceleration( runTime, mesh)

    
    #dimensionedVector g0(g);

    #Read the data file and initialise the interpolation table
    #interpolationTable timeSeriesAcceleration( runTime.path()/runTime.caseConstant()/"acceleration.dat" );
    
    ext_Info() << "Calculating field g.h\n" << nl
    gh = volScalarField( word( "gh" ), g & mesh.C() )
    ghf = surfaceScalarField( word( "ghf" ), g & mesh.Cf() )
    
    p = volScalarField( IOobject( word( "p" ),
                        fileName( runTime.timeName() ),
                        mesh,
                        IOobject.NO_READ,
                        IOobject.AUTO_WRITE ),
                        p_rgh + rho * gh )

    pRefCell = 0
    pRefValue = 0.0
    
    from Foam.finiteVolume import setRefCell, getRefCellValue
    pRefCell, pRefValue = setRefCell( p, mesh.solutionDict().subDict( word( "PISO" ) ), pRefCell, pRefValue )
    
    if p_rgh.needReference():
       p.ext_assign( p + dimensionedScalar( word( "p" ),
                                           p.dimensions(),
                                           pRefValue - getRefCellValue(p, pRefCell) ) )
       p_rgh.ext_assign( p - rho * gh )
       pass

    return p_rgh, p, alpha1, U, phi, rho1, rho2, rho, rhoPhi, twoPhaseProperties, pRefCell, pRefValue, interface, turbulence, g, gh, ghf
    

#--------------------------------------------------------------------------------------
def setDeltaT( runTime, adjustTimeStep, maxCo, CoNum, maxAlphaCo, alphaCoNum, maxDeltaT ):
    from Foam.OpenFOAM import SMALL
    if adjustTimeStep:
       maxDeltaTFact = min( maxCo / ( CoNum + SMALL ), maxAlphaCo / ( alphaCoNum + SMALL ) )
       deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
       
       runTime.setDeltaT( min( deltaTFact * runTime.deltaT().value(), maxDeltaT ) )
       from Foam.OpenFOAM import ext_Info, nl
       ext_Info() << "deltaT = " <<  runTime.deltaT().value() << nl
       pass
    return runTime


#--------------------------------------------------------------------------------------
def alphaCourantNo( runTime, mesh, alpha1, phi ):
    from Foam.OpenFOAM import readScalar, word
    maxAlphaCo = readScalar( runTime.controlDict().lookup(word( "maxAlphaCo" ) ) )
    
    alphaCoNum = 0.0
    meanAlphaCoNum = 0.0
    
    if mesh.nInternalFaces():
       from Foam import fvc
       alphaf = fvc.interpolate(alpha1)

       SfUfbyDelta = (alphaf - 0.01).pos() * ( 0.99 - alphaf ).pos() * mesh.deltaCoeffs() * phi.mag()
       alphaCoNum = ( SfUfbyDelta/mesh.magSf() ).ext_max().value() * runTime.deltaT().value()
       meanAlphaCoNum = ( SfUfbyDelta.sum() / mesh.magSf().sum() ).value() * runTime.deltaT().value()
       pass
    
    from Foam.OpenFOAM import ext_Info, nl
    ext_Info() << "Interface Courant Number mean: " << meanAlphaCoNum << " max: " << alphaCoNum << nl
    
    return maxAlphaCo, alphaCoNum, meanAlphaCoNum


#----------------------------------------------------------------------------------------
def correctPhi( runTime, mesh, phi, p, p_rgh, rho, U, cumulativeContErr, nNonOrthCorr, pRefCell, pRefValue ):
    
    from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
    cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )
    from Foam.OpenFOAM import wordList
    from Foam.finiteVolume import zeroGradientFvPatchScalarField
    pcorrTypes = wordList( p_rgh.ext_boundaryField().size(), zeroGradientFvPatchScalarField.typeName )
    
    from Foam.finiteVolume import fixedValueFvPatchScalarField
    for i in range( p.ext_boundaryField().size() ):
       if p_rgh.ext_boundaryField()[i].fixesValue():
          pcorrTypes[i] = fixedValueFvPatchScalarField.typeName
          pass
       pass
    
    from Foam.OpenFOAM import IOdictionary, IOobject, word, fileName, dimensionedScalar
    from Foam.finiteVolume import volScalarField
    pcorr = volScalarField( IOobject( word( "pcorr" ),
                                      fileName( runTime.timeName() ),
                                      mesh,
                                      IOobject.NO_READ,
                                      IOobject.NO_WRITE ),
                            mesh,
                            dimensionedScalar( word( "pcorr" ), p_rgh.dimensions(), 0.0 ),
                            pcorrTypes  )
    
    from Foam.OpenFOAM import dimTime
    rAUf = dimensionedScalar( word( "(1|A(U))" ), dimTime / rho.dimensions(), 1.0)
    
    from Foam.finiteVolume import adjustPhi
    adjustPhi(phi, U, pcorr)
    
    from Foam import fvc, fvm
    for nonOrth in range( nNonOrthCorr + 1 ):
        pcorrEqn = fvm.laplacian( rAUf, pcorr ) == fvc.div( phi )

        pcorrEqn.setReference(pRefCell, pRefValue)
        pcorrEqn.solve()

        if nonOrth == nNonOrthCorr:
           phi.ext_assign( phi  - pcorrEqn.flux() )
           pass
    
    from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
    cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )
    
    pass
    

#--------------------------------------------------------------------------------------
def alphaEqn( mesh, phi, alpha1, rhoPhi, rho1, rho2, interface, nAlphaCorr ):
    from Foam.OpenFOAM import word 
    alphaScheme = word( "div(phi,alpha)" )
    alpharScheme = word( "div(phirb,alpha)" )
    
    from Foam.finiteVolume import surfaceScalarField
    phic = surfaceScalarField( ( phi / mesh.magSf() ).mag() )
    phic.ext_assign( ( interface.cAlpha() * phic ).ext_min( phic.ext_max() ) )
    phir = phic * interface.nHatf()
    
    from Foam import fvc
    from Foam import MULES

    for aCorr in range( nAlphaCorr ):
       phiAlpha = fvc.flux( phi, alpha1, alphaScheme ) + fvc.flux( -fvc.flux( -phir, 1.0 - alpha1, alpharScheme ), alpha1, alpharScheme )
       MULES.explicitSolve( alpha1, phi, phiAlpha, 1, 0 )
       
       rhoPhi.ext_assign( phiAlpha * ( rho1 - rho2 ) + phi * rho2 )

       pass
    from Foam.OpenFOAM import ext_Info, nl
    ext_Info() << "Liquid phase volume fraction = " << alpha1.weightedAverage( mesh.V() ).value() \
               << "  Min(alpha1) = " << alpha1.ext_min().value() \
               << "  Max(alpha1) = " << alpha1.ext_max().value() < 1):
       totalDeltaT = runTime.deltaT()
       rhoPhiSum = 0.0 * rhoPhi
       from Foam.finiteVolume import subCycle_volScalarField
       alphaSubCycle = subCycle_volScalarField(alpha1, nAlphaSubCycles)
       for item in alphaSubCycle: 
           alphaEqn( mesh, phi, alpha1, rhoPhi, rho1, rho2, interface, nAlphaCorr )
           rhoPhiSum.ext_assign( rhoPhiSum + ( runTime.deltaT() / totalDeltaT ) * rhoPhi )
           pass
       # To make sure that variable in the local scope will be destroyed
       # - during destruction of this variable it performs some important actions
       # - there is a difference between C++ and Python memory management, namely
       # if C++ automatically destroys stack variables when they exit the scope,
       # Python relay its memory management of some "garbage collection" algorithm
       # that do not provide predictable behavior on the "exit of scope"
       del alphaSubCycle
       
       rhoPhi.ext_assign( rhoPhiSum )
    else:
       alphaEqn( mesh, phi, alpha1, rhoPhi, rho1, rho2, interface, nAlphaCorr )
       pass
    interface.correct()
    
    rho == alpha1 * rho1 + ( 1.0 - alpha1 ) * rho2

    pass


#--------------------------------------------------------------------------------------
def _UEqn( mesh, alpha1, U, p, p_rgh, ghf, rho, rhoPhi, turbulence, g, twoPhaseProperties, interface, momentumPredictor ):
    from Foam.OpenFOAM import word
    from Foam.finiteVolume import surfaceScalarField
    from Foam import fvc
    muEff = surfaceScalarField( word( "muEff" ),
                                twoPhaseProperties.muf() + fvc.interpolate( rho * turbulence.ext_nut() ) )
    from Foam import fvm

    UEqn = fvm.ddt( rho, U ) + fvm.div( rhoPhi, U ) - fvm.laplacian( muEff, U ) - ( fvc.grad( U ) & fvc.grad( muEff ) )
    
    UEqn.relax()

    if momentumPredictor:
       from Foam.finiteVolume import solve
       solve( UEqn == \
                   fvc.reconstruct( ( fvc.interpolate( interface.sigmaK() ) * fvc.snGrad( alpha1 ) - ghf * fvc.snGrad( rho ) \
                                                                                                 - fvc.snGrad( p_rgh ) ) * mesh.magSf() ) )
       pass
    
    return UEqn


#--------------------------------------------------------------------------------------
def _pEqn( runTime, mesh, UEqn, U, p, p_rgh, gh, ghf, phi, alpha1, rho, g, interface, corr, nCorr, nNonOrthCorr, pRefCell, pRefValue, cumulativeContErr ):
    rAU = 1.0/UEqn.A()
     
    from Foam import fvc
    rAUf = fvc.interpolate( rAU )
    
    U.ext_assign( rAU * UEqn.H() )
    
    from Foam.finiteVolume import surfaceScalarField
    from Foam.OpenFOAM import word
    phiU = surfaceScalarField( word( "phiU" ),
                               ( fvc.interpolate( U ) & mesh.Sf() ) + fvc.ddtPhiCorr( rAU, rho, U, phi ) )
                               
    from Foam.finiteVolume import adjustPhi
    adjustPhi(phiU, U, p)
    
    phi.ext_assign( phiU + ( fvc.interpolate( interface.sigmaK() ) * fvc.snGrad( alpha1 ) - ghf * fvc.snGrad( rho ) )*rAUf*mesh.magSf() )

    from Foam import fvm
    for nonOrth in range( nNonOrthCorr + 1 ):
        p_rghEqn = fvm.laplacian( rAUf, p_rgh ) == fvc.div( phi ) 
        p_rghEqn.setReference( pRefCell, pRefValue )

        p_rghEqn.solve( mesh.solver( p_rgh.select(corr == nCorr-1 and nonOrth == nNonOrthCorr) ) )
        
        if nonOrth == nNonOrthCorr:
           phi.ext_assign( phi - p_rghEqn.flux() )
           pass
        pass
    
    U.ext_assign( U + rAU * fvc.reconstruct( ( phi - phiU ) / rAUf ) )
    U.correctBoundaryConditions()

    from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
    cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )
    
    p == p_rgh + rho * gh

    if p_rgh.needReference():
       from Foam.OpenFOAM import pRefValue
       p.ext_assign( p + dimensionedScalar( word( "p" ),
                                            p.dimensions(),
                                            pRefValue - getRefCellValue(p, pRefCell) ) )
       p_rgh.ext_assign( p - rho * gh )
       pass
    pass
    
    
#--------------------------------------------------------------------------------------
def main_standalone( argc, argv ):

    from Foam.OpenFOAM.include import setRootCase
    args = setRootCase( argc, argv )

    from Foam.OpenFOAM.include import createTime
    runTime = createTime( args )

    from Foam.OpenFOAM.include import createMesh
    mesh = createMesh( runTime )
    
    from Foam.finiteVolume.cfdTools.general.include import readPISOControls
    piso, nCorr, nNonOrthCorr, momentumPredictor, transonic, nOuterCorr = readPISOControls( mesh )
    
    from Foam.finiteVolume.cfdTools.general.include import initContinuityErrs
    cumulativeContErr = initContinuityErrs()
    
    p_rgh, p, alpha1, U, phi, rho1, rho2, rho, rhoPhi, twoPhaseProperties, pRefCell, \
                                    pRefValue, interface, turbulence, g, gh, ghf = _createFields( runTime, mesh )

    from Foam.finiteVolume.cfdTools.general.include import readTimeControls
    adjustTimeStep, maxCo, maxDeltaT = readTimeControls( runTime )
    
    correctPhi( runTime, mesh, phi, p, p_rgh, rho, U, cumulativeContErr, nNonOrthCorr, pRefCell, pRefValue)
    
    from Foam.finiteVolume.cfdTools.incompressible import CourantNo
    CoNum, meanCoNum = CourantNo( mesh, phi, runTime )
    
    from Foam.finiteVolume.cfdTools.general.include import setInitialDeltaT
    runTime = setInitialDeltaT( runTime, adjustTimeStep, maxCo, maxDeltaT, CoNum )
    
    from Foam.OpenFOAM import ext_Info, nl
    ext_Info() << "\nStarting time loop\n" << nl
    
    while runTime.run() :
                
        piso, nCorr, nNonOrthCorr, momentumPredictor, transonic, nOuterCorr = readPISOControls( mesh )
        adjustTimeStep, maxCo, maxDeltaT = readTimeControls( runTime )
        CoNum, meanCoNum = CourantNo( mesh, phi, runTime )        
        maxAlphaCo, alphaCoNum, meanAlphaCoNum = alphaCourantNo( runTime, mesh, alpha1, phi )
        runTime = setDeltaT(  runTime, adjustTimeStep, maxCo, CoNum, maxAlphaCo, alphaCoNum, maxDeltaT )
        
        runTime.increment()
        ext_Info() << "Time = " << runTime.timeName() << nl << nl
        
        twoPhaseProperties.correct()
     
        alphaEqnSubCycle( runTime, piso, mesh, phi, alpha1, rho, rhoPhi, rho1, rho2, interface )
        
        UEqn = _UEqn( mesh, alpha1, U, p, p_rgh, ghf, rho, rhoPhi, turbulence, g, twoPhaseProperties, interface, momentumPredictor )

        # --- PISO loop
        for corr in range( nCorr ):
            _pEqn( runTime, mesh, UEqn, U, p, p_rgh, gh, ghf, phi, alpha1, rho, g, interface, corr, nCorr, nNonOrthCorr, pRefCell, pRefValue, cumulativeContErr )
            pass
        
        from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
        cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )

        turbulence.correct()

        runTime.write()

        ext_Info() << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << \
              "  ClockTime = " << runTime.elapsedClockTime() << " s" << nl << nl
        
        pass

    ext_Info() << "End\n" <=", "010700" ):
   if __name__ == "__main__" :
      argv = sys.argv
      if len( argv ) > 1 and argv[ 1 ] == "-test":
         argv = None
         test_dir= os.path.join( os.environ[ "PYFOAM_TESTING_DIR" ],'cases', 'propogated', 'r1.7.0',  'multiphase','interFoam', 'laminar', 'damBreak' )
         argv = [ __file__, "-case", test_dir ]
         pass
      os._exit( main_standalone( len( argv ), argv ) )


#--------------------------------------------------------------------------------------

then convert the openfoam data into VTK data by typing following command

foamToVTK

It will create a VTK files. then

cd VTK
mayavi2 damBreak_0.vtk

It will load vtk files into mayavi2. Then add surface module

select Visualize->Modules->Surface

Use Mayavi object editor to change various timestep and see the instance changes

Followings are the images at different timestep

15 Comments leave one →
  1. April 19, 2011 9:18 am

    Hy,

    I’m using interFoam in a relatively simple setup. I have a tank filled with water. The top of the tank is open to atmosphere. Close to the bottom the tank has an outlet into a large resevoir of water. The free surface lavel of the reservoir is at -1m. As far as I know the BC for the pressure at the outlet has to be set to 0, as no hydrostatic pressure has to be considered. BUT!!! How can I make the code understand, that my free surface level outside the domain is at -1m? This must be given for sure but I don’t know where.

    Could anybody give me a hint?

    Thanks in advance, Tony

  2. H.T.Phuc permalink
    June 10, 2011 2:30 am

    Hi Dhastha !
    I copied python code and ran it as “python dambreak.py”. However, appeared an error :
    “…..
    root@phuc:~/OpenFOAM/root-1.7.1/run/tutorials/multiphase/interFoam/laminar/damBreak/VTK# mayavi2 damBreak_0.vtk
    ERROR|2011-06-10 09:18:55,179|invalid syntax (dambreak.py, line 214)
    Traceback (most recent call last):
    File “/usr/lib/python2.6/dist-packages/enthought/mayavi/action/save_load.py”, line 116, in perform
    execfile(dialog.path, g, g)
    File “/home/phucdaigia/OpenFOAMPhuc/root-1.7.1/run/tutorials/multiphase/interFoam/laminar/damBreak/dambreak.py”, line 214
    << " Max(alpha1) = " << alpha1.ext_max().value() < 1):
    ^
    SyntaxError: invalid syntax
    ……"
    Please, to help me !
    Thanks !

  3. davide.conti86@libero.it permalink
    June 10, 2011 10:05 am

    Hello to everyone,
    I’m interested in understanding what MULES is (it should solve the gamma (or alpha) equation) and how it works. can anyone give me a hint? I couldn’t find any info on the web and no references.
    thanks a lot,
    Davide

  4. May 3, 2013 3:27 am

    Most property managers also collect a fee for lease renewal.
    Moreover, the perfect advisor about Florida and Orlando dealings will be “Property Management Companies Florida”.
    Upload Logo (Same logo would be published on contact information, voucher, itinerary & invoice).

  5. May 16, 2013 8:34 pm

    However a bigger surprise was waiting for him whenever for pennies and end up with a profit that’s beyond the belief of other people in additional types of companies. Ensure that a titles are effective plus catchy by utilizing terms which can evoke 1,100-square foot ranch. The employ and overuse of antibiotics is the biggest contributing element in the creation of yeast overgrowth In fact, this system has been responsible for more individual fortunes over the last 80 years than any additional financial practice!. Jobs of all kinds plus sizes need a method instant riches haven’t you?

    Some of dumpster sizes these sound tempting.

  6. June 12, 2013 2:13 am

    Terrific article! This is the kind of info that are meant to be
    shared around the internet. Shame on Google for no longer positioning this
    publish higher! Come on over and talk over with my website .
    Thank you =)

  7. July 12, 2013 10:58 am

    Undeniably believe that which you said. Your favorite justification
    seemed to be on the web the simplest factor to keep in mind of.
    I say to you, I definitely get irked whilst folks consider issues that they plainly do not know about.

    You controlled to hit the nail upon the highest and defined out the whole thing without having side effect
    , people can take a signal. Will likely be back
    to get more. Thanks

  8. July 12, 2013 8:24 pm

    Way cool! Some very valid points! I appreciate you penning this post plus
    the rest of the website is also really good.

  9. July 22, 2013 1:24 pm

    It is the best time to make some plans for the future
    and it’s time to be happy. I’ve read this post and if
    I could I desire to suggest you few interesting things or tips.
    Maybe you could write next articles referring to this article.

    I want to read more things about it!

  10. December 14, 2013 4:15 pm

    Lobeline has the exact same stress-free top qualities which
    are observed in nicotine. It could be gotten from any
    sort of medicine establishment or online. Where can I find garcinia cambogia draw out?

  11. June 13, 2014 12:44 pm

    Hi there, I do think your blog could be having browser compatibility problems.

    When I look at your website in Safari, it looks fine however when opening
    in IE, it has some overlapping issues. I simply wanted to provide you with a quick heads up!
    Apart from that, great website!

  12. June 21, 2014 11:29 am

    Great web site. Plenty of useful info here. I am sending it to
    some friends ans also sharing in delicious. And obviously,
    thanks on your sweat!

  13. January 27, 2016 1:01 pm

    The main problem lies in where the fungus is living:
    under the nail. Our therapists will provide hands-on therapy
    to your plantar fascia and calf muscles. The i – Pad 2S aka i – Pad HD aka
    i – Pad 3 is now known as just the i – Pad.

Trackbacks

  1. mayaviFoam: A postprocessing script for OpenFOAM « Li(G)NUx
  2. Two Liquid Mixing simulation using OpenFoam, PythonFlu, Mayavi2 « Scientific Visual Coder

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: