--- /dev/null
+//============================================================================
+// Author : Michael NDJINGA
+// Date : November 2020
+// Description : 2D linear wave system
+//============================================================================
+
+#include <iostream>
+#include <string>
+#include <cmath>
+
+#include "Mesh.hxx"
+#include "Cell.hxx"
+#include "Face.hxx"
+#include "Field.hxx"
+#include "SparseMatrixPetsc.hxx"
+#include "CdmathException.hxx"
+
+using namespace std;
+
+double p0 =155e5; //reference pressure in a pressurised nuclear vessel
+double c0 =700.; //reference sound speed for water at 155 bars, 600K
+double rho0=p0/c0*c0;//reference density
+double precision=1e-5;
+
+void initial_conditions_shock(Mesh my_mesh,Field& pressure_field,Field& velocity_field)
+{
+ double rayon=0.35;
+ double xcentre=0.;
+ double ycentre=0;
+
+ int dim =my_mesh.getMeshDimension();
+ int nbCells=my_mesh.getNumberOfCells();
+
+ for (int j=0 ; j<nbCells ; j++)
+ {
+ double x = my_mesh.getCell(j).x() ;
+ double y = my_mesh.getCell(j).y() ;
+ double valX=(x-xcentre)*(x-xcentre);
+ double valY=(y-ycentre)*(y-ycentre);
+ double val=sqrt(valX+valY);
+
+ for(int idim=0; idim<dim; idim++)
+ velocity_field[j,idim]=0;
+
+ if (val<rayon)
+ pressure_field(j) = 155e5;
+ else
+ pressure_field(j) = 70e5;
+ }
+}
+
+Matrix jacobianMatrices(Vector normal, double coeff)
+{
+ int dim=normal.size();
+ Matrix A(dim+1,dim+1);
+ Matrix absA(dim+1,dim+1);
+
+ absA(0,0)=c0*coeff;
+ for(int i=0 ; i<dim; i++)
+ {
+ A(i+1,0)= normal[i]*coeff;
+ A(0,i+1)=c0*c0*normal[i]*coeff;
+ for( int j=0 ; j<dim; j++)
+ absA(i+1,j+1)=c0*normal[i]*normal[j]*coeff;
+ }
+ return (A - absA)*(1./2);
+}
+
+SparseMatrixPetsc computeDivergenceMatrix(Mesh my_mesh, int nbVoisinsMax, double dt)
+{
+ int nbCells = my_mesh.getNumberOfCells();
+ int dim=my_mesh.getMeshDimension();
+ int nbComp=dim+1;
+ Vector normal(dim);
+
+ SparseMatrixPetsc implMat(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp);
+
+ Matrix idMoinsJacCL(nbComp);
+
+ for(int j=0; j<nbCells; j++)//On parcourt les cellules
+ {
+ Cell Cj = my_mesh.getCell(j);
+ int nbFaces = Cj.getNumberOfFaces();
+
+ for(int k=0; k<nbFaces; k++)
+ {
+ int indexFace = Cj.getFacesId()[k];
+ Face Fk = my_mesh.getFace(indexFace);
+ for(int i =0; i<dim ; i++)
+ normal[i] = Cj.getNormalVector(k, i);//normale sortante
+
+ Matrix Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
+
+ int cellAutre =-1;
+ if ( not Fk.isBorder())
+ {
+ /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
+ if (Fk.getCellsId()[0] == j)
+ // hypothese verifiée
+ cellAutre = Fk.getCellsId()[1];
+ else if(Fk.getCellsId()[1] == j)
+ // hypothese non verifiée
+ cellAutre = Fk.getCellsId()[0];
+ else
+ throw CdmathException("computeDivergenceMatrix: problem with mesh, unknown cell number");
+
+ implMat.addValue(j*nbComp,cellAutre*nbComp,Am);
+ implMat.addValue(j*nbComp, j*nbComp,Am*(-1.));
+ }
+ else
+ {
+ if( Fk.getGroupName() != "Periodic" && Fk.getGroupName() != "Neumann")//Wall boundary condition unless Periodic/Neumann specified explicitly
+ {
+ Vector v(dim+1);
+ for(int i=0; i<dim; i++)
+ v[i+1]=normal[i];
+ Matrix idMoinsJacCL=v.tensProduct(v)*2;
+
+ implMat.addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL);
+ }
+ else if( Fk.getGroupName() == "Periodic")//Periodic boundary condition
+ {
+ int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
+ Face Fp = my_mesh.getFace(indexFP);
+ cellAutre = Fp.getCellsId()[0];
+
+ implMat.addValue(j*nbComp,cellAutre*nbComp,Am);
+ implMat.addValue(j*nbComp, j*nbComp,Am*(-1.));
+ }
+ else if(Fk.getGroupName() != "Neumann")//Nothing to do for Neumann boundary condition
+ {
+ cout<< Fk.getGroupName() <<endl;
+ throw CdmathException("computeDivergenceMatrix: Unknown boundary condition name");
+ }
+ }
+ }
+ }
+ return implMat;
+}
+
+void WaveSystem2D(double tmax, double ntmax, double cfl, int output_freq, const Mesh& my_mesh, const string file)
+{
+ /* Retrieve mesh data */
+ int dim=my_mesh.getMeshDimension();
+ int nbCells = my_mesh.getNumberOfCells();
+ std::string meshName="HEXAGON800";//my_mesh.getName();
+ int nbVoisinsMax=my_mesh.getMaxNbNeighbours(CELLS);
+ double dx_min=my_mesh.minRatioVolSurf();
+
+
+ /* Initial conditions */
+ cout << "Construction of the initial condition …" << endl;
+
+ Field pressure_field("Pressure",CELLS,my_mesh,1) ;
+ Field velocity_field("Velocity",CELLS,my_mesh,1) ;
+ initial_conditions_shock(my_mesh,pressure_field, velocity_field);
+
+ /* iteration vectors */
+ Vector Un(nbCells*(dim+1));
+ Vector dUn(nbCells*(dim+1));
+
+ for(int k =0; k<nbCells; k++){
+ Un[k*(dim+1)+0] = pressure_field[k];
+ Un[k*(dim+1)+1] =rho0*velocity_field[k,0];
+ if(dim>=2)
+ Un[k + 2*nbCells] = rho0*velocity_field[k,1] ;
+ if(dim==3)
+ Un[k + 3*nbCells] = rho0*velocity_field[k,2];
+ }
+
+ /*
+ * MED output of the initial condition at t=0 and iter = 0
+ */
+ int it=0;
+ bool isStationary=false;
+ double time=0.;
+ double dt = cfl * dx_min / c0;
+
+ cout << "Saving the solution at T=" << time << "…" << endl;
+ pressure_field.setTime(time,it);
+ pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure");
+ velocity_field.setTime(time,it);
+ velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity");
+ /* --------------------------------------------- */
+
+ SparseMatrixPetsc divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt);
+
+ /* Time loop */
+ cout<< "Starting computation of the linear wave system with an explicit UPWIND scheme …" << endl;
+ while (it<ntmax && time <= tmax )
+ {
+ dUn=divMat*Un;
+ Un-=dUn;
+
+ time=time+dt;
+ it=it+1;
+
+ isStationary = dUn.norm()<precision;
+ /* Sauvegardes */
+ if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax)
+ {
+ cout<<"-- Iteration: " << it << ", Time: " << time << ", dt: " << dt<<endl;
+
+ for(int k=0; k<nbCells; k++)
+ {
+ pressure_field[k] =Un[k*(dim+1)+0];
+ velocity_field[k,0]=Un[k*(dim+1)+1]/rho0;
+ if(dim>1)
+ {
+ velocity_field[k,1]=Un[k*(dim+1)+2]/rho0;
+ if(dim>2)
+ velocity_field[k,2]=Un[k*(dim+1)+3]/rho0;
+ }
+ }
+ pressure_field.setTime(time,it);
+ pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure",false);
+ velocity_field.setTime(time,it);
+ velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity",false);
+ }
+ }
+ cout<<"End of calculation -- Iteration: " << it << ", Time: "<< time<< ", dt: " << dt<<endl;
+
+ if(it>=ntmax)
+ cout<< "Nombre de pas de temps maximum ntmax= "<< ntmax<< " atteint"<<endl;
+ else if(isStationary)
+ cout<< "Régime stationnaire atteint au pas de temps "<< it<< ", t= "<< time<<endl;
+ else
+ cout<< "Temps maximum Tmax= "<< tmax<< " atteint"<<endl;
+}
+
+int main(int argc, char *argv[])
+{
+ cout << "RESOLUTION OF THE 2D WAVE SYSTEM: Upwind explicit scheme" << endl;
+ cout << "- WALL BC" << endl;
+
+ // Problem data
+ double cfl=0.49;
+ double tmax=1.;
+ double ntmax=20000;
+ int freqSortie=100;
+ string fileOutPut="SphericalWave";
+
+ if(argc<2)
+ {
+ cout << "- DOMAIN: SQUARE" << endl;
+ cout << "- MESH: CARTESIAN, GENERATED INTERNALLY WITH CDMATH" << endl;
+ cout << "Construction of a cartesian mesh …" << endl;
+ double xinf=0.0;
+ double xsup=1.0;
+ double yinf=0.0;
+ double ysup=1.0;
+ int nx=300;
+ int ny=300;
+ Mesh myMesh(xinf,xsup,nx,yinf,ysup,ny);
+ double eps=1.E-10;
+ myMesh.setGroupAtPlan(xsup,0,eps,"RightEdge");
+ myMesh.setGroupAtPlan(xinf,0,eps,"LeftEdge");
+ myMesh.setGroupAtPlan(yinf,1,eps,"BottomEdge");
+ myMesh.setGroupAtPlan(ysup,1,eps,"TopEdge");
+
+ WaveSystem2D(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut);
+ }
+ else
+ {
+ cout << "- MESH: GENERATED EXTERNALLY WITH SALOME" << endl;
+ cout << "Loading of a mesh …" << endl;
+ string filename = argv[1];
+ Mesh myMesh(filename);
+
+ WaveSystem2D(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut);
+ }
+
+ cout << "Simulation complete." << endl;
+
+ return 0;
+}
--- /dev/null
+#!/usr/bin/env python
+# -*-coding:utf-8 -*
+
+from math import sin, cos, pi, sqrt
+import cdmath
+import PV_routines
+import VTK_routines
+import sys
+
+c0=700.#reference sound speed
+p0=155e5#reference pressure
+rho0=p0/c0*c0#reference density
+precision=1e-5
+
+def initial_conditions_wave_system(my_mesh):
+ print "Spherical wave initial data"
+ dim = my_mesh.getMeshDimension()
+ nbCells = my_mesh.getNumberOfCells()
+
+ rayon = 0.15
+ xcentre = 0.5
+ ycentre = 0.5
+ zcentre = 0.5
+
+ pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
+ velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
+
+ for i in range(nbCells):
+ x = my_mesh.getCell(i).x()
+ y = my_mesh.getCell(i).y()
+
+ velocity_field[i,0] = 0
+ velocity_field[i,1] = 0
+ velocity_field[i,2] = 0
+
+ valX = (x - xcentre) * (x - xcentre)
+ valY = (y - ycentre) * (y - ycentre)
+ if(dim==2):
+ val = sqrt(valX + valY)
+ if(dim==3):
+ z = my_mesh.getCell(i).z()
+ valZ = (z - zcentre) * (z - zcentre)
+ val = sqrt(valX + valY + valZ)
+
+ if val < rayon:
+ pressure_field[i] = p0
+ pass
+ else:
+ pressure_field[i] = p0/2
+ pass
+ pass
+
+ return pressure_field, velocity_field
+
+def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling):
+ nbCells = my_mesh.getNumberOfCells()
+ dim=my_mesh.getMeshDimension()
+ nbComp=dim+1
+
+ if(not my_mesh.isStructured()):
+ raise ValueError("WaveSystemStaggered: the mesh should be structured");
+
+ NxNyNz=my_mesh.getCellGridStructure()
+ DxDyDz=my_mesh.getDXYZ()
+
+ implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
+
+ idMoinsJacCL=cdmath.Matrix(nbComp)
+
+ if( dim == 1) :
+ nx=NxNyNz[0]
+ dx=DxDyDz[0]
+
+ if( scaling==0 ):
+ for k in range(nbCells):
+ implMat.addValue(k,1*nbCells + k , -c0*c0*dt/dx)
+ implMat.addValue(k,1*nbCells + (k+1)%nx, c0*c0*dt/dx)
+
+ implMat.addValue( 1*nbCells + k ,k, dt/dx)
+ implMat.addValue( 1*nbCells + (k+1)%nx,k, -dt/dx)
+ else : # scaling >0
+ for k in range(nbCells):
+ implMat.addValue(k,1*nbCells + k , -c0*dt/dx)
+ implMat.addValue(k,1*nbCells + (k+1)%nx, c0*dt/dx)
+
+ implMat.addValue( 1*nbCells + k ,k, c0*dt/dx)
+ implMat.addValue( 1*nbCells + (k+1)%nx,k, -c0*dt/dx)
+
+ elif( dim == 2) :# k = j*nx+i
+ nx=NxNyNz[0]
+ ny=NxNyNz[1]
+ dx=DxDyDz[0]
+ dy=DxDyDz[1]
+
+ if( scaling==0 ):
+ for k in range(nbCells):
+ i = k % nx
+ j = k //nx
+
+ implMat.addValue(k,1*nbCells + j*nx + i , -c0*c0*dt/dx)
+ implMat.addValue(k,1*nbCells + j*nx + (i+1)%nx, c0*c0*dt/dx)
+
+ implMat.addValue(k,2*nbCells + j *nx + i, -c0*c0*dt/dy)
+ implMat.addValue(k,2*nbCells + ((j+1)%ny)*nx + i, c0*c0*dt/dy)
+
+ implMat.addValue( 1*nbCells + j*nx + i , k, dt/dx)
+ implMat.addValue( 1*nbCells + j*nx + (i+1)%nx, k, -dt/dx)
+
+ implMat.addValue( 2*nbCells + j *nx + i,k, dt/dy)
+ implMat.addValue( 2*nbCells + ((j+1)%ny)*nx + i,k, -dt/dy)
+
+ else :# scaling >0
+ for k in range(nbCells):
+ i = k % nx
+ j = k //nx
+
+ implMat.addValue(k,1*nbCells + j*nx + i , -c0*dt/dx)
+ implMat.addValue(k,1*nbCells + j*nx + (i+1)%nx, c0*dt/dx)
+
+ implMat.addValue(k,2*nbCells + j *nx + i, -c0*dt/dy)
+ implMat.addValue(k,2*nbCells + ((j+1)%ny)*nx + i, c0*dt/dy)
+
+ implMat.addValue( 1*nbCells + j*nx + i , k, c0*dt/dx)
+ implMat.addValue( 1*nbCells + j*nx + (i+1)%nx, k, -c0*dt/dx)
+
+ implMat.addValue( 2*nbCells + j *nx + i,k, c0*dt/dy)
+ implMat.addValue( 2*nbCells + ((j+1)%ny)*nx + i,k, -c0*dt/dy)
+
+ elif( dim == 3) :# k = l*nx*ny+j*nx+i
+ nx=NxNyNz[0]
+ ny=NxNyNz[1]
+ nz=NxNyNz[2]
+ dx=DxDyDz[0]
+ dy=DxDyDz[1]
+ dz=DxDyDz[2]
+
+ if( scaling==0 ):
+ for k in range(nbCells):
+ i = k % nx
+ j = (k //nx)%ny
+ l = k //(nx*ny)
+
+ implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + i , -c0*c0*dt/dx)
+ implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + (i+1)%nx, c0*c0*dt/dx)
+
+ implMat.addValue(k,2*nbCells + l*nx*ny + j *nx + i, -c0*c0*dt/dy)
+ implMat.addValue(k,2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i, c0*c0*dt/dy)
+
+ implMat.addValue(k,3*nbCells + l*nx*ny + j*nx + i, -c0*c0*dt/dz)
+ implMat.addValue(k,3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i, c0*c0*dt/dz)
+
+ implMat.addValue( 1*nbCells + l*nx*ny + j*nx + i , k, dt/dx)
+ implMat.addValue( 1*nbCells + l*nx*ny + j*nx + (i+1)%nx, k, -dt/dx)
+
+ implMat.addValue( 2*nbCells + l*nx*ny + j *nx + i,k, dt/dy)
+ implMat.addValue( 2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i,k, -dt/dy)
+
+ implMat.addValue( 3*nbCells + l*nx*ny + j*nx + i,k, dt/dz)
+ implMat.addValue( 3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i,k, -dt/dz)
+
+ else:# scaling >0
+ for k in range(nbCells):
+ i = k % nx
+ j = (k //nx)%ny
+ l = k //(nx*ny)
+
+ implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + i , -c0*dt/dx)
+ implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + (i+1)%nx, c0*dt/dx)
+
+ implMat.addValue(k,2*nbCells + l*nx*ny + j *nx + i, -c0*dt/dy)
+ implMat.addValue(k,2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i, c0*dt/dy)
+
+ implMat.addValue(k,3*nbCells + l*nx*ny + j*nx + i, -c0*dt/dz)
+ implMat.addValue(k,3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i, c0*dt/dz)
+
+ implMat.addValue( 1*nbCells + l*nx*ny + j*nx + i , k, c0*dt/dx)
+ implMat.addValue( 1*nbCells + l*nx*ny + j*nx + (i+1)%nx, k, -c0*dt/dx)
+
+ implMat.addValue( 2*nbCells + l*nx*ny + j *nx + i,k, c0*dt/dy)
+ implMat.addValue( 2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i,k, -c0*dt/dy)
+
+ implMat.addValue( 3*nbCells + l*nx*ny + j*nx + i,k, c0*dt/dz)
+ implMat.addValue( 3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i,k, -c0*dt/dz)
+
+ return implMat
+
+def WaveSystemStaggered(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution):
+ dim=my_mesh.getMeshDimension()
+ nbCells = my_mesh.getNumberOfCells()
+
+ dt = 0.
+ time = 0.
+ it=0;
+ isStationary=False;
+
+ scaling=0
+
+ nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
+ iterGMRESMax=50
+
+ #iteration vectors
+ Un=cdmath.Vector(nbCells*(dim+1))
+ dUn=cdmath.Vector(nbCells*(dim+1))
+
+ # Initial conditions #
+ print("Construction of the initial condition …")
+ pressure_field, velocity_field = initial_conditions_wave_system(my_mesh)
+
+ for k in range(nbCells):
+ Un[k + 0*nbCells] = pressure_field[k]
+ Un[k + 1*nbCells] = rho0*velocity_field[k,0] # value on the left face
+ Un[k + 2*nbCells] = rho0*velocity_field[k,1] # value on the bottom face
+ if(dim==3):
+ Un[k + 3*nbCells] = rho0*initial_velocity[k,2]
+ if( scaling>0):
+ Vn = Un.deepCopy()
+ for k in range(nbCells):
+ Vn[k] = Vn[k]/c0
+
+ #sauvegarde de la donnée initiale
+ pressure_field.setTime(time,it);
+ pressure_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure");
+ velocity_field.setTime(time,it);
+ velocity_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity");
+ #Postprocessing : save 2D picture
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_initial")
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_initial")
+
+ dx_min=my_mesh.minRatioVolSurf()
+
+ dt = cfl * dx_min / c0
+
+ divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling)
+
+ #Add the identity matrix on the diagonal
+ for j in range(nbCells*(dim+1)):
+ divMat.addValue(j,j,1)
+
+ LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","LU")
+
+ print("Starting computation of the linear wave system with an pseudo staggered scheme …")
+
+ # Starting time loop
+ while (it<ntmax and time <= tmax and not isStationary):
+ dUn=Un.deepCopy()
+ LS.setSndMember(Un)
+ Un=LS.solve();
+ cvgceLS=LS.getStatus();
+ iterGMRES=LS.getNumberOfIter();
+ if(not cvgceLS):
+ print "Linear system did not converge ", iterGMRES, " GMRES iterations"
+ raise ValueError("Pas de convergence du système linéaire");
+ dUn-=Un
+
+ max_dp=0 ; max_dq=0
+ for k in range(nbCells):
+ max_dp = max(max_dp,abs(dUn[k]))
+ for i in range(dim):
+ max_dq=max(max_dq,abs(dUn[k+(1+i)*nbCells]))
+
+ isStationary= max_dp/p0<precision and max_dq/rho0<precision
+
+ time=time+dt;
+ it=it+1;
+
+ #Sauvegardes
+ if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
+ print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
+ print "Variation temporelle relative : pressure ", max_dp/p0 ,", velocity ", max_dq/rho0
+ print "Linear system converged in ", iterGMRES, " GMRES iterations"
+
+ for k in range(nbCells):
+ pressure_field[k]=Un[k]
+ velocity_field[k,0]=Un[k+1*nbCells]/rho0
+ if(dim>1):
+ velocity_field[k,1]=Un[k+2*nbCells]/rho0
+ if(dim>2):
+ velocity_field[k,2]=Un[k+3*nbCells]/rho0
+
+ pressure_field.setTime(time,it);
+ pressure_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure",False);
+ velocity_field.setTime(time,it);
+ velocity_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity",False);
+
+ print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
+ print "Variation temporelle relative : pressure ", max_dp/p0 ,", velocity ", max_dq/rho0
+ print
+
+ if(it>=ntmax):
+ print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
+ elif(isStationary):
+ print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
+ print "------------------------------------------------------------------------------------"
+
+ pressure_field.setTime(time,0);
+ pressure_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_Stat");
+ velocity_field.setTime(time,0);
+ velocity_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_Stat");
+
+ #Postprocessing : save 2D picture
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_Stat")
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_Stat")
+
+ else:
+ print "Temps maximum Tmax= ", tmax, " atteint"
+
+
+def solve(my_mesh,meshName,resolution):
+ print "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension()
+ print "Numerical method : staggered scheme"
+ print "Initial data : Spherical wave"
+ print "Periodic boundary conditions"
+ print "Mesh name : ",meshName , my_mesh.getNumberOfCells(), " cells"
+
+ # Problem data
+ tmax = 1000.
+ ntmax = 500
+ cfl = 1./my_mesh.getSpaceDimension()
+ output_freq = 1
+
+ WaveSystemStaggered(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution)
+
+def solve_file( filename,meshName, resolution):
+ my_mesh = cdmath.Mesh(filename+".med")
+
+ return solve(my_mesh, meshName+str(my_mesh.getNumberOfCells()),resolution)
+
+
+if __name__ == """__main__""":
+ if len(sys.argv) >1 :
+ my_mesh = cdmath.Mesh(sys.argv[1])
+ solve(my_mesh,my_mesh.getName(),100)
+ else :
+ nx=50
+ my_mesh = cdmath.Mesh(0,1,nx,0,1,nx)
+ solve(my_mesh,my_mesh.getName(),100)
--- /dev/null
+#!/usr/bin/env python
+# -*-coding:utf-8 -*
+
+from math import sin, cos, pi, sqrt
+import time, json
+import cdmath
+import PV_routines
+import VTK_routines
+
+p0=155e5#reference pressure in a pressurised nuclear vessel
+c0=700.#reference sound speed for water at 155 bars
+rho0=p0/c0*c0#reference density
+precision=1e-5
+
+def initial_conditions_square_shock(my_mesh):
+ print "Spherical wave in a square : initial data"
+ dim = my_mesh.getMeshDimension()
+ nbCells = my_mesh.getNumberOfCells()
+
+ rayon = 0.15
+ xcentre = 0.5
+ ycentre = 0.5
+ zcentre = 0.5
+
+
+ pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
+ velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
+
+ for i in range(nbCells):
+ x = my_mesh.getCell(i).x()
+ y = my_mesh.getCell(i).y()
+
+ velocity_field[i,0] = 0
+ velocity_field[i,1] = 0
+ velocity_field[i,2] = 0
+
+ valX = (x - xcentre) * (x - xcentre)
+ valY = (y - ycentre) * (y - ycentre)
+
+ if(dim==2):
+ val = sqrt(valX + valY)
+ if(dim==3):
+ z = my_mesh.getCell(i).z()
+ valZ = (z - zcentre) * (z - zcentre)
+ val = sqrt(valX + valY + valZ)
+
+ if val < rayon:
+ pressure_field[i] = p0
+ pass
+ else:
+ pressure_field[i] = p0/2
+ pass
+ pass
+
+ return pressure_field, velocity_field
+
+
+def jacobianMatrices(normal, coeff,scaling):
+ dim=normal.size()
+ A=cdmath.Matrix(dim+1,dim+1)
+ absA=cdmath.Matrix(dim+1,dim+1)
+
+ absA[0,0]=c0*coeff
+ for i in range(dim):
+ for j in range(dim):
+ absA[i+1,j+1]=c0*normal[i]*normal[j]*coeff
+ if( scaling==0):
+ A[0,i+1]=c0*c0*normal[i]*coeff
+ A[i+1,0]= normal[i]*coeff
+ else:
+ A[0,i+1]= c0*normal[i]*coeff
+ A[i+1,0]= c0*normal[i]*coeff
+
+ return (A-absA)/2
+
+
+def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling,test_bc):
+ nbCells = my_mesh.getNumberOfCells()
+ dim=my_mesh.getMeshDimension()
+ nbComp=dim+1
+ normal=cdmath.Vector(dim)
+
+ implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
+
+ idMoinsJacCL=cdmath.Matrix(nbComp)
+
+ for j in range(nbCells):#On parcourt les cellules
+ Cj = my_mesh.getCell(j)
+ nbFaces = Cj.getNumberOfFaces();
+
+ for k in range(nbFaces) :
+ indexFace = Cj.getFacesId()[k];
+ Fk = my_mesh.getFace(indexFace);
+ for i in range(dim) :
+ normal[i] = Cj.getNormalVector(k, i);#normale sortante
+
+ Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure(),scaling);
+
+ cellAutre =-1
+ if ( not Fk.isBorder()) :
+ # hypothese: La cellule d'index indexC1 est la cellule courante index j */
+ if (Fk.getCellsId()[0] == j) :
+ # hypothese verifiée
+ cellAutre = Fk.getCellsId()[1];
+ elif(Fk.getCellsId()[1] == j) :
+ # hypothese non verifiée
+ cellAutre = Fk.getCellsId()[0];
+ else :
+ raise ValueError("computeFluxes: problem with mesh, unknown cell number")
+
+ implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
+ implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
+ else :
+ indexFP = my_mesh.getIndexFacePeriodic(indexFace, my_mesh.getName()== "squareWithBrickWall", my_mesh.getName()== "squareWithHexagons")
+ Fp = my_mesh.getFace(indexFP)
+ cellAutre = Fp.getCellsId()[0]
+
+ implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
+ implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
+
+ return implMat
+
+def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution,scaling,test_bc,with_source=False):
+ dim=my_mesh.getMeshDimension()
+ nbCells = my_mesh.getNumberOfCells()
+
+ dt = 0.
+ time = 0.
+ it=0;
+ isStationary=False
+
+ nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
+ isImplicit=False
+
+ #iteration vectors
+ Un =cdmath.Vector(nbCells*(dim+1))
+ dUn=cdmath.Vector(nbCells*(dim+1))
+
+ # Initial conditions #
+ print("Construction of the initial condition …")
+ if(meshName.find("square")>-1 or meshName.find("Square")>-1 or meshName.find("cube")>-1 or meshName.find("Cube")>-1):
+ pressure_field, velocity_field = initial_conditions_square_shock(my_mesh)
+ stat_pressure, stat_velocity = initial_conditions_square_shock(my_mesh)
+ elif(meshName.find("disk")>-1 or meshName.find("Disk")>-1):
+ pressure_field, velocity_field = initial_conditions_disk_vortex(my_mesh)
+ stat_pressure, stat_velocity = initial_conditions_disk_vortex(my_mesh)
+ else:
+ print "Mesh name : ", meshName
+ raise ValueError("Mesh name should contain substring square, cube or disk")
+
+ for k in range(nbCells):
+ Un[k*(dim+1)+0] = pressure_field[k]
+ Un[k*(dim+1)+1] =rho0*velocity_field[k,0]
+ Un[k*(dim+1)+2] =rho0*velocity_field[k,1]
+ if(dim==3):
+ Un[k*(dim+1)+3] =rho0*velocity_field[k,2]
+
+ #sauvegarde de la donnée initiale
+ pressure_field.setTime(time,it);
+ pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure");
+ velocity_field.setTime(time,it);
+ velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity");
+ #Postprocessing : save 2D picture
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_initial")
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_initial")
+
+ dx_min=my_mesh.minRatioVolSurf()
+
+ dt = cfl * dx_min / c0
+
+ divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling,test_bc)
+ if( isImplicit):
+ #Adding the identity matrix on the diagonal
+ divMat.diagonalShift(1)#only after filling all coefficients
+ #for j in range(nbCells*(dim+1)):
+ # divMat.addValue(j,j,1)
+
+ iterGMRESMax=50
+ LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","ILU")
+
+ LS.setComputeConditionNumber()
+
+ print("Starting computation of the linear wave system with an Upwind scheme …")
+
+ # Starting time loop
+ while (it<ntmax and time <= tmax and not isStationary):
+ if(isImplicit):
+ dUn=Un.deepCopy()
+ LS.setSndMember(Un)
+ Un=LS.solve();
+ if(not LS.getStatus()):
+ print "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations"
+ raise ValueError("Pas de convergence du système linéaire");
+ dUn-=Un
+
+ else:
+ dUn=divMat*Un
+ Un-=dUn
+
+ time=time+dt;
+ it=it+1;
+
+ #Sauvegardes
+ if(it==1 or it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
+ print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
+ if(isImplicit):
+ print "Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations"
+
+ for k in range(nbCells):
+ pressure_field[k] =Un[k*(dim+1)+0]
+ velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
+ if(dim>1):
+ velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
+ if(dim>2):
+ velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
+
+ pressure_field.setTime(time,it);
+ pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure",False);
+ velocity_field.setTime(time,it);
+ velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity",False);
+
+ print
+ print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
+
+ if(it>=ntmax):
+ print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
+ raise ValueError("Maximum number of time steps reached : Stationary state not found !!!!!!!")
+ elif(isStationary):
+ print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
+ print "------------------------------------------------------------------------------------"
+
+ pressure_field.setTime(time,0);
+ pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat");
+ velocity_field.setTime(time,0);
+ velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat");
+
+ #Postprocessing : Extraction of the diagonal data
+ if(dim==2):
+ diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,1,0],[1,0,0], resolution)
+ diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,1,0],[1,0,0], resolution)
+ elif(dim==3):
+ diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,0,0],[1,1,1], resolution)
+ diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,0,0],[1,1,1], resolution)
+ #Postprocessing : save 2D picture
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat")
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat")
+
+ else:
+ print "Temps maximum Tmax= ", tmax, " atteint"
+ raise ValueError("Maximum time reached : Stationary state not found !!!!!!!")
+
+
+def solve(my_mesh, meshName, resolution, scaling, meshType, testColor, cfl, test_bc, with_source=False):
+ print "Resolution of the Wave system in dimension " +str( my_mesh.getSpaceDimension())+" on "+str(my_mesh.getNumberOfCells())+ " cells"
+ print "Initial data : ", "Riemann problem"
+ print "Boundary conditions : ", "Neumann"
+ print "Mesh name : ",meshName , ", ", my_mesh.getNumberOfCells(), " cells"
+
+ # Problem data
+ tmax = 10000.
+ ntmax = 500
+ output_freq = 1
+
+ WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution, scaling,test_bc,with_source)
+
+ return
+
+def solve_file( filename,meshName, resolution, scaling, meshType, testColor,cfl,test_bc,with_source=False):
+ my_mesh = cdmath.Mesh(filename+".med")
+
+ return solve(my_mesh, meshName+str(my_mesh.getNumberOfCells()),resolution, scaling, meshType, testColor,cfl,test_bc,with_source)
+
+
+if __name__ == """__main__""":
+ M1=cdmath.Mesh(0.,1.,50,0.,1.,50)
+ cfl=0.5
+ scaling=0
+ solve(M1,"SquareRegularSquares",100,scaling,"RegularSquares","Green",cfl,"Periodic",False)
--- /dev/null
+#!/usr/bin/env python
+# -*-coding:utf-8 -*
+
+from math import sin, cos, pi, sqrt
+import time, json
+import cdmath
+import PV_routines
+import VTK_routines
+
+p0=155e5#reference pressure in a pressurised nuclear vessel
+c0=700.#reference sound speed for water at 155 bars
+rho0=p0/c0*c0#reference density
+precision=1e-5
+
+def initial_conditions_square_shock(my_mesh):
+ print "Spherical wave in a square : initial data"
+ dim = my_mesh.getMeshDimension()
+ nbCells = my_mesh.getNumberOfCells()
+
+ rayon = 0.15
+ xcentre = 0.5
+ ycentre = 0.5
+ zcentre = 0.5
+
+
+ pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
+ velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
+
+ for i in range(nbCells):
+ x = my_mesh.getCell(i).x()
+ y = my_mesh.getCell(i).y()
+
+ velocity_field[i,0] = 0
+ velocity_field[i,1] = 0
+ velocity_field[i,2] = 0
+
+ valX = (x - xcentre) * (x - xcentre)
+ valY = (y - ycentre) * (y - ycentre)
+
+ if(dim==2):
+ val = sqrt(valX + valY)
+ if(dim==3):
+ z = my_mesh.getCell(i).z()
+ valZ = (z - zcentre) * (z - zcentre)
+ val = sqrt(valX + valY + valZ)
+
+ if val < rayon:
+ pressure_field[i] = p0
+ pass
+ else:
+ pressure_field[i] = p0/2
+ pass
+ pass
+
+ return pressure_field, velocity_field
+
+
+def jacobianMatrices(normal, coeff,scaling):
+ dim=normal.size()
+ A=cdmath.Matrix(dim+1,dim+1)
+ absA=cdmath.Matrix(dim+1,dim+1)
+
+ absA[0,0]=c0*coeff
+ for i in range(dim):
+ for j in range(dim):
+ absA[i+1,j+1]=c0*normal[i]*normal[j]*coeff
+ if( scaling==0):
+ A[0,i+1]=c0*c0*normal[i]*coeff
+ A[i+1,0]= normal[i]*coeff
+ else:
+ A[0,i+1]= c0*normal[i]*coeff
+ A[i+1,0]= c0*normal[i]*coeff
+
+ return (A-absA)/2
+
+
+def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling,test_bc):
+ nbCells = my_mesh.getNumberOfCells()
+ dim=my_mesh.getMeshDimension()
+ nbComp=dim+1
+ normal=cdmath.Vector(dim)
+
+ implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
+
+ idMoinsJacCL=cdmath.Matrix(nbComp)
+
+ for j in range(nbCells):#On parcourt les cellules
+ Cj = my_mesh.getCell(j)
+ nbFaces = Cj.getNumberOfFaces();
+
+ for k in range(nbFaces) :
+ indexFace = Cj.getFacesId()[k];
+ Fk = my_mesh.getFace(indexFace);
+ for i in range(dim) :
+ normal[i] = Cj.getNormalVector(k, i);#normale sortante
+
+ Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure(),scaling);
+
+ cellAutre =-1
+ if ( not Fk.isBorder()) :
+ # hypothese: La cellule d'index indexC1 est la cellule courante index j */
+ if (Fk.getCellsId()[0] == j) :
+ # hypothese verifiée
+ cellAutre = Fk.getCellsId()[1];
+ elif(Fk.getCellsId()[1] == j) :
+ # hypothese non verifiée
+ cellAutre = Fk.getCellsId()[0];
+ else :
+ raise ValueError("computeFluxes: problem with mesh, unknown cell number")
+
+ implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
+ implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
+ else :
+ indexFP = my_mesh.getIndexFacePeriodic(indexFace, my_mesh.getName()== "squareWithBrickWall", my_mesh.getName()== "squareWithHexagons")
+ Fp = my_mesh.getFace(indexFP)
+ cellAutre = Fp.getCellsId()[0]
+
+ implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
+ implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
+
+ return implMat
+
+def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution,scaling,test_bc,with_source=False):
+ dim=my_mesh.getMeshDimension()
+ nbCells = my_mesh.getNumberOfCells()
+
+ dt = 0.
+ time = 0.
+ it=0;
+ isStationary=False
+
+ nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
+ isImplicit=True
+
+ #iteration vectors
+ Un =cdmath.Vector(nbCells*(dim+1))
+ dUn=cdmath.Vector(nbCells*(dim+1))
+
+ # Initial conditions #
+ print("Construction of the initial condition …")
+ if(meshName.find("square")>-1 or meshName.find("Square")>-1 or meshName.find("cube")>-1 or meshName.find("Cube")>-1):
+ pressure_field, velocity_field = initial_conditions_square_shock(my_mesh)
+ stat_pressure, stat_velocity = initial_conditions_square_shock(my_mesh)
+ elif(meshName.find("disk")>-1 or meshName.find("Disk")>-1):
+ pressure_field, velocity_field = initial_conditions_disk_vortex(my_mesh)
+ stat_pressure, stat_velocity = initial_conditions_disk_vortex(my_mesh)
+ else:
+ print "Mesh name : ", meshName
+ raise ValueError("Mesh name should contain substring square, cube or disk")
+
+ for k in range(nbCells):
+ Un[k*(dim+1)+0] = pressure_field[k]
+ Un[k*(dim+1)+1] =rho0*velocity_field[k,0]
+ Un[k*(dim+1)+2] =rho0*velocity_field[k,1]
+ if(dim==3):
+ Un[k*(dim+1)+3] =rho0*velocity_field[k,2]
+
+ #sauvegarde de la donnée initiale
+ pressure_field.setTime(time,it);
+ pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure");
+ velocity_field.setTime(time,it);
+ velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity");
+ #Postprocessing : save 2D picture
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_initial")
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_initial")
+
+ dx_min=my_mesh.minRatioVolSurf()
+
+ dt = cfl * dx_min / c0
+
+ divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling,test_bc)
+ if( isImplicit):
+ #Adding the identity matrix on the diagonal
+ divMat.diagonalShift(1)#only after filling all coefficients
+ #for j in range(nbCells*(dim+1)):
+ # divMat.addValue(j,j,1)
+
+ iterGMRESMax=50
+ LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","ILU")
+
+ LS.setComputeConditionNumber()
+
+ print("Starting computation of the linear wave system with an Upwind scheme …")
+
+ # Starting time loop
+ while (it<ntmax and time <= tmax and not isStationary):
+ if(isImplicit):
+ dUn=Un.deepCopy()
+ LS.setSndMember(Un)
+ Un=LS.solve();
+ if(not LS.getStatus()):
+ print "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations"
+ raise ValueError("Pas de convergence du système linéaire");
+ dUn-=Un
+
+ else:
+ dUn=divMat*Un-S*dt
+ Un-=dUn
+
+ time=time+dt;
+ it=it+1;
+
+ #Sauvegardes
+ if(it==1 or it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
+ print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
+ if(isImplicit):
+ print "Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations"
+
+ for k in range(nbCells):
+ pressure_field[k] =Un[k*(dim+1)+0]
+ velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
+ if(dim>1):
+ velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
+ if(dim>2):
+ velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
+
+ pressure_field.setTime(time,it);
+ pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure",False);
+ velocity_field.setTime(time,it);
+ velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity",False);
+
+ print
+ print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
+
+ if(it>=ntmax):
+ print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
+ raise ValueError("Maximum number of time steps reached : Stationary state not found !!!!!!!")
+ elif(isStationary):
+ print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
+ print "------------------------------------------------------------------------------------"
+
+ pressure_field.setTime(time,0);
+ pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat");
+ velocity_field.setTime(time,0);
+ velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat");
+
+ #Postprocessing : Extraction of the diagonal data
+ if(dim==2):
+ diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,1,0],[1,0,0], resolution)
+ diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,1,0],[1,0,0], resolution)
+ elif(dim==3):
+ diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,0,0],[1,1,1], resolution)
+ diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,0,0],[1,1,1], resolution)
+ #Postprocessing : save 2D picture
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat")
+ PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat")
+
+ else:
+ print "Temps maximum Tmax= ", tmax, " atteint"
+ raise ValueError("Maximum time reached : Stationary state not found !!!!!!!")
+
+
+def solve(my_mesh, meshName, resolution, scaling, meshType, testColor, cfl, test_bc, with_source=False):
+ print "Resolution of the Wave system in dimension " +str( my_mesh.getSpaceDimension())+" on "+str(my_mesh.getNumberOfCells())+ " cells"
+ print "Initial data : ", "Riemann problem"
+ print "Boundary conditions : ", "Neumann"
+ print "Mesh name : ",meshName , ", ", my_mesh.getNumberOfCells(), " cells"
+
+ # Problem data
+ tmax = 10000.
+ ntmax = 500
+ output_freq = 1
+
+ WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution, scaling,test_bc,with_source)
+
+ return
+
+def solve_file( filename,meshName, resolution, scaling, meshType, testColor,cfl,test_bc,with_source=False):
+ my_mesh = cdmath.Mesh(filename+".med")
+
+ return solve(my_mesh, meshName+str(my_mesh.getNumberOfCells()),resolution, scaling, meshType, testColor,cfl,test_bc,with_source)
+
+
+if __name__ == """__main__""":
+ M1=cdmath.Mesh(0.,1.,50,0.,1.,50)
+ cfl=0.5
+ scaling=0
+ solve(M1,"SquareRegularSquares",100,scaling,"RegularSquares","Green",cfl,"Periodic",False)
+++ /dev/null
-#!/usr/bin/env python
-# -*-coding:utf-8 -*
-
-from math import sin, cos, pi, sqrt
-import cdmath
-import PV_routines
-import VTK_routines
-import sys
-
-c0=700.#reference sound speed
-p0=155e5#reference pressure
-rho0=p0/c0*c0#reference density
-precision=1e-5
-
-def initial_conditions_wave_system(my_mesh):
- print "Spherical wave initial data"
- dim = my_mesh.getMeshDimension()
- nbCells = my_mesh.getNumberOfCells()
-
- rayon = 0.15
- xcentre = 0.5
- ycentre = 0.5
- zcentre = 0.5
-
- pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
- velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
-
- for i in range(nbCells):
- x = my_mesh.getCell(i).x()
- y = my_mesh.getCell(i).y()
-
- velocity_field[i,0] = 0
- velocity_field[i,1] = 0
- velocity_field[i,2] = 0
-
- valX = (x - xcentre) * (x - xcentre)
- valY = (y - ycentre) * (y - ycentre)
- if(dim==2):
- val = sqrt(valX + valY)
- if(dim==3):
- z = my_mesh.getCell(i).z()
- valZ = (z - zcentre) * (z - zcentre)
- val = sqrt(valX + valY + valZ)
-
- if val < rayon:
- pressure_field[i] = p0
- pass
- else:
- pressure_field[i] = p0/2
- pass
- pass
-
- return pressure_field, velocity_field
-
-def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling):
- nbCells = my_mesh.getNumberOfCells()
- dim=my_mesh.getMeshDimension()
- nbComp=dim+1
-
- if(not my_mesh.isStructured()):
- raise ValueError("WaveSystemStaggered: the mesh should be structured");
-
- NxNyNz=my_mesh.getCellGridStructure()
- DxDyDz=my_mesh.getDXYZ()
-
- implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
-
- idMoinsJacCL=cdmath.Matrix(nbComp)
-
- if( dim == 1) :
- nx=NxNyNz[0]
- dx=DxDyDz[0]
-
- if( scaling==0 ):
- for k in range(nbCells):
- implMat.addValue(k,1*nbCells + k , -c0*c0*dt/dx)
- implMat.addValue(k,1*nbCells + (k+1)%nx, c0*c0*dt/dx)
-
- implMat.addValue( 1*nbCells + k ,k, dt/dx)
- implMat.addValue( 1*nbCells + (k+1)%nx,k, -dt/dx)
- else : # scaling >0
- for k in range(nbCells):
- implMat.addValue(k,1*nbCells + k , -c0*dt/dx)
- implMat.addValue(k,1*nbCells + (k+1)%nx, c0*dt/dx)
-
- implMat.addValue( 1*nbCells + k ,k, c0*dt/dx)
- implMat.addValue( 1*nbCells + (k+1)%nx,k, -c0*dt/dx)
-
- elif( dim == 2) :# k = j*nx+i
- nx=NxNyNz[0]
- ny=NxNyNz[1]
- dx=DxDyDz[0]
- dy=DxDyDz[1]
-
- if( scaling==0 ):
- for k in range(nbCells):
- i = k % nx
- j = k //nx
-
- implMat.addValue(k,1*nbCells + j*nx + i , -c0*c0*dt/dx)
- implMat.addValue(k,1*nbCells + j*nx + (i+1)%nx, c0*c0*dt/dx)
-
- implMat.addValue(k,2*nbCells + j *nx + i, -c0*c0*dt/dy)
- implMat.addValue(k,2*nbCells + ((j+1)%ny)*nx + i, c0*c0*dt/dy)
-
- implMat.addValue( 1*nbCells + j*nx + i , k, dt/dx)
- implMat.addValue( 1*nbCells + j*nx + (i+1)%nx, k, -dt/dx)
-
- implMat.addValue( 2*nbCells + j *nx + i,k, dt/dy)
- implMat.addValue( 2*nbCells + ((j+1)%ny)*nx + i,k, -dt/dy)
-
- else :# scaling >0
- for k in range(nbCells):
- i = k % nx
- j = k //nx
-
- implMat.addValue(k,1*nbCells + j*nx + i , -c0*dt/dx)
- implMat.addValue(k,1*nbCells + j*nx + (i+1)%nx, c0*dt/dx)
-
- implMat.addValue(k,2*nbCells + j *nx + i, -c0*dt/dy)
- implMat.addValue(k,2*nbCells + ((j+1)%ny)*nx + i, c0*dt/dy)
-
- implMat.addValue( 1*nbCells + j*nx + i , k, c0*dt/dx)
- implMat.addValue( 1*nbCells + j*nx + (i+1)%nx, k, -c0*dt/dx)
-
- implMat.addValue( 2*nbCells + j *nx + i,k, c0*dt/dy)
- implMat.addValue( 2*nbCells + ((j+1)%ny)*nx + i,k, -c0*dt/dy)
-
- elif( dim == 3) :# k = l*nx*ny+j*nx+i
- nx=NxNyNz[0]
- ny=NxNyNz[1]
- nz=NxNyNz[2]
- dx=DxDyDz[0]
- dy=DxDyDz[1]
- dz=DxDyDz[2]
-
- if( scaling==0 ):
- for k in range(nbCells):
- i = k % nx
- j = (k //nx)%ny
- l = k //(nx*ny)
-
- implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + i , -c0*c0*dt/dx)
- implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + (i+1)%nx, c0*c0*dt/dx)
-
- implMat.addValue(k,2*nbCells + l*nx*ny + j *nx + i, -c0*c0*dt/dy)
- implMat.addValue(k,2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i, c0*c0*dt/dy)
-
- implMat.addValue(k,3*nbCells + l*nx*ny + j*nx + i, -c0*c0*dt/dz)
- implMat.addValue(k,3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i, c0*c0*dt/dz)
-
- implMat.addValue( 1*nbCells + l*nx*ny + j*nx + i , k, dt/dx)
- implMat.addValue( 1*nbCells + l*nx*ny + j*nx + (i+1)%nx, k, -dt/dx)
-
- implMat.addValue( 2*nbCells + l*nx*ny + j *nx + i,k, dt/dy)
- implMat.addValue( 2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i,k, -dt/dy)
-
- implMat.addValue( 3*nbCells + l*nx*ny + j*nx + i,k, dt/dz)
- implMat.addValue( 3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i,k, -dt/dz)
-
- else:# scaling >0
- for k in range(nbCells):
- i = k % nx
- j = (k //nx)%ny
- l = k //(nx*ny)
-
- implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + i , -c0*dt/dx)
- implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + (i+1)%nx, c0*dt/dx)
-
- implMat.addValue(k,2*nbCells + l*nx*ny + j *nx + i, -c0*dt/dy)
- implMat.addValue(k,2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i, c0*dt/dy)
-
- implMat.addValue(k,3*nbCells + l*nx*ny + j*nx + i, -c0*dt/dz)
- implMat.addValue(k,3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i, c0*dt/dz)
-
- implMat.addValue( 1*nbCells + l*nx*ny + j*nx + i , k, c0*dt/dx)
- implMat.addValue( 1*nbCells + l*nx*ny + j*nx + (i+1)%nx, k, -c0*dt/dx)
-
- implMat.addValue( 2*nbCells + l*nx*ny + j *nx + i,k, c0*dt/dy)
- implMat.addValue( 2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i,k, -c0*dt/dy)
-
- implMat.addValue( 3*nbCells + l*nx*ny + j*nx + i,k, c0*dt/dz)
- implMat.addValue( 3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i,k, -c0*dt/dz)
-
- return implMat
-
-def WaveSystemStaggered(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution):
- dim=my_mesh.getMeshDimension()
- nbCells = my_mesh.getNumberOfCells()
-
- dt = 0.
- time = 0.
- it=0;
- isStationary=False;
-
- scaling=0
-
- nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
- iterGMRESMax=50
-
- #iteration vectors
- Un=cdmath.Vector(nbCells*(dim+1))
- dUn=cdmath.Vector(nbCells*(dim+1))
-
- # Initial conditions #
- print("Construction of the initial condition …")
- pressure_field, velocity_field = initial_conditions_wave_system(my_mesh)
-
- for k in range(nbCells):
- Un[k + 0*nbCells] = pressure_field[k]
- Un[k + 1*nbCells] = rho0*velocity_field[k,0] # value on the left face
- Un[k + 2*nbCells] = rho0*velocity_field[k,1] # value on the bottom face
- if(dim==3):
- Un[k + 3*nbCells] = rho0*initial_velocity[k,2]
- if( scaling>0):
- Vn = Un.deepCopy()
- for k in range(nbCells):
- Vn[k] = Vn[k]/c0
-
- #sauvegarde de la donnée initiale
- pressure_field.setTime(time,it);
- pressure_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure");
- velocity_field.setTime(time,it);
- velocity_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity");
- #Postprocessing : save 2D picture
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_initial")
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_initial")
-
- dx_min=my_mesh.minRatioVolSurf()
-
- dt = cfl * dx_min / c0
-
- divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling)
-
- #Add the identity matrix on the diagonal
- for j in range(nbCells*(dim+1)):
- divMat.addValue(j,j,1)
-
- LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","LU")
-
- print("Starting computation of the linear wave system with an pseudo staggered scheme …")
-
- # Starting time loop
- while (it<ntmax and time <= tmax and not isStationary):
- dUn=Un.deepCopy()
- LS.setSndMember(Un)
- Un=LS.solve();
- cvgceLS=LS.getStatus();
- iterGMRES=LS.getNumberOfIter();
- if(not cvgceLS):
- print "Linear system did not converge ", iterGMRES, " GMRES iterations"
- raise ValueError("Pas de convergence du système linéaire");
- dUn-=Un
-
- max_dp=0 ; max_dq=0
- for k in range(nbCells):
- max_dp = max(max_dp,abs(dUn[k]))
- for i in range(dim):
- max_dq=max(max_dq,abs(dUn[k+(1+i)*nbCells]))
-
- isStationary= max_dp/p0<precision and max_dq/rho0<precision
-
- time=time+dt;
- it=it+1;
-
- #Sauvegardes
- if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
- print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
- print "Variation temporelle relative : pressure ", max_dp/p0 ,", velocity ", max_dq/rho0
- print "Linear system converged in ", iterGMRES, " GMRES iterations"
-
- for k in range(nbCells):
- pressure_field[k]=Un[k]
- velocity_field[k,0]=Un[k+1*nbCells]/rho0
- if(dim>1):
- velocity_field[k,1]=Un[k+2*nbCells]/rho0
- if(dim>2):
- velocity_field[k,2]=Un[k+3*nbCells]/rho0
-
- pressure_field.setTime(time,it);
- pressure_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure",False);
- velocity_field.setTime(time,it);
- velocity_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity",False);
-
- print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
- print "Variation temporelle relative : pressure ", max_dp/p0 ,", velocity ", max_dq/rho0
- print
-
- if(it>=ntmax):
- print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
- elif(isStationary):
- print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
- print "------------------------------------------------------------------------------------"
-
- pressure_field.setTime(time,0);
- pressure_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_Stat");
- velocity_field.setTime(time,0);
- velocity_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_Stat");
-
- #Postprocessing : save 2D picture
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_Stat")
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_Stat")
-
- else:
- print "Temps maximum Tmax= ", tmax, " atteint"
-
-
-def solve(my_mesh,meshName,resolution):
- print "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension()
- print "Numerical method : staggered scheme"
- print "Initial data : Spherical wave"
- print "Periodic boundary conditions"
- print "Mesh name : ",meshName , my_mesh.getNumberOfCells(), " cells"
-
- # Problem data
- tmax = 1000.
- ntmax = 500
- cfl = 1./my_mesh.getSpaceDimension()
- output_freq = 1
-
- WaveSystemStaggered(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution)
-
-def solve_file( filename,meshName, resolution):
- my_mesh = cdmath.Mesh(filename+".med")
-
- return solve(my_mesh, meshName+str(my_mesh.getNumberOfCells()),resolution)
-
-
-if __name__ == """__main__""":
- if len(sys.argv) >1 :
- my_mesh = cdmath.Mesh(sys.argv[1])
- solve(my_mesh,my_mesh.getName(),100)
- else :
- nx=50
- my_mesh = cdmath.Mesh(0,1,nx,0,1,nx)
- solve(my_mesh,my_mesh.getName(),100)
+++ /dev/null
-#!/usr/bin/env python
-# -*-coding:utf-8 -*
-
-from math import sin, cos, pi, sqrt
-import time, json
-import cdmath
-import PV_routines
-import VTK_routines
-
-p0=155e5#reference pressure in a pressurised nuclear vessel
-c0=700.#reference sound speed for water at 155 bars
-rho0=p0/c0*c0#reference density
-precision=1e-5
-
-def initial_conditions_square_shock(my_mesh):
- print "Spherical wave in a square : initial data"
- dim = my_mesh.getMeshDimension()
- nbCells = my_mesh.getNumberOfCells()
-
- rayon = 0.15
- xcentre = 0.5
- ycentre = 0.5
- zcentre = 0.5
-
-
- pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
- velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
-
- for i in range(nbCells):
- x = my_mesh.getCell(i).x()
- y = my_mesh.getCell(i).y()
-
- velocity_field[i,0] = 0
- velocity_field[i,1] = 0
- velocity_field[i,2] = 0
-
- valX = (x - xcentre) * (x - xcentre)
- valY = (y - ycentre) * (y - ycentre)
-
- if(dim==2):
- val = sqrt(valX + valY)
- if(dim==3):
- z = my_mesh.getCell(i).z()
- valZ = (z - zcentre) * (z - zcentre)
- val = sqrt(valX + valY + valZ)
-
- if val < rayon:
- pressure_field[i] = p0
- pass
- else:
- pressure_field[i] = p0/2
- pass
- pass
-
- return pressure_field, velocity_field
-
-
-def jacobianMatrices(normal, coeff,scaling):
- dim=normal.size()
- A=cdmath.Matrix(dim+1,dim+1)
- absA=cdmath.Matrix(dim+1,dim+1)
-
- absA[0,0]=c0*coeff
- for i in range(dim):
- for j in range(dim):
- absA[i+1,j+1]=c0*normal[i]*normal[j]*coeff
- if( scaling==0):
- A[0,i+1]=c0*c0*normal[i]*coeff
- A[i+1,0]= normal[i]*coeff
- else:
- A[0,i+1]= c0*normal[i]*coeff
- A[i+1,0]= c0*normal[i]*coeff
-
- return (A-absA)/2
-
-
-def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling,test_bc):
- nbCells = my_mesh.getNumberOfCells()
- dim=my_mesh.getMeshDimension()
- nbComp=dim+1
- normal=cdmath.Vector(dim)
-
- implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
-
- idMoinsJacCL=cdmath.Matrix(nbComp)
-
- for j in range(nbCells):#On parcourt les cellules
- Cj = my_mesh.getCell(j)
- nbFaces = Cj.getNumberOfFaces();
-
- for k in range(nbFaces) :
- indexFace = Cj.getFacesId()[k];
- Fk = my_mesh.getFace(indexFace);
- for i in range(dim) :
- normal[i] = Cj.getNormalVector(k, i);#normale sortante
-
- Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure(),scaling);
-
- cellAutre =-1
- if ( not Fk.isBorder()) :
- # hypothese: La cellule d'index indexC1 est la cellule courante index j */
- if (Fk.getCellsId()[0] == j) :
- # hypothese verifiée
- cellAutre = Fk.getCellsId()[1];
- elif(Fk.getCellsId()[1] == j) :
- # hypothese non verifiée
- cellAutre = Fk.getCellsId()[0];
- else :
- raise ValueError("computeFluxes: problem with mesh, unknown cell number")
-
- implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
- implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
- else :
- indexFP = my_mesh.getIndexFacePeriodic(indexFace, my_mesh.getName()== "squareWithBrickWall", my_mesh.getName()== "squareWithHexagons")
- Fp = my_mesh.getFace(indexFP)
- cellAutre = Fp.getCellsId()[0]
-
- implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
- implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
-
- return implMat
-
-def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution,scaling,test_bc,with_source=False):
- dim=my_mesh.getMeshDimension()
- nbCells = my_mesh.getNumberOfCells()
-
- dt = 0.
- time = 0.
- it=0;
- isStationary=False
-
- nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
- isImplicit=False
-
- #iteration vectors
- Un =cdmath.Vector(nbCells*(dim+1))
- dUn=cdmath.Vector(nbCells*(dim+1))
-
- # Initial conditions #
- print("Construction of the initial condition …")
- if(meshName.find("square")>-1 or meshName.find("Square")>-1 or meshName.find("cube")>-1 or meshName.find("Cube")>-1):
- pressure_field, velocity_field = initial_conditions_square_shock(my_mesh)
- stat_pressure, stat_velocity = initial_conditions_square_shock(my_mesh)
- elif(meshName.find("disk")>-1 or meshName.find("Disk")>-1):
- pressure_field, velocity_field = initial_conditions_disk_vortex(my_mesh)
- stat_pressure, stat_velocity = initial_conditions_disk_vortex(my_mesh)
- else:
- print "Mesh name : ", meshName
- raise ValueError("Mesh name should contain substring square, cube or disk")
-
- for k in range(nbCells):
- Un[k*(dim+1)+0] = pressure_field[k]
- Un[k*(dim+1)+1] =rho0*velocity_field[k,0]
- Un[k*(dim+1)+2] =rho0*velocity_field[k,1]
- if(dim==3):
- Un[k*(dim+1)+3] =rho0*velocity_field[k,2]
-
- #sauvegarde de la donnée initiale
- pressure_field.setTime(time,it);
- pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure");
- velocity_field.setTime(time,it);
- velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity");
- #Postprocessing : save 2D picture
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_initial")
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_initial")
-
- dx_min=my_mesh.minRatioVolSurf()
-
- dt = cfl * dx_min / c0
-
- divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling,test_bc)
- if( isImplicit):
- #Adding the identity matrix on the diagonal
- divMat.diagonalShift(1)#only after filling all coefficients
- #for j in range(nbCells*(dim+1)):
- # divMat.addValue(j,j,1)
-
- iterGMRESMax=50
- LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","ILU")
-
- LS.setComputeConditionNumber()
-
- print("Starting computation of the linear wave system with an Upwind scheme …")
-
- # Starting time loop
- while (it<ntmax and time <= tmax and not isStationary):
- if(isImplicit):
- dUn=Un.deepCopy()
- LS.setSndMember(Un)
- Un=LS.solve();
- if(not LS.getStatus()):
- print "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations"
- raise ValueError("Pas de convergence du système linéaire");
- dUn-=Un
-
- else:
- dUn=divMat*Un
- Un-=dUn
-
- time=time+dt;
- it=it+1;
-
- #Sauvegardes
- if(it==1 or it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
- print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
- if(isImplicit):
- print "Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations"
-
- for k in range(nbCells):
- pressure_field[k] =Un[k*(dim+1)+0]
- velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
- if(dim>1):
- velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
- if(dim>2):
- velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
-
- pressure_field.setTime(time,it);
- pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure",False);
- velocity_field.setTime(time,it);
- velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity",False);
-
- print
- print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
-
- if(it>=ntmax):
- print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
- raise ValueError("Maximum number of time steps reached : Stationary state not found !!!!!!!")
- elif(isStationary):
- print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
- print "------------------------------------------------------------------------------------"
-
- pressure_field.setTime(time,0);
- pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat");
- velocity_field.setTime(time,0);
- velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat");
-
- #Postprocessing : Extraction of the diagonal data
- if(dim==2):
- diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,1,0],[1,0,0], resolution)
- diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,1,0],[1,0,0], resolution)
- elif(dim==3):
- diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,0,0],[1,1,1], resolution)
- diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,0,0],[1,1,1], resolution)
- #Postprocessing : save 2D picture
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat")
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat")
-
- else:
- print "Temps maximum Tmax= ", tmax, " atteint"
- raise ValueError("Maximum time reached : Stationary state not found !!!!!!!")
-
-
-def solve(my_mesh, meshName, resolution, scaling, meshType, testColor, cfl, test_bc, with_source=False):
- print "Resolution of the Wave system in dimension " +str( my_mesh.getSpaceDimension())+" on "+str(my_mesh.getNumberOfCells())+ " cells"
- print "Initial data : ", "Riemann problem"
- print "Boundary conditions : ", "Neumann"
- print "Mesh name : ",meshName , ", ", my_mesh.getNumberOfCells(), " cells"
-
- # Problem data
- tmax = 10000.
- ntmax = 500
- output_freq = 1
-
- WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution, scaling,test_bc,with_source)
-
- return
-
-def solve_file( filename,meshName, resolution, scaling, meshType, testColor,cfl,test_bc,with_source=False):
- my_mesh = cdmath.Mesh(filename+".med")
-
- return solve(my_mesh, meshName+str(my_mesh.getNumberOfCells()),resolution, scaling, meshType, testColor,cfl,test_bc,with_source)
-
-
-if __name__ == """__main__""":
- M1=cdmath.Mesh(0.,1.,50,0.,1.,50)
- cfl=0.5
- scaling=0
- solve(M1,"SquareRegularSquares",100,scaling,"RegularSquares","Green",cfl,"Periodic",False)
+++ /dev/null
-#!/usr/bin/env python
-# -*-coding:utf-8 -*
-
-from math import sin, cos, pi, sqrt
-import time, json
-import cdmath
-import PV_routines
-import VTK_routines
-
-p0=155e5#reference pressure in a pressurised nuclear vessel
-c0=700.#reference sound speed for water at 155 bars
-rho0=p0/c0*c0#reference density
-precision=1e-5
-
-def initial_conditions_square_shock(my_mesh):
- print "Spherical wave in a square : initial data"
- dim = my_mesh.getMeshDimension()
- nbCells = my_mesh.getNumberOfCells()
-
- rayon = 0.15
- xcentre = 0.5
- ycentre = 0.5
- zcentre = 0.5
-
-
- pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
- velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
-
- for i in range(nbCells):
- x = my_mesh.getCell(i).x()
- y = my_mesh.getCell(i).y()
-
- velocity_field[i,0] = 0
- velocity_field[i,1] = 0
- velocity_field[i,2] = 0
-
- valX = (x - xcentre) * (x - xcentre)
- valY = (y - ycentre) * (y - ycentre)
-
- if(dim==2):
- val = sqrt(valX + valY)
- if(dim==3):
- z = my_mesh.getCell(i).z()
- valZ = (z - zcentre) * (z - zcentre)
- val = sqrt(valX + valY + valZ)
-
- if val < rayon:
- pressure_field[i] = p0
- pass
- else:
- pressure_field[i] = p0/2
- pass
- pass
-
- return pressure_field, velocity_field
-
-
-def jacobianMatrices(normal, coeff,scaling):
- dim=normal.size()
- A=cdmath.Matrix(dim+1,dim+1)
- absA=cdmath.Matrix(dim+1,dim+1)
-
- absA[0,0]=c0*coeff
- for i in range(dim):
- for j in range(dim):
- absA[i+1,j+1]=c0*normal[i]*normal[j]*coeff
- if( scaling==0):
- A[0,i+1]=c0*c0*normal[i]*coeff
- A[i+1,0]= normal[i]*coeff
- else:
- A[0,i+1]= c0*normal[i]*coeff
- A[i+1,0]= c0*normal[i]*coeff
-
- return (A-absA)/2
-
-
-def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling,test_bc):
- nbCells = my_mesh.getNumberOfCells()
- dim=my_mesh.getMeshDimension()
- nbComp=dim+1
- normal=cdmath.Vector(dim)
-
- implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
-
- idMoinsJacCL=cdmath.Matrix(nbComp)
-
- for j in range(nbCells):#On parcourt les cellules
- Cj = my_mesh.getCell(j)
- nbFaces = Cj.getNumberOfFaces();
-
- for k in range(nbFaces) :
- indexFace = Cj.getFacesId()[k];
- Fk = my_mesh.getFace(indexFace);
- for i in range(dim) :
- normal[i] = Cj.getNormalVector(k, i);#normale sortante
-
- Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure(),scaling);
-
- cellAutre =-1
- if ( not Fk.isBorder()) :
- # hypothese: La cellule d'index indexC1 est la cellule courante index j */
- if (Fk.getCellsId()[0] == j) :
- # hypothese verifiée
- cellAutre = Fk.getCellsId()[1];
- elif(Fk.getCellsId()[1] == j) :
- # hypothese non verifiée
- cellAutre = Fk.getCellsId()[0];
- else :
- raise ValueError("computeFluxes: problem with mesh, unknown cell number")
-
- implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
- implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
- else :
- indexFP = my_mesh.getIndexFacePeriodic(indexFace, my_mesh.getName()== "squareWithBrickWall", my_mesh.getName()== "squareWithHexagons")
- Fp = my_mesh.getFace(indexFP)
- cellAutre = Fp.getCellsId()[0]
-
- implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
- implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
-
- return implMat
-
-def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution,scaling,test_bc,with_source=False):
- dim=my_mesh.getMeshDimension()
- nbCells = my_mesh.getNumberOfCells()
-
- dt = 0.
- time = 0.
- it=0;
- isStationary=False
-
- nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
- isImplicit=True
-
- #iteration vectors
- Un =cdmath.Vector(nbCells*(dim+1))
- dUn=cdmath.Vector(nbCells*(dim+1))
-
- # Initial conditions #
- print("Construction of the initial condition …")
- if(meshName.find("square")>-1 or meshName.find("Square")>-1 or meshName.find("cube")>-1 or meshName.find("Cube")>-1):
- pressure_field, velocity_field = initial_conditions_square_shock(my_mesh)
- stat_pressure, stat_velocity = initial_conditions_square_shock(my_mesh)
- elif(meshName.find("disk")>-1 or meshName.find("Disk")>-1):
- pressure_field, velocity_field = initial_conditions_disk_vortex(my_mesh)
- stat_pressure, stat_velocity = initial_conditions_disk_vortex(my_mesh)
- else:
- print "Mesh name : ", meshName
- raise ValueError("Mesh name should contain substring square, cube or disk")
-
- for k in range(nbCells):
- Un[k*(dim+1)+0] = pressure_field[k]
- Un[k*(dim+1)+1] =rho0*velocity_field[k,0]
- Un[k*(dim+1)+2] =rho0*velocity_field[k,1]
- if(dim==3):
- Un[k*(dim+1)+3] =rho0*velocity_field[k,2]
-
- #sauvegarde de la donnée initiale
- pressure_field.setTime(time,it);
- pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure");
- velocity_field.setTime(time,it);
- velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity");
- #Postprocessing : save 2D picture
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_initial")
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_initial")
-
- dx_min=my_mesh.minRatioVolSurf()
-
- dt = cfl * dx_min / c0
-
- divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling,test_bc)
- if( isImplicit):
- #Adding the identity matrix on the diagonal
- divMat.diagonalShift(1)#only after filling all coefficients
- #for j in range(nbCells*(dim+1)):
- # divMat.addValue(j,j,1)
-
- iterGMRESMax=50
- LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","ILU")
-
- LS.setComputeConditionNumber()
-
- print("Starting computation of the linear wave system with an Upwind scheme …")
-
- # Starting time loop
- while (it<ntmax and time <= tmax and not isStationary):
- if(isImplicit):
- dUn=Un.deepCopy()
- LS.setSndMember(Un)
- Un=LS.solve();
- if(not LS.getStatus()):
- print "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations"
- raise ValueError("Pas de convergence du système linéaire");
- dUn-=Un
-
- else:
- dUn=divMat*Un-S*dt
- Un-=dUn
-
- time=time+dt;
- it=it+1;
-
- #Sauvegardes
- if(it==1 or it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
- print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
- if(isImplicit):
- print "Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations"
-
- for k in range(nbCells):
- pressure_field[k] =Un[k*(dim+1)+0]
- velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
- if(dim>1):
- velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
- if(dim>2):
- velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
-
- pressure_field.setTime(time,it);
- pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure",False);
- velocity_field.setTime(time,it);
- velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity",False);
-
- print
- print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
-
- if(it>=ntmax):
- print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
- raise ValueError("Maximum number of time steps reached : Stationary state not found !!!!!!!")
- elif(isStationary):
- print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
- print "------------------------------------------------------------------------------------"
-
- pressure_field.setTime(time,0);
- pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat");
- velocity_field.setTime(time,0);
- velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat");
-
- #Postprocessing : Extraction of the diagonal data
- if(dim==2):
- diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,1,0],[1,0,0], resolution)
- diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,1,0],[1,0,0], resolution)
- elif(dim==3):
- diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,0,0],[1,1,1], resolution)
- diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,0,0],[1,1,1], resolution)
- #Postprocessing : save 2D picture
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat")
- PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat")
-
- else:
- print "Temps maximum Tmax= ", tmax, " atteint"
- raise ValueError("Maximum time reached : Stationary state not found !!!!!!!")
-
-
-def solve(my_mesh, meshName, resolution, scaling, meshType, testColor, cfl, test_bc, with_source=False):
- print "Resolution of the Wave system in dimension " +str( my_mesh.getSpaceDimension())+" on "+str(my_mesh.getNumberOfCells())+ " cells"
- print "Initial data : ", "Riemann problem"
- print "Boundary conditions : ", "Neumann"
- print "Mesh name : ",meshName , ", ", my_mesh.getNumberOfCells(), " cells"
-
- # Problem data
- tmax = 10000.
- ntmax = 500
- output_freq = 1
-
- WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution, scaling,test_bc,with_source)
-
- return
-
-def solve_file( filename,meshName, resolution, scaling, meshType, testColor,cfl,test_bc,with_source=False):
- my_mesh = cdmath.Mesh(filename+".med")
-
- return solve(my_mesh, meshName+str(my_mesh.getNumberOfCells()),resolution, scaling, meshType, testColor,cfl,test_bc,with_source)
-
-
-if __name__ == """__main__""":
- M1=cdmath.Mesh(0.,1.,50,0.,1.,50)
- cfl=0.5
- scaling=0
- solve(M1,"SquareRegularSquares",100,scaling,"RegularSquares","Green",cfl,"Periodic",False)
--- /dev/null
+#!/usr/bin/env python
+# -*-coding:utf-8 -*
+
+import CoreFlows as cf
+import cdmath
+
+
+def SinglePhase_2DSphericalExplosion_HEXAGON():
+
+ inputfile="./meshHexagonWithTriangles300.med";
+ my_mesh=cdmath.Mesh(inputfile);
+ spaceDim=2
+
+ # Initial field data
+ nVar=2+spaceDim;
+ radius=0.5;
+ Center=cdmath.Vector(spaceDim);#default value is (0,0,0)
+ Vout=cdmath.Vector(nVar)
+ Vin =cdmath.Vector(nVar)
+ Vin[0]=155e5;
+ Vin[1]=0;
+ Vin[2]=0;
+ Vin[3]=563;
+ Vout[0]=154e5;
+ Vout[1]=0;
+ Vout[2]=0;
+ Vout[3]=563;
+
+ myProblem = cf.SinglePhase(cf.Liquid,cf.around155bars600K,spaceDim);
+
+ # Initial field creation
+ print ("Setting mesh and initial data" ) ;
+ myProblem.setInitialFieldSphericalStepFunction( my_mesh, Vout, Vin, radius, Center);
+
+ # set the boundary conditions
+ wallVelocityX=0;
+ wallVelocityY=0;
+ wallTemperature=563;
+
+ myProblem.setWallBoundaryCondition("boundaries", wallTemperature, wallVelocityX, wallVelocityY);
+
+ # set the numerical method
+ myProblem.setNumericalScheme(cf.upwind, cf.Explicit);
+
+ # name file save
+ fileName = "2DSphericalExplosion_HEXAGON";
+
+ # parameters calculation
+ MaxNbOfTimeStep = 3000 ;
+ freqSave = 20;
+ cfl = 0.5;
+ maxTime = 5;
+ precision = 1e-6;
+
+ myProblem.setCFL(cfl);
+ myProblem.setPrecision(precision);
+ myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
+ myProblem.setTimeMax(maxTime);
+ myProblem.setFreqSave(freqSave);
+ myProblem.setFileName(fileName);
+ myProblem.setNewtonSolver(precision,20);
+ myProblem.saveConservativeField(False);
+ if(spaceDim>1):
+ myProblem.saveVelocity();
+ pass
+
+
+ # evolution
+ myProblem.initialize();
+
+ ok = myProblem.run();
+ if (ok):
+ print( "Simulation python " + fileName + " is successful !" );
+ pass
+ else:
+ print( "Simulation python " + fileName + " failed ! " );
+ pass
+
+ print( "------------ End of calculation !!! -----------" );
+
+ myProblem.terminate();
+ return ok
+
+if __name__ == """__main__""":
+ SinglePhase_2DSphericalExplosion_HEXAGON()
