main.cpp 29.2 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "inmost.h"
using namespace INMOST;

#ifndef M_PI
#define M_PI 3.141592653589
#endif

#if defined(USE_MPI)
#define BARRIER MPI_Barrier(MPI_COMM_WORLD);
#else
#define BARRIER
#endif

//shortcuts
Kirill Terekhov's avatar
Kirill Terekhov committed
20
typedef Storage::bulk bulk;
Kirill Terekhov's avatar
Kirill Terekhov committed
21 22 23 24 25 26
typedef Storage::real real;
typedef Storage::integer integer;
typedef Storage::enumerator enumerator;
typedef Storage::real_array real_array;
typedef Storage::var_array var_array;

ahmadabushaikha's avatar
ahmadabushaikha committed
27
bool rt0 = true;
Kirill Terekhov's avatar
Kirill Terekhov committed
28

Kirill Terekhov's avatar
Kirill Terekhov committed
29
//#define OPTIMIZATION
Kirill Terekhov's avatar
Kirill Terekhov committed
30 31 32

int main(int argc,char ** argv)
{
Kirill Terekhov's avatar
Kirill Terekhov committed
33
    Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
Kirill Terekhov's avatar
Kirill Terekhov committed
34
#if defined(USE_PARTITIONER)
Kirill Terekhov's avatar
Kirill Terekhov committed
35
    Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
36
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
37
    if( argc > 1 )
Kirill Terekhov's avatar
Kirill Terekhov committed
38
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
39
        /*
Kirill Terekhov's avatar
Kirill Terekhov committed
40 41 42 43 44
        std::cout << "Row: " << sizeof(Sparse::Row) << std::endl;
        std::cout << "HessianRow: " << sizeof(Sparse::HessianRow) << std::endl;
        std::cout << "Matrix: " << sizeof(Sparse::Matrix) << std::endl;
        std::cout << "HessianMatrix: " << sizeof(Sparse::HessianMatrix) << std::endl;
        std::cout << "variable: " << sizeof(multivar_expression) << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
45 46 47 48 49 50 51 52 53 54 55 56
        Sparse::Row J;
        Sparse::HessianRow H;
        unknown a = unknown(0.5,1);
        unknown b = unknown(0.1,2);
        (a*a*b*(a+b)).GetHessian(1.0,J,1.0,H);
        std::cout << "J: "; J.Print();
        std::cout << "H: "; H.Print();
        J.Clear();
        (a*a*b*(a+b)).GetJacobian(1.0,J);
        std::cout << "J: "; J.Print();
        return 0;
         */
ahmadabushaikha's avatar
ahmadabushaikha committed
57 58 59 60 61 62 63 64
		if( argc > 2 ) 
		{
			if( std::string(argv[2]) == "MFD" )
				rt0 = false;
			else if ( std::string(argv[2]) == "MHFE" )
				rt0 = true;
		}

Kirill Terekhov's avatar
Kirill Terekhov committed
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
        double ttt; // Variable used to measure timing
        bool repartition = false; // Is it required to redistribute the mesh?
        Mesh * m = new Mesh(); // Create an empty mesh
        { // Load the mesh
            ttt = Timer();
            m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
            if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl;
            if( m->isParallelFileFormat(argv[1]) ) //The format is
            {
                m->Load(argv[1]); // Load mesh from the parallel file format
                repartition = true; // Ask to repartition the mesh
            }
            else if( m->GetProcessorRank() == 0 ) m->Load(argv[1]); // Load mesh from the serial file format
            BARRIER
            if( m->GetProcessorRank() == 0 ) std::cout << "Load the mesh: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
80
        }
ahmadabushaikha's avatar
ahmadabushaikha committed
81 82 83 84 85

		if( rt0 )
			std::cout << "Running MHFE RT0" << std::endl;
		else
			std::cout << "Running MFD" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
86 87 88 89
        
        
#if defined(USE_PARTITIONER)
        if (m->GetProcessorsNumber() > 1 )//&& !repartition) // Currently only non-distributed meshes are supported by Inner_RCM partitioner
Kirill Terekhov's avatar
Kirill Terekhov committed
90
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
            { // Compute mesh partitioning
                ttt = Timer();
                Partitioner p(m); //Create Partitioning object
                p.SetMethod(Partitioner::Inner_RCM,repartition ? Partitioner::Repartition : Partitioner::Partition); // Specify the partitioner
                p.Evaluate(); // Compute the partitioner and store new processor ID in the mesh
                BARRIER
                if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl;
            }
            
            { //Distribute the mesh
                ttt = Timer();
                m->Redistribute(); // Redistribute the mesh data
                m->ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring
                BARRIER
                if( m->GetProcessorRank() == 0 ) std::cout << "Redistribute: " << Timer()-ttt << std::endl;
            }
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
108
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
109 110 111 112 113 114 115 116 117 118 119 120
        
        { // prepare geometrical data on the mesh
            ttt = Timer();
            Mesh::GeomParam table;
            table[CENTROID]    = CELL | FACE; //Compute averaged center of mass
            table[NORMAL]      = FACE;        //Compute normals
            table[ORIENTATION] = FACE;        //Check and fix normal orientation
            table[MEASURE]     = CELL | FACE; //Compute volumes and areas
            //table[BARYCENTER]  = CELL | FACE; //Compute volumetric center of mass
            m->PrepareGeometricData(table); //Ask to precompute the data
            BARRIER
            if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
121
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
122 123 124 125 126 127 128
        
        // data tags for
        Tag tag_P;  // Pressure
        Tag tag_K;  // Diffusion tensor
        Tag tag_F;  // Forcing term
        Tag tag_BC; // Boundary conditions
        Tag tag_W;  // Gradient matrix acting on harmonic points on faces and returning gradient on faces
Kirill Terekhov's avatar
Kirill Terekhov committed
129
        Tag tag_DMP; // Indicates weather local W matrix satisfy DMP condition
Kirill Terekhov's avatar
Kirill Terekhov committed
130 131 132 133 134 135 136
        
        if( m->GetProcessorsNumber() > 1 ) //skip for one processor job
        { // Exchange ghost cells
            ttt = Timer();
            m->ExchangeGhost(1,FACE); // Produce layer of ghost cells
            BARRIER
            if( m->GetProcessorRank() == 0 ) std::cout << "Exchange ghost: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
137
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188
        
        { //initialize data
            if( m->HaveTag("PERM") ) // is diffusion tensor already defined on the mesh? (PERM from permeability)
                tag_K = m->GetTag("PERM"); // get the diffusion tensor
            
            if( !tag_K.isValid() || !tag_K.isDefined(CELL) ) // diffusion tensor was not initialized or was not defined on cells.
            {
                tag_K = m->CreateTag("PERM",DATA_REAL,CELL,NONE,6); // create a new tag for symmetric diffusion tensor K
                for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) // loop over mesh cells
                {
                    Cell cell = m->CellByLocalID(q);
                    real_array K = cell->RealArray(tag_K);
                    // assign a symmetric positive definite tensor K
                    K[0] = 1.0; //XX
                    K[1] = 0.0; //XY
                    K[2] = 0.0; //XZ
                    K[3] = 1.0; //YY
                    K[4] = 0.0; //YZ
                    K[5] = 1.0; //ZZ
                }
                
                m->ExchangeData(tag_K,CELL,0); //Exchange diffusion tensor
            }
            
            if( m->HaveTag("PRESSURE") ) //Is there a pressure on the mesh?
                tag_P = m->GetTag("PRESSURE"); //Get the pressure
            
            if( !tag_P.isValid() || !tag_P.isDefined(CELL) ) // Pressure was not initialized or was not defined on nodes
            {
                srand(1); // Randomization
                tag_P = m->CreateTag("PRESSURE",DATA_REAL,CELL|FACE,NONE,1); // Create a new tag for the pressure
                for(Mesh::iteratorElement e = m->BeginElement(CELL|FACE); e != m->EndElement(); ++e) //Loop over mesh cells
                    e->Real(tag_P) = 0;//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
            }
            
            if( !tag_P.isDefined(FACE) )
            {
                tag_P = m->CreateTag("PRESSURE",DATA_REAL,FACE,NONE,1);
                for(Mesh::iteratorElement e = m->BeginElement(FACE); e != m->EndElement(); ++e) //Loop over mesh cells
                    e->Real(tag_P) = 0;//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
            }
            
            
            
            if( m->HaveTag("BOUNDARY_CONDITION") ) //Is there boundary condition on the mesh?
            {
                tag_BC = m->GetTag("BOUNDARY_CONDITION");
                //initialize unknowns at boundary
            }
            m->ExchangeData(tag_P,CELL|FACE,0); //Synchronize initial solution with boundary unknowns
            tag_W = m->CreateTag("nKGRAD",DATA_REAL,CELL,NONE);
Kirill Terekhov's avatar
Kirill Terekhov committed
189
            tag_DMP = m->CreateTag("isDMP",DATA_BULK,CELL,NONE);
Kirill Terekhov's avatar
Kirill Terekhov committed
190 191
            ttt = Timer();
            //Assemble gradient matrix W on cells
Kirill Terekhov's avatar
Kirill Terekhov committed
192
            int total = 0, dmp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
193
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
194
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
Kirill Terekhov's avatar
Kirill Terekhov committed
195
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
196
            for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
Kirill Terekhov's avatar
Kirill Terekhov committed
197
            {
ahmadabushaikha's avatar
ahmadabushaikha committed
198
				Mesh * mesh = m;
Kirill Terekhov's avatar
Kirill Terekhov committed
199
                Cell cell = m->CellByLocalID(q);
ahmadabushaikha's avatar
ahmadabushaikha committed
200 201 202 203
				ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
				int NF = (int)faces.size(); //number of faces;	
				rMatrix W(NF,NF);
				
ahmadabushaikha's avatar
ahmadabushaikha committed
204
				if( cell->GetGeometricType() == Element::Tet && rt0 ) // RT0 consturction of W matrix 
ahmadabushaikha's avatar
ahmadabushaikha committed
205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379
				{
					double V = cell.Volume();
					double dN[12];
					double J[9];
					double J_T[9];
					double K_inv[9];
					double K_inv_ref1[9];
					double W_RT0[12];
					double B_RT0[16];
					double gauss_pt_xyz[12];
					double gauss_wei[4];
					double tetra[12];
					double xyz_i[3];
					double xyz_j[3];
					int i, j,k,Z;
					double w,l,m, J_det,x,y,z;
					
						
					rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),cell->RealArrayDF(tag_K).size()); //get permeability for the cell
					rMatrix K_inverse (3,3); // inverse of permeability 
					K_inverse= K.Invert(true).first;
					// gauss points for intgration
					gauss_pt_xyz[0]=(5.0+3.0*(sqrt(5.0)))/20.0;
					gauss_pt_xyz[1]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[2]=(5.0-(sqrt(5.0)))/20.0;

					gauss_pt_xyz[3]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[4]=(5.0+3.0*(sqrt(5.0)))/20.0;
					gauss_pt_xyz[5]=(5.0-(sqrt(5.0)))/20.0;

					gauss_pt_xyz[6]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[7]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[8]=(5.0+3.0*(sqrt(5.0)))/20.0;

					gauss_pt_xyz[9]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[10]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[11]=(5.0-(sqrt(5.0)))/20.0;

					// gauss weights 
					gauss_wei[0]=  1.0/24.0;
					gauss_wei[1]=  1.0/24.0;
					gauss_wei[2]=  1.0/24.0;
					gauss_wei[3]= 1.0/24.0;

					// derivative of the finite element shape functions for tetraherdon (local) 
					dN[0]=-1; 
					dN[1]=-1;
					dN[2]=-1;

					dN[3]=1; 
					dN[4]=0;
					dN[5]=0;

					dN[6]=0; 
					dN[7]=1;
					dN[8]=0;

					dN[9]=0; 
					dN[10]=0;
					dN[11]=1;


					// local to reference tranformation 
					ElementArray<Node> node(mesh);
					for(j = 0; j < 4; ++j)
					{
						ElementArray<Node> cell_nodes = cell->getNodes();
						cell_nodes.Subtract(faces[j].getNodes());
						node.Unite(cell_nodes);
					}
					
					for ( j = 0; j < 4; j++ )
					{  

						for ( i = 0; i < 3; i++ )
							tetra[i+j*3]= node[j].Coords()[i]; // <-- X
					}

					J[0]=dN[0]*tetra[0]+dN[3]*tetra[3]+dN[6]*tetra[6]+dN[9]*tetra[9];
					J[1]=dN[0]*tetra[1]+dN[3]*tetra[4]+dN[6]*tetra[7]+dN[9]*tetra[10];
					J[2]=dN[0]*tetra[2]+dN[3]*tetra[5]+dN[6]*tetra[8]+dN[9]*tetra[11];

					J[3]=dN[1]*tetra[0]+dN[4]*tetra[3]+dN[7]*tetra[6]+dN[10]*tetra[9];
					J[4]=dN[1]*tetra[1]+dN[4]*tetra[4]+dN[7]*tetra[7]+dN[10]*tetra[10];
					J[5]=dN[1]*tetra[2]+dN[4]*tetra[5]+dN[7]*tetra[8]+dN[10]*tetra[11];

					J[6]=dN[2]*tetra[0]+dN[5]*tetra[3]+dN[8]*tetra[6]+dN[11]*tetra[9];
					J[7]=dN[2]*tetra[1]+dN[5]*tetra[4]+dN[8]*tetra[7]+dN[11]*tetra[10];
					J[8]=dN[2]*tetra[2]+dN[5]*tetra[5]+dN[8]*tetra[8]+dN[11]*tetra[11];

					w = J[0] * (J[4] * J[8] -J [5] * J[7]);
					l = J[1] * (J[3] * J[8] -J [5] * J[6]);
					m = J[2] * (J[3] * J[7] -J [4] * J[6]);
					J_det = w - l + m;
					// cout << " Determinent equals " << J_det << endl;
					if (J_det == 0.0)
					{
						std::cout << "As Determinent=0 so it is singular matrix and its inverse cannot exist for element "
							<< std::endl;

						exit(0);
					}
					// transpose of jacobian
					J_T[0]=J[0]; J_T[1]=J[3]; J_T[2]=J[6];
					J_T[3]=J[1]; J_T[4]=J[4]; J_T[5]=J[7];
					J_T[6]=J[2]; J_T[7]=J[5]; J_T[8]=J[8];

					// inverse of K tensor 
					for ( j = 0; j < 3; j++ )
					{  

						for ( i = 0; i < 3; i++ )
						{
							K_inv[i+j*3]=K_inverse (i,j);
						}

					}

					for ( i = 0; i < 3; i++) { // row number of output
						for ( j = 0; j < 3; j++) { // column number of output
							K_inv_ref1[3*i+j] = 0.0;
							for ( k = 0; k < 3; k++) { // three elements are added for this output
								K_inv_ref1[3*i+j] += J[3*i+k] * K_inv[3*k+j];
							}
						}
					}

					for ( i = 0; i < 16; i++) 
					{ 
						B_RT0[i]=0.0;
					}

					for ( i = 0; i < 4; i++) { // row number of output

						for ( j = 0; j < 4; j++) { // column number of output


							for ( Z = 0; Z < 4; Z++) { // GAUSS POINT IN Z DIRECTION 
								x=gauss_pt_xyz[Z*3];
								y=gauss_pt_xyz[Z*3+1];		
								z=gauss_pt_xyz[Z*3+2];



								W_RT0[0]=  2*x;///(sqrt(3.0));
								W_RT0[1]=  2*y;///(sqrt(3.0));
								W_RT0[2]=  2*z;///(sqrt(3.0));

								W_RT0[3]=  2*(-1+x);
								W_RT0[4]=  2*y;
								W_RT0[5]=  2*z;  

								W_RT0[6]=  (2*x);///(sqrt(3.0));
								W_RT0[7]=  2*(-1+y);;///(sqrt(3.0));
								W_RT0[8]=  (2*z);///(sqrt(3.0));

								W_RT0[9]=  2*x;
								W_RT0[10]= 2*y; 
								W_RT0[11]= 2*(-1+z);





								xyz_i[0]=W_RT0[i*3]*J[0]+W_RT0[i*3+1]*J[3]+W_RT0[i*3+2]*J[6];
								xyz_i[1]=W_RT0[i*3]*J[1]+W_RT0[i*3+1]*J[4]+W_RT0[i*3+2]*J[7];
								xyz_i[2]=W_RT0[i*3]*J[2]+W_RT0[i*3+1]*J[5]+W_RT0[i*3+2]*J[8];



								xyz_j[0]=W_RT0[j*3]*K_inv_ref1[0]+W_RT0[j*3+1]*K_inv_ref1[3]+W_RT0[j*3+2]*K_inv_ref1[6];
								xyz_j[1]=W_RT0[j*3]*K_inv_ref1[1]+W_RT0[j*3+1]*K_inv_ref1[4]+W_RT0[j*3+2]*K_inv_ref1[7];
								xyz_j[2]=W_RT0[j*3]*K_inv_ref1[2]+W_RT0[j*3+1]*K_inv_ref1[5]+W_RT0[j*3+2]*K_inv_ref1[8];


Kirill Terekhov's avatar
Kirill Terekhov committed
380
								B_RT0[i+j*4]+=gauss_wei[Z]*(xyz_i[0]*xyz_j[0]+ xyz_i[1]*xyz_j[1]+xyz_i[2]*xyz_j[2])/std::abs(J_det);
ahmadabushaikha's avatar
ahmadabushaikha committed
381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467


							}

						}
					}	
					for ( j = 0; j < 4; j++) { // row number of output
						double sum = 0;
						for ( i = 0; i < 4; i++) { // column number of output

							sum += W(i,j);
							W(i,j)= B_RT0[i+j*4];
						}
						//std::cout << "row " << j << " sum " << sum << std::endl;
					}


					//std::cout << "W" << std::endl;
					//W.Print();

					W = W.Invert(true).first ;


					for ( j = 0; j < 4; j++) { // row number of output
						double sum = 0;
						for ( i = 0; i < 4; i++) { // column number of output

							sum += W(i,j);
						}
						//std::cout << "row inverted " << j << " sum " << sum << std::endl;
					}

					//std::cout << "W inverted" << std::endl;

					//W.Print();

					

					//scanf("%*c");

				}
				else
				{
					real xP[3]; //center of the cell
					real yF[3]; //center of the face
					real nF[3]; //normal to the face
					real aF; //area of the face
					real vP = cell->Volume(); //volume of the cell
					cell->Centroid(xP);
					rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),cell->RealArrayDF(tag_K).size()); //get permeability for the cell
					//rMatrix U,S,V;
					//K0.SVD(U,S,V);
					//for(int k = 0; k < 3; ++k) S(k,k) = sqrt(S(k,k));
					//rMatrix K = U*S*V;
					rMatrix NK(NF,3), R(NF,3), D(NF,NF), U(NF,NF), Areas(NF,NF); //big gradient matrix, co-normals, directions
					Areas.Zero();
					for(int k = 0; k < NF; ++k) //loop over faces
					{
						aF = faces[k].Area();
						faces[k].Centroid(yF);
						faces[k].OrientedUnitNormal(cell->self(),nF);
						// assemble matrix of directions
						R(k,0) = (yF[0]-xP[0])*aF;
						R(k,1) = (yF[1]-xP[1])*aF;
						R(k,2) = (yF[2]-xP[2])*aF;
						// assemble matrix of co-normals
						rMatrix nK = rMatrix::FromVector(nF,3).Transpose()*K;
						NK(k,0) = nK(0,0);
						NK(k,1) = nK(0,1);
						NK(k,2) = nK(0,2);

						Areas(k,k) = aF;
					} //end of loop over faces
					rMatrix SU,SS,SV;
					W = NK*(NK.Transpose()*R).Invert(true).first*NK.Transpose(); //stability part
					/*
					 std::cout << "W" << std::endl;
					 nKGRAD.Print();
					 nKGRAD.SVD(SU,SS,SV);
					 std::cout << "U" << std::endl;
					 SU.Print();
					 std::cout << "S" << std::endl;
					 SS.Print();
					 std::cout << "V" << std::endl;
					 SV.Print();
					 std::cout << "Check " << (nKGRAD - SU*SS*SV.Transpose()).FrobeniusNorm() << std::endl;
					 */
Kirill Terekhov's avatar
Kirill Terekhov committed
468
                
ahmadabushaikha's avatar
ahmadabushaikha committed
469
					int rank = 0; //size of matrix U
Kirill Terekhov's avatar
Kirill Terekhov committed
470
                
ahmadabushaikha's avatar
ahmadabushaikha committed
471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500
					{ //Retrive orthogonal to R matrix D
						//Symmetric orthogonal matrix
						rMatrix DUD = (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert(true).first*R.Transpose());
						//perfrom singular value decomposition
						//S should be unity matrix with rank NF-3
						rMatrix DUD_U,DUD_S,DUD_V;
						DUD.SVD(DUD_U,DUD_S,DUD_V);
						//compute the rank
						for(int q = 0; q < NF; ++q)
							if( DUD_S(q,q) > 1.0e-2 )
								rank++;
						rank = NF-3;
						if( rank != NF-3)
						{
							std::cout << "rank: " << rank << " expected " << NF-3 << std::endl;
							DUD_S.Print();
						}
						//chop matrix to the full rank
						DUD_S.RemoveSubset(rank,NF,rank,NF);
						DUD_V.RemoveColumns(rank,NF);
						//assign the matrix
						D = DUD_V;
						U = DUD_S;
					}
					//std::cout << "D" << std::endl;
					//D.Print();
					//std::cout << "U" << std::endl;
					//U.Print();
					//std::cout << "DtR" << std::endl;
					//(D.Transpose()*R).Print();
Kirill Terekhov's avatar
Kirill Terekhov committed
501
                
ahmadabushaikha's avatar
ahmadabushaikha committed
502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531
					U *=(2.0/(static_cast<real>(NF)*vP)*(NK*K.Invert(true).first*NK.Transpose()).Trace());

					//std::cout << "UDtR" << std::endl;
					//(U*D.Transpose()*R).Print();
                
					W += D*U*D.Transpose();


					W = Areas*W*Areas;
				}
					//std::cout << "W: " << std::endl;
					//W.Print();
				bulk & isDMP = cell->Bulk(tag_DMP);
				isDMP = 1;
				for(int k = 0; k < NF; ++k)
				{
					real row_sum = 0;
					if( W(k,k) < 0.0 ) isDMP = 0;
					for(int j = 0; j < NF; ++j)
						row_sum += W(k,j);
					if( row_sum < 0.0 ) isDMP = 0;
					for(int j = k+1; j < NF; ++j)
						if( W(k,j) > 0.0 )
							isDMP = 0;
				}
				++total;
				if( isDMP ) ++dmp;
				real_array store_W = cell->RealArrayDV(tag_W); //access data structure for gradient matrix in mesh
				store_W.resize(NF*NF); //resize the structure
				std::copy(W.data(),W.data()+NF*NF,store_W.data()); //write down the gradient matrix
Kirill Terekhov's avatar
Kirill Terekhov committed
532 533
            } //end of loop over cells
            std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
534
            std::cout << "Satisfy DMP: " << dmp << " out of " << total << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
535 536 537 538 539 540 541 542 543 544 545 546 547 548 549
            
            if( m->HaveTag("FORCE") ) //Is there force on the mesh?
            {
                tag_F = m->GetTag("FORCE"); //initial force
                assert(tag_F.isDefined(CELL)); //assuming it was defined on cells
            } // end of force
        } //end of initialize data
        
        std::cout << "Initialization done" << std::endl;
        
        
        integer nit = 0;
        ttt = Timer();
        
        { //Main loop for problem solution
Kirill Terekhov's avatar
Kirill Terekhov committed
550
            Automatizator aut; // declare class to help manage unknowns
Kirill Terekhov's avatar
Kirill Terekhov committed
551
            Automatizator::MakeCurrent(&aut);
Kirill Terekhov's avatar
Kirill Terekhov committed
552 553
            dynamic_variable P(aut,aut.RegisterTag(tag_P,CELL|FACE)); //register pressure as primary unknown
            aut.EnumerateTags(); //enumerate all primary variables
Kirill Terekhov's avatar
Kirill Terekhov committed
554
            
ahmadabushaikha's avatar
ahmadabushaikha committed
555
            std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579
            
            
            Residual R("",aut.GetFirstIndex(),aut.GetLastIndex());
            Sparse::LockService Locks(aut.GetFirstIndex(),aut.GetLastIndex());
            Sparse::AnnotationService Text(aut.GetFirstIndex(),aut.GetLastIndex());
            Sparse::Vector Update  ("",aut.GetFirstIndex(),aut.GetLastIndex()); //vector for update
            {//Annotate matrix
                for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
                {
                    Cell cell = m->CellByLocalID(q);
                    if( cell.GetStatus() != Element::Ghost )
                        Text.SetAnnotation(P.Index(cell),"Cell-centered pressure value");
                }
                for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                {
                    Face face = m->FaceByLocalID(q);
                    if( face.GetStatus() != Element::Ghost )
                    {
                        if( tag_BC.isValid() && face.HaveData(tag_BC) )
                            Text.SetAnnotation(P.Index(face),"Pressure guided by boundary condition");
                        else
                            Text.SetAnnotation(P.Index(face),"Interface pressure");
                    }
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
580
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
581 582 583 584 585 586 587
            
            std::cout << "Matrix was annotated" << std::endl;
            
            do
            {
                R.Clear(); //clean up the residual
                double tttt = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
588
                int total = 0, dmp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
589
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
590
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
Kirill Terekhov's avatar
Kirill Terekhov committed
591
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
592 593 594 595 596 597
                for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) //loop over cells
                {
                    Cell cell = m->CellByLocalID(q);
                    ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
                    int NF = (int)faces.size();
                    rMatrix nKGRAD(cell->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
Kirill Terekhov's avatar
Kirill Terekhov committed
598
                    bulk & isDMP = cell->Bulk(tag_DMP);
Kirill Terekhov's avatar
Kirill Terekhov committed
599 600
                    vMatrix pF(NF,1); //vector of pressure differences on faces
                    vMatrix FLUX(NF,1); //computed flux on faces
ahmadabushaikha's avatar
ahmadabushaikha committed
601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663
                    vMatrix FLUX_MASS(NF,1); //computed flux on faces for mass balance on cell
					vMatrix POT (2,1); //vector of pressure differences on faces

					for(int k = 0; k < NF; ++k)
					{
					
					
					if (faces[k]->FrontCell().isValid())
					{
						Cell cell_n = cell.Neighbour(faces[k]); 
						ElementArray<Face> faces_n = cell_n->getFaces(); //obtain faces of the cell neighoubr
					    rMatrix nKGRAD_n(cell_n->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
						double B_1, B_2, L_1, L_2;
						int face_n_k;

                        L_1 =0.0; 
						B_1= nKGRAD(k,k) ;
						for(int j = 0; j < NF; ++j)
						L_1 +=nKGRAD(k,j) ;
						for(int j = 0; j < NF; ++j) if (faces_n[j] == faces[k])
						{
							face_n_k= j ;
						}
						L_2 =0.0; 
						for(int j = 0; j < NF; ++j)
						L_2 +=nKGRAD_n(face_n_k,j) ;

						B_2 = nKGRAD_n(face_n_k,face_n_k) ;

						POT(0,0)=L_1*P(cell) ;
						POT(1,0)=L_2*P(cell_n) ;

						for(int j = 0; j < NF; ++j)if (j != k)
						POT(0,0) -= nKGRAD(k,j)* P(faces[j]);

						for(int j = 0; j < NF; ++j) if (j != face_n_k)
						POT(1,0) -= nKGRAD_n(face_n_k,j)* P(faces_n[j]);

                        FLUX_MASS(k,0) = -(B_2*POT(0,0)-B_1 *POT(1,0))/(B_1+B_2) ;


					}
					else {
						FLUX_MASS(k,0)=0.0;
						
					for(int j = 0; j < NF; ++j)
						FLUX_MASS(k,0)+= nKGRAD(k,j)* (P(faces[j]) - P(cell)) ;
					}
					}
				

				//	for(int k = 0; k < NF; ++k)
				//	{
				//	FLUX_MASS(k,0)=0.0;
				////	Cell cell_l =faces[k].BackCell();
				//	Cell cell_l= cell.Neighbour(faces[k]);
				//	Cell cell_n= cell_l.Neighbour(faces[k]);
				//	for(int j = 0; j < NF; ++j)
				//	FLUX_MASS(k,0)+= nKGRAD(k,j)* (P(faces[j]) - P(cell_n)) ;
				//	}
					
	                    for(int k = 0; k < NF; ++k)
                    pF(k,0) = (P(faces[k]) - P(cell));//*faces[k].Area();
Kirill Terekhov's avatar
Kirill Terekhov committed
664
                    FLUX = nKGRAD*pF; //fluxes on faces
Kirill Terekhov's avatar
Kirill Terekhov committed
665 666
                    //Change matrix by nonlinear correction for DMP
                    ++total;
Kirill Terekhov's avatar
Kirill Terekhov committed
667 668 669
                    if( cell.GetStatus() != Element::Ghost )
                    {
                        for(int k = 0; k < NF; ++k) //loop over faces of current cell
ahmadabushaikha's avatar
ahmadabushaikha committed
670 671 672
                        { //R[P.Index(cell)] += FLUX(k,0);//faces[k].Area();
						   R[P.Index(cell)] +=  FLUX_MASS(k,0);//faces[k].Area();
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689
                    }
                    for(int k = 0; k < NF; ++k) //loop over faces of current cell
                    {
                        if( faces[k].GetStatus() == Element::Ghost ) continue;
                        int index = P.Index(faces[k]);
                        Locks.Lock(index);
                        if( tag_BC.isValid() && faces[k].HaveData(tag_BC) )
                        {
                            real_array BC = faces[k].RealArray(tag_BC);
                            R[index] -= BC[0]*P(faces[k]) + BC[1]*FLUX(k,0) - BC[2];
                        }
                        else
                            R[index] -= FLUX(k,0);
                        Locks.UnLock(index);
                    }
                } //end of loop over cells
                
Kirill Terekhov's avatar
Kirill Terekhov committed
690 691
                std::cout << "Satisfy DMP: " << dmp << " out of " << total << std::endl;
                
Kirill Terekhov's avatar
Kirill Terekhov committed
692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709
                
                if( tag_F.isValid() )
                {
#if defined(USE_OMP)
#pragma omp parallel for
#endif
                    for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
                    {
                        Cell cell = m->CellByLocalID(q);
                        if( cell.GetStatus() == Element::Ghost ) continue;
                        if( cell->HaveData(tag_F) ) R[P.Index(cell)] += cell->Real(tag_F)*cell->Volume();
                    }
                }
                
                std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl;
                
                
                R.Rescale();
Kirill Terekhov's avatar
Kirill Terekhov committed
710 711
                //R.GetJacobian().Save("jacobian.mtx",&Text);
                //R.GetResidual().Save("residual.mtx");
Kirill Terekhov's avatar
Kirill Terekhov committed
712 713 714 715 716 717
                
                
                std::cout << "Nonlinear residual: " << R.Norm() << "\t\t" << std::endl;
                
                if( R.Norm() < 1.0e-4 ) break;
                
ahmadabushaikha's avatar
ahmadabushaikha committed
718
				//Solver S(Solver::INNER_ILU2);
719 720
                //Solver S(Solver::INNER_MPTILUC);
				Solver S(Solver::SUPERLU);
Kirill Terekhov's avatar
Kirill Terekhov committed
721
				S.SetParameterReal("relative_tolerance", 1.0e-14);
Kirill Terekhov's avatar
Kirill Terekhov committed
722
                S.SetParameterReal("absolute_tolerance", 1.0e-12);
ahmadabushaikha's avatar
ahmadabushaikha committed
723 724
                S.SetParameterReal("drop_tolerance", 1.0e-1);
                S.SetParameterReal("reuse_tolerance", 1.0e-2);
Kirill Terekhov's avatar
Kirill Terekhov committed
725
                S.SetMatrix(R.GetJacobian());
Kirill Terekhov's avatar
Kirill Terekhov committed
726 727 728
                //std::fill(Update.Begin(),Update.End(),0.0);
                if( S.Solve(R.GetResidual(),Update) )
                {
Kirill Terekhov's avatar
Kirill Terekhov committed
729 730 731
#if defined(USE_OMP)
#pragma omp parallel for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
732 733 734 735 736
                    for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
                    {
                        Cell cell = m->CellByLocalID(q);
                        cell->Real(tag_P) -= Update[P.Index(cell)];
                    }
Kirill Terekhov's avatar
Kirill Terekhov committed
737 738 739
#if defined(USE_OMP)
#pragma omp parallel for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762
                    for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                    {
                        Face face = m->FaceByLocalID(q);
                        face->Real(tag_P) -= Update[P.Index(face)];
                    }
                    m->ExchangeData(tag_P, CELL|FACE, 0);
                    {
                        std::stringstream str;
                        str << "iter" << nit;
                        if( m->GetProcessorsNumber() == 1 )
                            str << ".vtk";
                        else
                            str << ".pvtk";
                        m->Save(str.str());
                    }
                }
                else
                {
                    std::cout << "Unable to solve: " << S.GetReason() << std::endl;
                    break;
                }
                ++nit;
            } while( R.Norm() > 1.0e-4 && nit < 10); //check the residual norm
Kirill Terekhov's avatar
Kirill Terekhov committed
763
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
764 765 766
        std::cout << "Solved problem in " << Timer() - ttt << " seconds with " << nit << " iterations " << std::endl;
        
        if( m->HaveTag("REFERENCE_SOLUTION") )
Kirill Terekhov's avatar
Kirill Terekhov committed
767
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801
            Tag tag_E = m->CreateTag("ERRROR",DATA_REAL,CELL,NONE,1);
            Tag tag_R = m->GetTag("REFERENCE_SOLUTION");
            real C, L2, volume;
            C = L2 = volume = 0.0;
            for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
            {
                Cell cell = m->CellByLocalID(q);
                real err = cell->Real(tag_P) - cell->Real(tag_R);
                real vol = cell->Volume();
                if( C < fabs(err) ) C = fabs(err);
                L2 += err*err*vol;
                volume += vol;
                cell->Real(tag_E) = err;
            }
            L2 = sqrt(L2/volume);
            std::cout << "Error on cells, C-norm " << C << " L2-norm " << L2 << std::endl;
            C = L2 = volume = 0.0;
            if( tag_R.isDefined(FACE) )
            {
                tag_E = m->CreateTag("ERRROR",DATA_REAL,FACE,NONE,1);
                for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                {
                    Face face = m->FaceByLocalID(q);
                    real err = face->Real(tag_P) - face->Real(tag_R);
                    real vol = (face->BackCell()->Volume() + (face->FrontCell().isValid() ? face->FrontCell()->Volume() : 0))*0.5;
                    if( C < fabs(err) ) C = fabs(err);
                    L2 += err*err*vol;
                    volume += vol;
                    face->Real(tag_E) = err;
                }
                L2 = sqrt(L2/volume);
                std::cout << "Error on faces, C-norm " << C << " L2-norm " << L2 << std::endl;
            }
            else std::cout << "Reference solution was not defined on faces" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
802
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
803 804 805 806 807 808 809
        
        if( m->GetProcessorsNumber() == 1 )
            m->Save("out.vtk");
        else
            m->Save("out.pvtk");
        
        delete m; //clean up the mesh
Kirill Terekhov's avatar
Kirill Terekhov committed
810
    }
Kirill Terekhov's avatar
Kirill Terekhov committed
811
    else
Kirill Terekhov's avatar
Kirill Terekhov committed
812
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
813
        std::cout << argv[0] << " mesh_file" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
814
    }
Kirill Terekhov's avatar
Kirill Terekhov committed
815
    
Kirill Terekhov's avatar
Kirill Terekhov committed
816
#if defined(USE_PARTITIONER)
Kirill Terekhov's avatar
Kirill Terekhov committed
817
    Partitioner::Finalize(); // Finalize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
818
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
819 820
    Solver::Finalize(); // Finalize solver and close MPI activity
    return 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
821
}