From: michael Date: Mon, 23 Nov 2020 11:52:26 +0000 (+0100) Subject: Saved the numerical results if the spherical explosion in an hexagonal cavity X-Git-Tag: V9_7_0~141 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=dcde305b057accb4e2a4a403d6242f331a13522a;p=tools%2Fsolverlab.git Saved the numerical results if the spherical explosion in an hexagonal cavity --- diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureStaggered.png b/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureStaggered.png deleted file mode 100644 index f689cff..0000000 Binary files a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureStaggered.png and /dev/null differ diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindExplicit.png b/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindExplicit.png deleted file mode 100644 index c910314..0000000 Binary files a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindExplicit.png and /dev/null differ diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindExplicitScaled.png b/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindExplicitScaled.png deleted file mode 100644 index a125c4b..0000000 Binary files a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindExplicitScaled.png and /dev/null differ diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindImplicit.png b/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindImplicit.png deleted file mode 100644 index db66c42..0000000 Binary files a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindImplicit.png and /dev/null differ diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindImplicitScaled.png b/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindImplicitScaled.png deleted file mode 100644 index fd8fbf9..0000000 Binary files a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalPressureUpwindImplicitScaled.png and /dev/null differ diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalUpwindExplicitImplicit.png b/CDMATH/tests/doc/2DWaveSystemShockWave/FinalUpwindExplicitImplicit.png deleted file mode 100644 index 46ef81c..0000000 Binary files a/CDMATH/tests/doc/2DWaveSystemShockWave/FinalUpwindExplicitImplicit.png and /dev/null differ diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/HEXAGON/mainCDMATH.cxx b/CDMATH/tests/doc/2DWaveSystemShockWave/HEXAGON/mainCDMATH.cxx new file mode 100755 index 0000000..82f70c2 --- /dev/null +++ b/CDMATH/tests/doc/2DWaveSystemShockWave/HEXAGON/mainCDMATH.cxx @@ -0,0 +1,276 @@ +//============================================================================ +// Author : Michael NDJINGA +// Date : November 2020 +// Description : 2D linear wave system +//============================================================================ + +#include +#include +#include + +#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=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 or isStationary or time >=tmax) + { + cout<<"-- Iteration: " << it << ", Time: " << time << ", dt: " << dt<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<=ntmax) + cout<< "Nombre de pas de temps maximum ntmax= "<< ntmax<< " atteint"<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 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) diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/SQUARE/WaveSystemUpwindExplicite.py b/CDMATH/tests/doc/2DWaveSystemShockWave/SQUARE/WaveSystemUpwindExplicite.py new file mode 100644 index 0000000..cd974e6 --- /dev/null +++ b/CDMATH/tests/doc/2DWaveSystemShockWave/SQUARE/WaveSystemUpwindExplicite.py @@ -0,0 +1,278 @@ +#!/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 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) diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/SQUARE/WaveSystemUpwindImplicite.py b/CDMATH/tests/doc/2DWaveSystemShockWave/SQUARE/WaveSystemUpwindImplicite.py new file mode 100644 index 0000000..a4b2333 --- /dev/null +++ b/CDMATH/tests/doc/2DWaveSystemShockWave/SQUARE/WaveSystemUpwindImplicite.py @@ -0,0 +1,278 @@ +#!/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 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) diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/WaveSystemStaggered.py b/CDMATH/tests/doc/2DWaveSystemShockWave/WaveSystemStaggered.py deleted file mode 100644 index cefdd1e..0000000 --- a/CDMATH/tests/doc/2DWaveSystemShockWave/WaveSystemStaggered.py +++ /dev/null @@ -1,336 +0,0 @@ -#!/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 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) diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/WaveSystemUpwindExplicite.py b/CDMATH/tests/doc/2DWaveSystemShockWave/WaveSystemUpwindExplicite.py deleted file mode 100644 index cd974e6..0000000 --- a/CDMATH/tests/doc/2DWaveSystemShockWave/WaveSystemUpwindExplicite.py +++ /dev/null @@ -1,278 +0,0 @@ -#!/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 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) diff --git a/CDMATH/tests/doc/2DWaveSystemShockWave/WaveSystemUpwindImplicite.py b/CDMATH/tests/doc/2DWaveSystemShockWave/WaveSystemUpwindImplicite.py deleted file mode 100644 index a4b2333..0000000 --- a/CDMATH/tests/doc/2DWaveSystemShockWave/WaveSystemUpwindImplicite.py +++ /dev/null @@ -1,278 +0,0 @@ -#!/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 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) diff --git a/CoreFlows/Documentation/Figures/Simulations/SphericalExploxion_HEXAGON/FinalState_HEXAGON300.png b/CoreFlows/Documentation/Figures/Simulations/SphericalExploxion_HEXAGON/FinalState_HEXAGON300.png new file mode 100644 index 0000000..22355b1 Binary files /dev/null and b/CoreFlows/Documentation/Figures/Simulations/SphericalExploxion_HEXAGON/FinalState_HEXAGON300.png differ diff --git a/CoreFlows/Documentation/Figures/Simulations/SphericalExploxion_HEXAGON/SinglePhase_2DSphericalExplosion_HEXAGON.py b/CoreFlows/Documentation/Figures/Simulations/SphericalExploxion_HEXAGON/SinglePhase_2DSphericalExplosion_HEXAGON.py new file mode 100755 index 0000000..526bcf8 --- /dev/null +++ b/CoreFlows/Documentation/Figures/Simulations/SphericalExploxion_HEXAGON/SinglePhase_2DSphericalExplosion_HEXAGON.py @@ -0,0 +1,85 @@ +#!/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() diff --git a/CoreFlows/Documentation/Figures/Simulations/SphericalExploxion_HEXAGON/pressureWaveHEXAGON300.ogv b/CoreFlows/Documentation/Figures/Simulations/SphericalExploxion_HEXAGON/pressureWaveHEXAGON300.ogv new file mode 100644 index 0000000..84ca9ac Binary files /dev/null and b/CoreFlows/Documentation/Figures/Simulations/SphericalExploxion_HEXAGON/pressureWaveHEXAGON300.ogv differ