Commit 7f8977af authored by David Vincent's avatar David Vincent Committed by Valentin Vallaeys
Browse files

stommel testcase is now implemented via the interface

a new testcase is implemented in order to test the function restart
the interface (slim.restart)

modifying slimPre.py and slim.py such as the user can set the polynomial order of
discretization in slimPre.mesh and in slim.domain (the default value is
1)
parent aa6f7c11
Pipeline #1710 failed with stage
in 30 minutes and 47 seconds
......@@ -53,7 +53,7 @@ import dgpy
class Domain:
"""Create the numerical domain"""
def __init__(self, mesh, bathy, g=None, density=1025, solve_on_sphere=False):
def __init__(self, mesh, bathy, g=None, density=1025, solve_on_sphere=False,order=1):
"""keyword arguments:
* mesh
......@@ -66,6 +66,8 @@ class Domain:
density of the liquid (default: 1025 kg/m^3)
* solve_on_sphere
boolean stating if you want to solve the shallow water equation on a sphere or on any 2d manifold instead of solving it on a plane (default: False)
* order
Polynomial order of the discretization
"""
self._solve_on_manifold = solve_on_sphere
if (dgpy.Msg.GetCommSize() > 1) and (not slim_private._check_if_partionned(mesh)):
......@@ -73,19 +75,18 @@ class Domain:
self._mesh = dgpy.dgMesh(mesh[:-4] + '_' + str(dgpy.Msg.GetCommSize()) + '.msh')
else:
self._mesh = dgpy.dgMesh(mesh)
#TODO add option to run in high order
order = 1
self._groups = dgpy.dgGroupCollection(self._mesh, order)
self._order = order
self._groups = dgpy.dgGroupCollection(self._mesh, self._order)
if self._solve_on_manifold:
self._space_transform = dgpy.dgSpaceTransformFlatTriangle(self._mesh);
self._mesh.setSpaceTransform(self._space_transform, self._groups)
self._bath = dgpy.dgDofContainer(self._groups, 1)
slim_private._load(self._bath, bathy)
self._bath_PC = dgpy.functionPrecomputed(self._groups, 3)
self._bath_PC = dgpy.functionPrecomputed(self._groups, 2*self._order+1)
self._bath_PC.compute(self._bath.getFunction())
self._bath_PC.compute(self._bath.getFunction(),dgpy.dataCacheMap.NODE_GROUP_MODE)
self._bath_gradient = self._bath.getFunctionGradient()
self._bath_gradient_PC = dgpy.functionPrecomputed(self._groups, 3)
self._bath_gradient_PC = dgpy.functionPrecomputed(self._groups, 2*self._order+1)
self._bath_gradient_PC.compute(self._bath.getFunctionGradient())
self._bath_function = self._bath_PC
self._bath_gradient_function = self._bath_gradient_PC
......@@ -211,7 +212,7 @@ class ShallowWater2dWD:
netcdf or .msh file containing the coriolis term for the whole domain.
"""
self._coriolis = slim_private._load_function(coriolis, self._domain._groups)
self._coriolis_PC = dgpy.functionPrecomputed(self._domain._groups, 3)
self._coriolis_PC = dgpy.functionPrecomputed(self._domain._groups, self._domain._order*2+1)
self._coriolis_PC.compute(self._coriolis)
self._equation.setCoriolisFactor(self._coriolis_PC)
......@@ -1272,7 +1273,7 @@ class Loop:
index = self._n_iter_offline
i._hydro_sol_dof.importMsh(i._datafile+"offline_sw2d-%06d_%06d/offline_sw2d-%06d_%06d" %(index,j,index,j))
i._hydro_sol_dof.scatter()
newton_converged = i._temporal_solver.subiterate(i._solution, self._dt, self._time)
newton_converged = i._temporal_solver.subiterate(i._solution, self._dt, self._time)
if not newton_converged:
dgpy.Msg.Fatal("Solver failed")
if offline and i._temporal_scheme == "explicit":
......
......@@ -80,7 +80,7 @@ import numpy as np
class Mesh :
"""Create the slimPre Mesh"""
def __init__(self, mesh_file_name, mesh_proj=None) :
def __init__(self, mesh_file_name, mesh_proj=None, order=1) :
"""keyword arguments:
* mesh_file_name
......@@ -89,13 +89,15 @@ class Mesh :
mesh projection system (Default: None).
This argument must not be given if the mesh is written in 3D (on a sphere), or if all the data are in the same projection system than the mesh.
* order
Polynomial order of the discretization
"""
isPart = slim_private._check_if_partionned(mesh_file_name)
if (dgpy.Msg.GetCommSize() > 1) and not isPart:
dgpy.Msg.Fatal("Calling slimPre.Mesh in multiprocessing requires providing a partitioned mesh (first use slimPre.partition_mesh).")
# elif (dgpy.Msg.GetCommSize() == 1) and isPart:
# dgpy.Msg.Fatal("Calling slimPre.Mesh on a single core requires providing a non-partitioned mesh.")
self._groups = dgpy.dgGroupCollection(mesh_file_name)
self._groups = dgpy.dgGroupCollection(mesh_file_name,order)
self._mesh_proj = mesh_proj
......
l = 5e5;
Point(1) = {-l, -l, 0};
Point(2) = {l, -l, 0};
Point(3) = {l, l, 0};
Point(4) = {-l, l, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {3, 4, 1, 2};
Field[1] = MathEval;
Field[1].F = "1e5";
Background Field = 1;
//Field[1] = MathEval;
//Field[1].F = "0.5e4+1*(abs(x)+abs(y))";
//Background Field = 1;
Plane Surface(6) = {5};
Physical Surface("Surface")={6};
Physical Line("Wall") = {1,2,3,4};
Mesh.CharacteristicLengthExtendFromBoundary=1;
Mesh.CharacteristicLengthFromPoints=1;
Mesh.SecondOrderIncomplete=0;
##################################################
########## PREPROCESSING ##########
#################################################
import slimPre
import numpy as np
from dgpy.scripts import Common
#Polynomial order of the discretization
#time step
dt = 2e5
nbSteps = 5
#mesh
Common.genMesh("stommel_square.geo", 2)
mesh = slimPre.Mesh("stommel_square.msh",order=1)
xyz = mesh._groups.getFunctionCoordinates()
print("start preprocessing")
#bathymetry
bathymetryDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
bathy = slimPre.dgpy.functionConstant(1000)
bathymetryDof.interpolate(bathy)
bathymetryDof.exportMsh("bathymetry")
#gravity
gravityDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
gravity = slimPre.dgpy.functionConstant(9.80616)
gravityDof.interpolate(gravity)
gravityDof.exportMsh("gravity")
#bottom friction
gammaDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
gamma = slimPre.dgpy.functionConstant(1e-6)
gammaDof.interpolate(gamma)
gammaDof.exportMsh("gamma")
#coriolis
def corioF(cm):
coord = cm.get(xyz)
nPt = coord.shape
val = np.zeros([nPt[0],1])
val[:,0] = 1e-4+2e-11*coord[:,1]
return val
corioDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
corio = slimPre.dgpy.functorNumpy(corioF)
corioDof.interpolate(corio)
corioDof.exportMsh("coriolis")
#wind
def windF(cm):
coord = cm.get(xyz)
nPt = coord.shape
val = np.zeros([nPt[0],2])
val[:,0] = 0.1* np.sin(np.pi*coord[:,1]/1e6)/1e6
return val
windXDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
windX = slimPre.dgpy.functorNumpy(windF)
windXDof.interpolate(windX)
windXDof.exportMsh("windX")
windYDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
windY = slimPre.dgpy.functionConstant(0.)
windYDof.interpolate(windY)
windYDof.exportMsh("windY")
print("End preprocessing")
#################################################
########## SIMULATION ##########
#################################################
import slim
print("start Simulation")
domain = slim.Domain("stommel_square.msh","bathymetry_COMP_0.msh", g= "gravity_COMP_0.msh")
xyz = domain._groups.getFunctionCoordinates()
eq = slim.ShallowWater2d(domain, "implicit", linear_equation=True, initial_time=0, final_time=nbSteps*dt)
eq.set_forcing(forcing_x="windX_COMP_0.msh", forcing_y="windY_COMP_0.msh")
eq.set_viscosity("constant", constant_viscosity=0.00)
eq.set_dissipation("linear","gamma_COMP_0.msh")
eq.set_coriolis("coriolis_COMP_0.msh")
eq.set_boundary_coast("Wall", slip=False)
loop=slim.Loop(maximum_time_step = dt, export_time = dt,path="./output")
loop.add_equation(eq)
loop._Rtol = 1e-6
loop._Atol = 1e-6
loop.run()
print("end Simulation")
sol = slimPre.dgpy.dgDofContainer(domain._groups,3)
sol.importIdx("output/sw2d/eta/eta-000005.idx", {0:0})
sol.importIdx("output/sw2d/uv/uv-000005.idx", {1:0,2:1})
loop.restart(2,2*dt)
solRestart = slimPre.dgpy.dgDofContainer(domain._groups,3)
solRestart.importIdx("output/sw2d/eta/eta-000005.idx", {0:0})
solRestart.importIdx("output/sw2d/uv/uv-000005.idx", {1:0,2:1})
def errorRestart(cm):
solution = cm.get(sol.getFunction())
solutionR = cm.get(solRestart.getFunction())
return (solution-solutionR)**2
print("start post processing")
errorRestartF = slimPre.dgpy.functorNumpy(errorRestart)
integrator = slimPre.dgpy.dgFunctionIntegrator(domain._groups, errorRestartF)
integral = slimPre.dgpy.fullMatrixDouble(3,1)
integrator.compute(integral)
errorU = [np.sqrt(integral(0,0))/1e12, np.sqrt(integral(1,0))/1e12, np.sqrt(integral(2,0))/1e12]
print ("L2 error on (eta, u, v): (%e, %e, %e)\n" % (errorU[0], errorU[1], errorU[2]))
#Abnormal termination if the error is too high
if (errorU[0] + errorU[1] + errorU[2] > 1e-12):
slimPre.dgpy.Msg.Error("Error is too high");
slimPre.dgpy.Msg.Exit(1)
slimPre.dgpy.Msg.Barrier()
if (slimPre.dgpy.Msg.GetCommRank() == 0):
print( '')
print('******************** End ********************')
print( '')
slimPre.dgpy.Msg.Exit(0)
\ No newline at end of file
#!/usr/bin/env rundgpy
#We import gmsg dg module to interact with the c++ library
#c++ functions from gmsh can now be called in the python code
from dgpy import *
#We import common python packages
import os, math
##################################################
########## PREPROCESSING ##########
#################################################
import slimPre
import numpy as np
from dgpy.scripts import Common
#A GModel is a mesh database
#The mesh "stommel_square.msh" is generated using the geometry file "stommel_square.geo"
#The generated mesh is 2d (second argument) and of order 1 (third argument)
#Common.genMesh('stommel_square', 2, 1)
#The mesh is partitionned for the number of processes
#model.load(gmshPartition.simple('stommel_square.msh'))
#Polynomial order of the discretization
order = 3
#Groups are used to differentiate different elements of the mesh
#Surface elements and edges belong to different groups
#Group can be split into sub-groups (i.e. split groups of elements that belong to searate physical entities)
Common.genMesh("stommel_square.geo", 2, 1, Msg.GetCommSize())
groups = dgGroupCollection("stommel_square.msh", order)
groups.splitGroupsRandom(5)
#Create directory for output
try : os.mkdir("output")
except: pass
order=3
#time step
dt = 2e5
nbSteps = 60
#mesh
Common.genMesh("stommel_square.geo", 2)
mesh = slimPre.Mesh("stommel_square.msh",order=order)
xyz = mesh._groups.getFunctionCoordinates()
#For functions that need to be fast (often called), c++ code can be included in the python files and compiled during runtime
#This is an example of c++ functions used to define the forcings (wind), coriolis and the analytical solution
#Not that the analytical solution could be written in python as it is called only once
CCode = """
#include "fullMatrix.h"
#include "function.h"
extern "C" {
void wind (dataCacheMap *, fullMatrix<double> &sol, fullMatrix<double> &xyz) {
for (size_t i = 0; i< sol.size1(); i++) {
sol.set(i,0,0.1* sin(M_PI*xyz(i,1)/1e6)/1e6);
sol.set(i,1,0);
}
}
void coriolis (dataCacheMap *, fullMatrix<double> &sol, fullMatrix<double> &xyz) {
for (size_t i = 0; i< sol.size1(); i++) {
sol.set(i,0,1e-4+2e-11*(xyz(i,1)));
}
}
void UVEtaAnalytic (dataCacheMap *, fullMatrix<double> &val, fullMatrix<double> &xyz) {
double L = 1e6;
double tau0 = 0.1;
double rho = 1000;
double gamm = 1e-6;
double delta = 1;
double g = 9.80616;
double h = 1000;
double f0 = 1e-4;
double beta = 2e-11;
double epsilon = gamm / (L * beta);
double z1 = (-1 + sqrt(1 + (2 * M_PI * delta * epsilon) * (2 * M_PI * delta * epsilon))) / (2 * epsilon);
double z2 = (-1 - sqrt(1 + (2 * M_PI * delta * epsilon) * (2 * M_PI * delta * epsilon))) / (2 * epsilon);
double D = ((exp(z2) - 1) * z1 + (1 - exp(z1)) * z2) / (exp(z1) - exp(z2));
for (size_t i = 0; i< val.size1(); i++) {
double x = xyz(i,0) / L + 0.5;
double y = xyz(i,1) / L;
double f1 = M_PI / D * (1 + ((exp(z2) - 1) * exp(x * z1) + (1 - exp(z1)) * exp(x * z2)) / (exp(z1) - exp(z2)));
double f2 = 1 / D * (((exp(z2) - 1) * z1 * exp(x * z1) + (1 - exp(z1)) * z2 * exp(x * z2)) / (exp(z1) - exp(z2)));
//eta
val(i, 0) = D * tau0 * f0 * L / (M_PI * gamm * rho * delta * g * h) * (-gamm / (f0 * delta * M_PI) * f2 * sin(M_PI * y) +
1 / M_PI * f1 * (cos(M_PI * y) * (1 + beta * y) - beta / M_PI * sin(M_PI * y)));
//u
val(i, 1) = D * tau0 / (M_PI * gamm * rho * h) * f1 * sin(M_PI * y);
//v
val(i, 2) = D * tau0 / (M_PI * gamm * rho * delta * h) * f2 * cos(M_PI * y);
//w
val(i, 3) = 0;
}
}
}
"""
print("start preprocessing")
#bathymetry
bathymetryDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
bathy = slimPre.dgpy.functionConstant(1000)
bathymetryDof.interpolate(bathy)
bathymetryDof.exportMsh("bathymetry")
#The c++ code is compiled
CCodeLib = "./CCode.dylib"
if (Msg.GetCommRank() == 0 ) :
functionC.buildLibrary (CCode, CCodeLib);
Msg.Barrier()
#gravity
gravityDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
gravity = slimPre.dgpy.functionConstant(9.80616)
gravityDof.interpolate(gravity)
gravityDof.exportMsh("gravity")
#Here we define the conservation law to use. Thera are several conservation laws
#available (see the files beginning with dgConservationLaw)
#Two-dimensional shallow water are considered here
claw = dgConservationLawShallowWater2d()
#The container that contains the degrees of freedom
#For the shallow water equations, claw.getNbFields() is 3 (u, v and eta)
solution = dgDofContainer(groups, claw.getNbFields())
solution.setAll(0.)
XYZ = groups.getFunctionCoordinates();
#In this example, we consider the linear shallow water equation
#(availability of the analytical solution)
claw.setIsLinear(True)
#We use the functions defined in c++ above and pass it to the
#conservation law
coriolis = functionC(CCodeLib,"coriolis",1,[XYZ])
claw.setCoriolisFactor(coriolis)
wind = functionC(CCodeLib,"wind",2,[XYZ])
claw.setSource(wind)
#Various parameters: gamma is the linear dissipation coefficient and bath is the water depth at rest
gamma = functionConstant([1e-6])
claw.setLinearDissipation(gamma)
bath = functionConstant([1000])
claw.setBathymetry(bath)
#We enforce the boundary conditions
claw.addBoundaryCondition('Wall',claw.newBoundaryWall())
#The linear system and its solver
#Here we choose a block structured system implemented in petsc
#We pass additional options to petsc (lu in serial, ilu in parallel)
sys = linearSystemPETScBlockDouble ()
if(Msg.GetCommSize() > 1):
sys.setParameter("petscOptions", "-pc_type bjacobi -sub_pc_type lu")
else:
sys.setParameter("petscOptions", "-pc_type lu")
dof = dgDofManager.newDGBlock (groups, 3, sys)
#We can choose severa time discretization. Here an implicit euler,
#as we as interested in a stationary solution
implicitEuler = dgDIRK(claw, dof, 1)
implicitEuler.getNewton().setAtol(1e-6)
implicitEuler.getNewton().setRtol(1e-6)
implicitEuler.getNewton().setVerb(10)
#bottom friction
gammaDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
gamma = slimPre.dgpy.functionConstant(1e-6)
gammaDof.interpolate(gamma)
gammaDof.exportMsh("gamma")
#The time integration loop is written in the python file
#Here, 60 iterations are considered
dt = 2e5
t = 0.0
solution.exportMsh("output/stommel-%06d" % (0), 0, 0) #export initial solution
for i in range (60) :
t = t+dt
print('|ITER|',i,'|DT|',dt, '|t|', t)
implicitEuler.iterate(solution, dt, t)
solution.exportMsh("output/stommel-%06d" % (i+1), t, i+1) #export solution after each time step
#coriolis
def corioF(cm):
coord = cm.get(xyz)
nPt = coord.shape
val = np.zeros([nPt[0],1])
val[:,0] = 1e-4+2e-11*coord[:,1]
return val
corioDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
corio = slimPre.dgpy.functorNumpy(corioF)
corioDof.interpolate(corio)
corioDof.exportMsh("coriolis")
#wind
def windF(cm):
coord = cm.get(xyz)
nPt = coord.shape
val = np.zeros([nPt[0],2])
val[:,0] = 0.1* np.sin(np.pi*coord[:,1]/1e6)/1e6
return val
windXDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
windX = slimPre.dgpy.functorNumpy(windF)
windXDof.interpolate(windX)
windXDof.exportMsh("windX")
windYDof = slimPre.dgpy.dgDofContainer(mesh._groups,1)
windY = slimPre.dgpy.functionConstant(0.)
windYDof.interpolate(windY)
windYDof.exportMsh("windY")
print("End preprocessing")
#################################################
########## SIMULATION ##########
#################################################
import slim
#We use a function to define the error versus the analytical solution (the latest being written above in c++)
#Then, we integrate this error over the domain
#Elevation eta is defined up to a constant. Hence, we do not compute error for eta.
print("start Simulation")
domain = slim.Domain("stommel_square.msh","bathymetry_COMP_0.msh", g= "gravity_COMP_0.msh",order=order)
eq = slim.ShallowWater2d(domain, "implicit", linear_equation=True, initial_time=0, final_time=nbSteps*dt)
eq.set_forcing(forcing_x="windX_COMP_0.msh", forcing_y="windY_COMP_0.msh")
eq.set_viscosity("constant", constant_viscosity=0.00)
eq.set_dissipation("linear","gamma_COMP_0.msh")
eq.set_coriolis("coriolis_COMP_0.msh")
eq.set_boundary_coast("Wall", slip=False)
loop=slim.Loop(maximum_time_step = dt, export_time = dt,path="./output")
loop.add_equation(eq)
loop._Rtol = 1e-6
loop._Atol = 1e-6
loop.run()
print("end Simulation")
xyz = domain._groups.getFunctionCoordinates()
#####################################################################################################################################
def UVEtaAnalytic(cm):
L = 1e6
tau0 = 0.1
rho = 1000
gamm = 1e-6
delta = 1
g = 9.80616
h = 1000
f0 = 1e-4
beta = 2e-11
epsilon = gamm / (L * beta)
z1 = (-1 + np.sqrt(1 + (2 * np.pi * delta * epsilon) * (2 * np.pi * delta * epsilon))) / (2 * epsilon)
z2 = (-1 - np.sqrt(1 + (2 * np.pi * delta * epsilon) * (2 * np.pi * delta * epsilon))) / (2 * epsilon)
D = ((np.exp(z2) - 1) * z1 + (1 - np.exp(z1)) * z2) / (np.exp(z1) - np.exp(z2))
xyzA= cm.get(xyz)
nPt = xyzA.shape
val = np.zeros([nPt[0],3])
f1 = np.zeros([nPt[0],1])
f2 = np.zeros([nPt[0],1])
x = xyzA[:,0] / L + 0.5
y = xyzA[:,1]/ L
f1[:,0] = np.pi / D * (1 + ((np.exp(z2) - 1) * np.exp(x * z1) + (1 - np.exp(z1)) * np.exp(x * z2)) / (np.exp(z1) - np.exp(z2)))
f2[:,0] = 1 / D * (((np.exp(z2) - 1) * z1 * np.exp(x * z1) + (1 - np.exp(z1)) * z2 * np.exp(x * z2)) / (np.exp(z1) - np.exp(z2)))
for i in range(nPt[0]):
#eta
val[i,0] = D * tau0 * f0 * L / (np.pi * gamm * rho * delta * g * h) * (-gamm / (f0 * delta * np.pi) * f2[i] * np.sin(np.pi * y[i]) + 1. / np.pi * f1[i] * (np.cos(np.pi * y[i]) * (1 + beta * y[i]) - beta / np.pi * np.sin(np.pi * y[i])))
#u
val[i,1] = D * tau0 / (np.pi * gamm * rho * h) * f1[i] * np.sin(np.pi * y[i])
#u
val[i,2] = D * tau0 / (np.pi * gamm * rho * delta * h) * f2[i] * np.cos(np.pi * y[i])
return val
def errorAnalytical(cm):
sol = cm.get(solution.getFunction())
sol = cm.get(eq._solution.getFunction())
ana = cm.get(UVEtaAnalytic)
return (sol - ana[:,:3])**2
UVEtaAnalytic = functionC(CCodeLib,"UVEtaAnalytic",4,[XYZ])
errorAnalyticalF = functorNumpy(errorAnalytical)
integrator = dgFunctionIntegrator(groups, errorAnalyticalF)
integral = fullMatrixDouble(3,1)
return (sol - ana)**2
print("start post processing")
UVEtaAnalytic = slimPre.dgpy.functorNumpy(UVEtaAnalytic)
errorAnalyticalF = slimPre.dgpy.functorNumpy(errorAnalytical)
integrator = slimPre.dgpy.dgFunctionIntegrator(domain._groups, errorAnalyticalF)
integral = slimPre.dgpy.fullMatrixDouble(3,1)
integrator.compute(integral)
errorU = [math.sqrt(integral(1,0))/1e12, math.sqrt(integral(2,0))/1e12]
errorU = [np.sqrt(integral(1,0))/1e12, np.sqrt(integral(2,0))/1e12]
print ("L2 error on (u, v): (%e, %e)\n" % (errorU[0], errorU[1]))
#Abnormal termination if the error is too high
if (errorU[0] + errorU[1] > 4e-10):
Msg.Error("Error is too high");
Msg.Exit(1)
Msg.Exit(0)
slimPre.dgpy.Msg.Error("Error is too high");
slimPre.dgpy.Msg.Exit(1)
slimPre.dgpy.Msg.Barrier()
if (slimPre.dgpy.Msg.GetCommRank() == 0):
print( '')
print('******************** End ********************')
print( '')
slimPre.dgpy.Msg.Exit(0)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment