在Docker里把OF里的完整矩阵抓出来了。


  • OpenFOAM讲师

    use OpenFOAM in docker

    OS: windows 10

    Install docker for windows

    site: https://www.docker.com/

    download: Docker for Windows Installer.exe

    edition: most recent stable, for me, it is Docker version 17.09.0, community edition

    installation instructions: https://docs.docker.com/docker-for-windows/

    test

    after completion of the installation, run Docker for Windows. It will cost some time to start the docker engine.

    right click “Start” button, choose “Windows Powershell”

    docker --version
    ## Docker version 17.09.0-ce, build afdb6d4
    docker-compose --version
    ## docker-compose version 1.16.1, build 6d1ac219
    docker-machine --version
    ## docker-machine.exe version 0.12.2, build 9371605
    docker run hello-world
    ## 
    ## Hello from Docker!
    ## This message shows that your installation appears to be working correctly.
    ## 
    ## To generate this message, Docker took the following steps:
    ##  1. The Docker client contacted the Docker daemon.
    ##  2. The Docker daemon pulled the "hello-world" image from the Docker Hub.
    ##     (amd64)
    ##  3. The Docker daemon created a new container from that image which runs the
    ##     executable that produces the output you are currently reading.
    ##  4. The Docker daemon streamed that output to the Docker client, which sent it
    ##     to your terminal.
    ## 
    ## To try something more ambitious, you can run an Ubuntu container with:
    ##  $ docker run -it ubuntu bash
    ## 
    ## Share images, automate workflows, and more with a free Docker ID:
    ##  https://cloud.docker.com/
    ## 
    ## For more examples and ideas, visit:
    ##  https://docs.docker.com/engine/userguide/
    docker ps -a
    ## CONTAINER ID        IMAGE               COMMAND             CREATED             STATUS                      PORTS               NAMES
    ## 32a621d46d34        hello-world         "/hello"            24 seconds ago      Exited (0) 23 seconds ago                       clever_agnesi
    docker rm 32a621d46d34 #change it to your container ID
    ## 32a621d46d34 
    docker ps -a
    ## CONTAINER ID        IMAGE               COMMAND             CREATED             STATUS                      PORTS               NAMES
    

    download OpenFOAM image

    right click “Start” button, choose “Windows Powershell”

    docker pull openfoamplus/of_v1706_centos73 # OF+
    # or
    docker pull openfoam/openfoam5-paraview54 # OF5
    # list images
    docker image ls
    

    run OpenFOAM

    Right click on the “Docker” icon in the system tray, choose “settings”

    Click “Shared drives”

    choose any drive you want to share, in my case, I choose C, then click “apply”

    you may need to input credentials

    use openfoam+1706 as example.

    press “WIN+R”, input cmd and enter

    rem test
    docker run --rm -v c:/Users:/data alpine ls /data
    
    rem ONLY `C:/Users` and its subfolders can be mounted
    docker run ^
    -i -t ^
    --name myOFplus_1706 ^
    -v c:/Users/dic17007/OpenFOAM:/OF ^
    openfoamplus/of_v1706_centos73 ^
    bash
    
    rem press `Ctrl+p`, `Ctrl+q` to return without stop the container
    docker attach myOFplus_1706
    
    rem press `Ctrl+c` or input `exit` to return and stop the container
    

    Note: you are root user in docker and the password is ofuser2017. Reference

    setup the environment:

    alias of1706="HOME=/OF source $DOCKER_OPENFOAM_PATH" 
    # add to `~/.bashrc`
    # echo 'alias of1706="HOME=/OF source $DOCKER_OPENFOAM_PATH"'>>~/.bashrc
    

    run testcase

    of1706 #activate openfoam
    mkdir -p $FOAM_RUN
    run
    pwd #/OF/OpenFOAM/-v1706/run
    cp $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity -r .
    cd cavity
    foamJob -screen blockMesh
    foamJob -screen icoFoam
    touch a.foam
    
    ## use Ctrl+P, Ctrl+Q to return to `cmd`
    ## use `docker attach myOFplus_1706`
    
    
    

    return to windows, use paraview to do post-process.

    0_1513906857278_cavity.png

    modify OpenFOAM solver

    replication

    run
    cd ..
    mkdir -p applications/solvers
    cd applications/solvers
    # I put my solver here.
    pwd # /OF/OpenFOAM/-v1706/applications/solvers
    cp $FOAM_SOLVERS/incompressible/icoFoam -r .
    mv icoFoam myIcoFoam
    cd myIcoFoam
    mv icoFoam.C myIcoFoam.C
    sed -i s/icoFoam/myIcoFoam/g myIcoFoam.C
    sed -i s/icoFoam/myIcoFoam/g Make/files
    sed -i s/FOAM_APPBIN/FOAM_USER_APPBIN Make/files
    

    make

    wmake
    

    test

    run
    cd cavity
    which myIcoFoam
    foamJob -screen myIcoFoam
    

    modification

    I am trying to output the matrix in COO format ( Coordinate Format). It will be consists of three parts:

    • AA: non-zero values
    • JR: row index
    • JC: column index

    According to the definition of Foam::lduMatrix::Amul() function, there are 4 part of scalar matrix:

    1. diagonal terms: JC, JR=1 … nCells, AA = matrix.diag();

    2. upper terms: JC > JR

      1. JR=matrix.lduAddr().upperAddr()[0...mFaces-1]+1
      2. JC=matrix.lduAddr().lowerAddr()[0...mFaces-1]+1
      3. AA=matrix.upper()
    3. lower terms: JR > JC

      1. JR=matrix.lduAddr().lowerAddr()[0...mFaces-1]+1
      2. JC=matrix.lduAddr().upperAddr()[0...mFaces-1]+1
      3. AA=matrix.lower()
    4. boundary term, only considering single processor problem here. There are 3 types of patch types:

      1. geometric (constraint) type.
      2. basic
      3. derived

      In the following program, I assume there is not coupled interface such as processor patch or cyclic patch.

    reference: Matrix coupling of different processors

    Here is the code.

    • A c++ library cnpy is used to generate npy or npz file for numpy.

      git clone https://github.com/rogersce/cnpy.git
      cd cnpy
      mkdir build
      cd build
      cmake .. -DENABLE_STATIC=ON
      make 
      make install
      
      
    • dumpFvMatrix.H

    #pragma once
    // added by CatDog
    #include<iostream>
    #include<fstream>
    #include<string>
    #include"cnpy.h"
    
    void dumpFvMatrix(string path, const fvScalarMatrix& EqnPtr)
    {
    
    	const label nCells = EqnPtr.diag().size();
    	const label nFaces = EqnPtr.lower().size();
    
    	const scalar* const __restrict__ diagPtr = EqnPtr.diag().begin();
    
    	const label* const __restrict__ uPtr = EqnPtr.lduAddr().upperAddr().begin();
    	const label* const __restrict__ lPtr = EqnPtr.lduAddr().lowerAddr().begin();
    
    	const scalar* const __restrict__ upperPtr = EqnPtr.upper().begin();
    	const scalar* const __restrict__ lowerPtr = EqnPtr.lower().begin();
    	
    	std::vector<scalar> AA(nCells+2*nFaces);
    	std::vector<label> JR(nCells+2*nFaces);
    	std::vector<label> JC(nCells+2*nFaces);
    	
    	// diag
    	for(label cell=0;cell<nCells;cell++)
    	{
    		AA[cell]=diagPtr[cell];
    		JR[cell]=cell;
    		JC[cell]=cell;
    	}
    	
    	for(label face=0;face<nFaces;face++)
    	{
    		AA[face]=upperPtr[face];
    		JR[face]=lPtr[face];
    		JC[face]=uPtr[face];
    	}
    	
    	for(label face=0;face<nFaces;face++)
    	{
    		AA[face]=lowerPtr[face];
    		JR[face]=uPtr[face];
    		JC[face]=lPtr[face];
    	}
    	
    	cnpy::npz_save(path,"nCells",&nCells,{1},"w");
    	cnpy::npz_save(path,"nFaces",&nFaces,{1},"a");
    	cnpy::npz_save(path,"AA",&AA[0],{nCells+2*nFaces},"a");
    	cnpy::npz_save(path,"JR",&JR[0],{nCells+2*nFaces},"a");
    	cnpy::npz_save(path,"JC",&JC[0],{nCells+2*nFaces},"a");
    	
    	return;
    }
    
    • in myIcoFoam.C
    // ...
    
    #include "fvCFD.H"
    #include "pisoControl.H"
    
    #include "dumpFvMatrix.H"
    
    // ...
    
    // Non-orthogonal pressure corrector loop
                while (piso.correctNonOrthogonal())
                {
                    // Pressure corrector
    
                    fvScalarMatrix pEqn
                    (
                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                    );
    
                    pEqn.setReference(pRefCell, pRefValue);
    
                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
    
                    if (piso.finalNonOrthogonalIter())
                    {
                        phi = phiHbyA - pEqn.flux();
                    }
    				
    				if (runTime.timeIndex()==2)
    				{
    					Info<< "TimeIndex = 2, output matrix pEqn"<<endl;
    					dumpFvMatrix("/OF/OpenFOAM/-v1706/run/cavity/pEqn.npz",pEqn);
    				}
    			}
    
    // ...
    
    
    • options file:
    EXE_INC = \
        -I$(LIB_SRC)/finiteVolume/lnInclude \
        -I$(LIB_SRC)/meshTools/lnInclude \
        -I/usr/local/include
    EXE_LIBS = \
        -lfiniteVolume \
        -lmeshTools \
        -Wl,-rpath -Wl,/usr/local/lib -lcnpy 
    
    
    • python script to read the matrix and find LU-SGS's effectiveness.
    import numpy as np
    import scipy as sp
    from scipy.sparse import coo_matrix,tril,triu,diags
    from scipy.linalg import norm 
    
    data=np.load('pEqn.npz')
    nCells=data["nCells"][0]
    nFaces=data["nFaces"][0]
    AA=data['AA']
    JR=data['JR']
    JC=data['JC']
    
    A=coo_matrix((AA,(JR,JC)),shape=(nCells,nCells))
    
    L=tril(A,-1)
    U=triu(A,1)
    D=diags(AA[0:nCells],0)
    Dinv=diags(1.0/AA[0:nCells],0)
    
    delta= L.dot(Dinv).dot(U)
    
    print "norm of error of LU-SGS: norm(L*Dinv*U)=", norm(delta.toarray())/norm(A.toarray())
        
    # a Laplacian operator
    # proof of diagonally dominance
    print abs(np.abs((L+U).toarray()).sum(1)/np.abs(D.toarray()).sum(1)-1).max() 
    # proof of symmetry
    print abs(L-U.T).sum() 
    
    
    ## reference: Jisheng Kou, Yitian Li, A uniparametric LU-SGS method for systems of nonlinear equations, In Applied Mathematics and Computation, Volume 189, Issue 1, 2007, Pages 235-240, ISSN 0096-3003, 
    f=lambda w:norm(((1-w)*(L+U)-w*w*L.dot(Dinv).dot(U)).toarray())/norm(A.toarray())
    for w in np.linspace(0,1,100):
    	print f(w)
    

    最后结果:

    >>>
    >>> print "norm of error of LU-SGS: norm(L*Dinv*U)=", norm(delta.toarray())/norm(A.toarray())
    relative norm of error of LU-SGS: norm(L*Dinv*U)= 0.141875980931
    >>>
    ... # a Laplacian operator
    ... # proof of diagonally dominance
    ... print abs(np.abs((L+U).toarray()).sum(1)/np.abs(D.toarray()).sum(1)-1).max()
    0.5
    >>> # proof of symmetry
    ... print abs(L-U.T).sum()
    0.0