Commit a5614153 authored by thomas's avatar thomas
Browse files

Add realistic GBR test case

git-svn-id: https://geuz.org/svn/gmsh/trunk@13428 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent c7d75204
# Script to generate smooth bathymetry.
# Filename is set in InputGBR.py
# CARTESIAN
#NB: bathy class defined in dgParticleTracker2D.h
from dgpy import *
from ProjectMesh import *
import os
#-------------
filename = "gbr10K"
#-------------
# Physical parameters for bathymetry smoothing
bNu=20 # Diffusion
bT=10000 # Time for smoothing
#-------------
print""
print"Generating smooth bathymetry for CARTESIAN PLANE mesh:", filename
#First project mesh in planar space
if(not os.path.exists("./Meshes/"+filename+"_tan.msh")):
print'Projecting mesh in planar space ...'
m = GModel()
m.load("./Meshes/"+filename+".msh")
projectMesh(m)
m.save("./Meshes/"+filename+"_tan.msh")
model=GModel()
model.load("Meshes/"+filename+"_tan.msh")
#Compile libraries
if (Msg.GetCommRank() == 0):
functionC.buildLibraryFromFile ("./Libraries/GBR.cc", "./Libraries/lib_gbr.so")
glib = "./Libraries/lib_gbr.so"
#Get started
if(Msg.GetCommSize() > 1):
Msg.Fatal("You need to run the script a first time on only one processor before you can run it on more!")
Msg.Info( "Diffuse bathymetry")
#Create output folder
try : os.mkdir('./Bath');
except: 0;
XYZ = function.getCoordinates()
lonLatDegrees = functionC(glib,"lonLatDegrees",3,[XYZ])
groups = dgGroupCollection(model,2,1)
bathDC = dgDofContainer(groups,1);
#OLD:
#bathS = functionC(blib, "bath", 1, [stereoToLonLatDegrees])
#bathDC.interpolate(bathS);
#NEW:
#1. Create a bathy object (initialises object and creates bathymetry function)
bathymetryObj = bathy("./Data/BathHR.bin", lonLatDegrees)
#2. Interpolate the function onto DC
bathDC.interpolate(bathymetryObj)
#3. Delete bathymetry object data
bathymetryObj.clearBathData()
#Diffuse & save:
nuB = functionConstant(bNu)
dlaw = dgConservationLawAdvectionDiffusion.diffusionLaw(nuB)
dlaw.addBoundaryCondition("Coast", dlaw.new0FluxBoundary())
dlaw.addBoundaryCondition("Islands", dlaw.new0FluxBoundary())
dlaw.addBoundaryCondition("Open Sea North", dlaw.new0FluxBoundary())
dlaw.addBoundaryCondition("Open Sea Central", dlaw.new0FluxBoundary())
dlaw.addBoundaryCondition("Open Sea South", dlaw.new0FluxBoundary())
sys = linearSystemPETScDouble()
dof = dgDofManager.newCG (groups, 1, sys)
implicitEuler = dgDIRK(dlaw, dof)
implicitEuler.iterate (bathDC, bT, 0)
exporterBath=dgIdxExporter(bathDC, "./Bath/"+filename+"_bath_smooth");
exporterBath.exportIdx()
sys = 0
dof = 0
Msg.Exit(0)
# Script to generate smooth bathymetry.
# Filename is set in InputGBR.py
# STEREOGRAPHIC
#NB: bathy class defined in dgParticleTracker2D.h
from dgpy import *
import os
#-------------
filename = "gbr20K_S_sgbr"
#-------------
# Physical parameters for bathymetry smoothing
bNu=20 # Diffusion
bT=10000 # Time for smoothing
#-------------
print""
print"Generating smooth bathymetry for STEREO mesh:", filename
model=GModel()
model.load("Meshes/"+filename+".msh")
#Compile libraries
if (Msg.GetCommRank() == 0):
functionC.buildLibraryFromFile ("./Libraries/GBR.cc", "./Libraries/lib_gbr.so")
glib = "./Libraries/lib_gbr.so"
#Get started
if(Msg.GetCommSize() > 1):
Msg.Fatal("You need to run the script a first time on only one processor before you can run it on more!")
Msg.Info( "Diffuse bathymetry")
XYZ = function.getCoordinates()
stereoToLonLatDegrees = functionC(glib,"stereoToLonLatDegrees",3,[XYZ])
groups = dgGroupCollection(model,2,1)
bathDC = dgDofContainer(groups,1);
#OLD:
#bathS = functionC(blib, "bath", 1, [stereoToLonLatDegrees])
#bathDC.interpolate(bathS);
#NEW:
#1. Create a bathy object (initialises object and creates bathymetry function)
bathymetryObj = bathy("./Data/BathHR.bin", stereoToLonLatDegrees)
#2. Interpolate the function onto DC
bathDC.interpolate(bathymetryObj)
#3. Delete bathymetry object data
bathymetryObj.clearBathData()
#Diffuse & save:
nuB = functionConstant(bNu)
dlaw = dgConservationLawAdvectionDiffusion.diffusionLaw(nuB)
dlaw.addBoundaryCondition("Coast", dlaw.new0FluxBoundary())
dlaw.addBoundaryCondition("Island", dlaw.new0FluxBoundary())
dlaw.addBoundaryCondition("OpenSeaNorth", dlaw.new0FluxBoundary())
dlaw.addBoundaryCondition("OpenSeaCentral", dlaw.new0FluxBoundary())
dlaw.addBoundaryCondition("OpenSeaSouth", dlaw.new0FluxBoundary())
sys = linearSystemPETScDouble()
dof = dgDofManager.newCG (groups, 1, sys)
implicitEuler = dgDIRK(dlaw, dof)
implicitEuler.iterate (bathDC, bT, 0)
exporterBath=dgIdxExporter(bathDC, "./Bath/"+filename+"_bath_smooth");
exporterBath.exportIdx()
sys = 0
dof = 0
Msg.Exit(0)
\ No newline at end of file
from dgpy import *
from ProjectMesh import *
import calendar, os, time
from math import sqrt
#================
# 1. CHOOSE MESH
#----------------
mesh = "10K"
#mesh = "100K"
#================
#================
# 2. SET FILE NAME & OUTPUT DIRECTORY
#----------------
filename = "gbr"+mesh
outputDir = "output"
#================
#================
#3. SET BATHYMETRY & REEF MAP DIRECTORIES
#----------------
bathdir=outputDir+"/Bath/"+filename+"_bath_smooth"
bathname=outputDir+"/Bath/"+filename+"_bath_smooth"+"/"+filename+"_bath_smooth.idx"
reefMapFile="Data/ReefMap_9200.dat"
#================
#================
#4. SET EXPORT PARAMETERS
#----------------
# Print every "iter" iterations, export every "export" iterations, export from t=t_exportStart
iter = 1
export = 1
export_fnEval = 1
t_exportStart = 0 # 6*24*3600.0
#================
# Mesh Parameters
order = 1
dimension = 2
#================
# 5. SET INITIAL AND FINAL TIME [year, month, day, hour, minute, second]
#----------------
Ti = calendar.timegm([2007, 11, 28, 0, 0, 0])
#Tf = calendar.timegm([2000, 4, 1, 1, 0, 0])
Tf = Ti + 36*24*3600.0
#================
#================
# 6. OPTIONAL: LOAD INITIAL SOLUTION FROM EXISTING FILEi
# NB: must also set, above, Ti = Ti of old simulation
#----------------
# startFromExistingSolution: 0 if no, 1 if yes
startFromExistingSolution = 0
stepNbOfInitialSolution = 53360
solutionFileName = '/gbr-%06d'%( stepNbOfInitialSolution )+'.idx'
old_dt = 1
if (startFromExistingSolution == 1):
Ti = Ti + stepNbOfInitialSolution*old_dt
Tf = Ti + 7*24*3600
#================
#================
# 7. OTHER PARAMETERS
#----------------
# Radius of the Earth (for stereo projection)
EarthRadius = 6371e3
# Multirate and Partition Parameters
RK_TYPE = ERK_22_C
mL = 1000
algo = 5
"""
Number of Refinement Levels for the Multirate Method
"""
mrMl=1000 #SingleRate: set mL = 1
# Physical parameters for bathymetry smoothing
bNu=20 # Diffusion
bT=10000 # Time for smoothing
#================
#================
# 8. PHYSICAL PARAMETERS FOR GBR SIMULATION
#----------------
"""
-gbrDiff= 0: Constant Diffusion diff || 1: Parametrization of Smagorinsky
"""
gbrDiff=1
diff=1
#================
#================
# 9. SET FORCINGS
#----------------
"""
-gbrWind= 0: NO wind || 1: netCDF data
"""
gbrWind=1
"""
gbrTide= 0: NO tide || 1: tpxo data
"""
gbrTide=1
"""
gbrResidual= 0: NO residual current || 1: standard NCJ forcing
"""
gbrResidual=1
#Define residual current fluxes (North, Central, South)
F_A = 3.0e5
F_B = -10.0e5
F_C = 7.0e5
#Define current direction vectors at boundary in stereographic coordinates
vx_A_stereo = 0.76
vy_A_stereo = -sqrt(1 - vx_A_stereo*vx_A_stereo)
vx_B_stereo = 0.06
vy_B_stereo = -sqrt(1 - vx_B_stereo*vx_B_stereo)
vx_C_stereo = -0.98
vy_C_stereo = -sqrt(1 - vx_C_stereo*vx_C_stereo)
#Define current direction vectors at boundary
vx_A = 0.0
vy_A = sqrt(1 - vx_A*vx_A)
vx_B = 0.76
vy_B = sqrt(1 - vx_B*vx_B)
vy_C = -0.7
vx_C = sqrt(1 - vy_C*vy_C)
#================
#================
#10. Define lat/lon points at which to evaluate current and convert to 2D, projected Cartesian
#----------------
evalPointsLatLon = [
[-14.7406, 145.4253, 'Lizard Island', 0],
[-19.6253, 147.9142, 'Cape Upstart', 1],
[-19.4071, 148.0197, 'Old Reef', 2],
[-18.8311, 148.2896, 'near shelf break', 3],
#
[-19.0600, 147.9597, 'Bowden Reef', 4],
[-19.9826, 148.5833, 'Rattray Island', 5],
[-19.9400, 149.1100, 'Hook', 6],
[-20.8900, 150.1600, 'Bushy', 7],
[-21.8200, 151.1400, 'Bell', 8] ]
[-18.241667, 147.366667, 'Myrmidon Reef', 9], #NB: WAS -18.2452, 147.4100
tidalPointsAbbotLatLon = [
[-19.8187007904, 148.013000488, 'Deepwater West', 0],
[-19.9050006866, 148.005004883, 'Coastal West', 1],
[-19.8815994263, 148.074996948, 'Coastal Footprint', 2],
[-19.8484992981, 148.04699707, 'Inshore Footprint', 3],
[-19.8791999817, 147.970993042, 'Inshore West', 4],
[-19.831199646, 148.115005493, 'Deepwater East', 5],
[-19.9340991974, 148.130004883, 'Coastal East', 6],
[-19.8999004364, 148.158996582, 'Inshore East', 7] ]
tidalPointsSeafarerLatLon = [
[-18.466667, 146.866667, 'Rib Reef', 0],
[-19.238889, 146.838889, 'Townsville', 1],
[-18.933333, 148.066667, 'Shrimp Reef', 2],
[-19.283333, 148.066667, 'Stanley Reef', 3],
[-19.85, 148.116667, 'Abbot Point', 4],
[-19.5, 149.143333, 'Kennedy Reef', 5],
[-20.266667,149.066667, 'Haslewood Island', 6],
[-20.083333, 150.3, 'Bugatti Reef', 7],
[-20.533333, 150.366667, 'Creal Reef', 8],
[-20.997222, 149.9, 'Penrith Island', 9] ]
#Points on reefs
reefPointsLatLon = [
[-19.701517,149.380458, 'Reef A', 0],
[-19.004364,148.100131, 'Reef B', 1],
[-19.419703,148.760181, 'Reef C', 2],
[-19.140206,147.636978, 'Reef D', 3],
[-18.587256,147.100464, 'Passage between reefs', 4],
[-19.921219,149.015717, 'Open sea', 5] ]
#================
This diff is collapsed.
/*
* General library of functions for GBR simulation.
* Part 1: "Numerical & computational" functions
* Part 2: "Physical" functions
*/
#include "fullMatrix.h"
#include "function.h"
#include "dgGroupOfElements.h"
#include "math.h"
#include "polynomialBasis.h"
#include "dgMeshJacobian.h"
static double pOx, pOy, pOz, pPhiX, pPhiY, pPhiZ, pThetaX, pThetaY, pThetaZ, R;
static bool init = false;
static void initialize() {
double phi = 146.5892/180*M_PI;
double theta = -17.147569/180*M_PI;
R = 6.371e6;
pOx = cos(theta)*cos(phi)*R;
pOy = cos(theta)*sin(phi)*R;
pOz = sin(theta)*R;
pPhiX = -sin(phi);
pPhiY = cos(phi);
pPhiZ = 0;
pThetaX = -sin(theta)*cos(phi);
pThetaY = -sin(theta)*sin(phi);
pThetaZ = cos(theta);
init = true;
}
extern "C" {
//Numerical & computational functions:
// void lonLatToStereo (dataCacheMap *, fullMatrix<double> &lonlat, fullMatrix<double> &stereo) {
// {
// double pi;
// pi = acos(-1.0);
//
// //geo coords --> S.P.C.
// double pointTheta = 90.0 - geoLat;
// double pointPhi = geoLon;
//
// //S.P.C. --> Cartesian(on sphere)
// pointCart[0] = R * sin(pointTheta*pi/180.0) * cos(pointPhi*pi/180.0);
// pointCart[1] = R * sin(pointTheta*pi/180.0) * sin(pointPhi*pi/180.0);
// pointCart[2] = R * cos(pointTheta*pi/180.0);
//
// //Cartesian(sphere) --> Stereographic
// pointStereo[0] = (2 * R * pointCart[0]) / (R + pointCart[2]);
// pointStereo[1] = (2 * R * pointCart[1]) / (R + pointCart[2]);
// pointStereo[2] = 0;
// }
void stereoToLonLatDegrees(dataCacheMap *, fullMatrix<double> &lonLat, fullMatrix<double> &XYZ){
if (! init) initialize();
for (size_t i = 0; i < lonLat.size1(); i++) {
double xi = XYZ(i,0);
double zeta = XYZ(i,1);
double x = 4*R*R*xi/(4*R*R +xi*xi + zeta*zeta);
double y = 4*R*R*zeta/(4*R*R +xi*xi + zeta*zeta);
double z = R *(4*R*R-xi*xi-zeta*zeta)/(4*R*R+xi*xi + zeta*zeta);
double lon = atan2(y,x)*180/M_PI;
double lat = asin(z/R)*180/M_PI;
lonLat.set(i,0,lon);
lonLat.set(i,1,lat);
lonLat.set(i,2,0.);
}
}
void stereoToLatLonDegrees(dataCacheMap *, fullMatrix<double> &lonLat, fullMatrix<double> &XYZ){
if (! init) initialize();
for (size_t i = 0; i < lonLat.size1(); i++) {
double xi = XYZ(i,0);
double zeta = XYZ(i,1);
double x = 4*R*R*xi/(4*R*R +xi*xi + zeta*zeta);
double y = 4*R*R*zeta/(4*R*R +xi*xi + zeta*zeta);
double z = R *(4*R*R-xi*xi-zeta*zeta)/(4*R*R+xi*xi + zeta*zeta);
double lon = atan2(y,x)*180/M_PI;
double lat = asin(z/R)*180/M_PI;
lonLat.set(i,1,lon);
lonLat.set(i,0,lat);
lonLat.set(i,2,0.);
}
}
void stereoToLonLat(dataCacheMap *, fullMatrix<double> &lonLat, fullMatrix<double> &XYZ){
if (! init) initialize();
for (size_t i = 0; i < lonLat.size1(); i++) {
double xi = XYZ(i,0);
double zeta = XYZ(i,1);
double x = 4*R*R*xi/(4*R*R +xi*xi + zeta*zeta);
double y = 4*R*R*zeta/(4*R*R +xi*xi + zeta*zeta);
double z = R *(4*R*R-xi*xi-zeta*zeta)/(4*R*R+xi*xi + zeta*zeta);
double lon = atan2(y,x);
double lat = asin(z/R);
lonLat.set(i,0,lon);
lonLat.set(i,1,lat);
lonLat.set(i,2,0.);
}
}
void lonLat (dataCacheMap *, fullMatrix<double> &lonlat, fullMatrix<double> &xyz) {
if (! init) initialize();
for (size_t i = 0; i< lonlat.size1(); i++) {
double x = pPhiX * xyz(i,0) + pThetaX * xyz(i,1) + pOx;
double y = pPhiY * xyz(i,0) + pThetaY * xyz(i,1) + pOy;
double z = pPhiZ * xyz(i,0) + pThetaZ * xyz(i,1) + pOz;
double r = sqrt(x*x+y*y+z*z);
lonlat.set(i, 0, atan2(y,x));
lonlat.set(i, 1, asin(z/r));
lonlat.set(i, 2, 0);
}
}
void latLon (dataCacheMap *, fullMatrix<double> &lonlat, fullMatrix<double> &xyz) {
if (! init) initialize();
for (size_t i = 0; i< lonlat.size1(); i++) {
double x = pPhiX * xyz(i,0) + pThetaX * xyz(i,1) + pOx;
double y = pPhiY * xyz(i,0) + pThetaY * xyz(i,1) + pOy;
double z = pPhiZ * xyz(i,0) + pThetaZ * xyz(i,1) + pOz;
double r = sqrt(x*x+y*y+z*z);
lonlat.set(i, 0, asin(z/r));
lonlat.set(i, 1, atan2(y,x));
lonlat.set(i, 2, 0);
}
}
void lonLatDegrees (dataCacheMap *, fullMatrix<double> &lonlat, fullMatrix<double> &xyz) {
if (! init) initialize();
for (size_t i = 0; i< lonlat.size1(); i++) {
double x = pPhiX * xyz(i,0) + pThetaX * xyz(i,1) + pOx;
double y = pPhiY * xyz(i,0) + pThetaY * xyz(i,1) + pOy;
double z = pPhiZ * xyz(i,0) + pThetaZ * xyz(i,1) + pOz;
double r = sqrt(x*x+y*y+z*z);
lonlat.set(i, 0, atan2(y,x)*180/M_PI);
lonlat.set(i, 1, asin(z/r)*180/M_PI);
lonlat.set(i, 2, 0);
}
}
void latLonDegrees (dataCacheMap *, fullMatrix<double> &lonlat, fullMatrix<double> &xyz) {
if (! init) initialize();
for (size_t i = 0; i< lonlat.size1(); i++) {
double x = pPhiX * xyz(i,0) + pThetaX * xyz(i,1) + pOx;
double y = pPhiY * xyz(i,0) + pThetaY * xyz(i,1) + pOy;
double z = pPhiZ * xyz(i,0) + pThetaZ * xyz(i,1) + pOz;
double r = sqrt(x*x+y*y+z*z);
lonlat.set(i, 1, atan2(y,x)*180/M_PI);
lonlat.set(i, 0, asin(z/r)*180/M_PI);
lonlat.set(i, 2, 0);
}
}
void lonLatVector(dataCacheMap *, fullMatrix<double> &uv, fullMatrix<double> &lonlat, fullMatrix<double> &u, fullMatrix<double> &v){
for (size_t i = 0; i< lonlat.size1(); i++) {
double lon=lonlat(i,0);
double lat=lonlat(i,1);
double uLon=u(i,0);
double uLat=v(i,0);
double lon0=146.5892*M_PI/180;
double lat0=-17.147569*M_PI/180;
double u1= -uLat*sin(lat)*cos(lon)-uLon*sin(lon);
double u2= -uLat*sin(lat)*sin(lon)+uLon*cos(lon);
double u3= uLat*cos(lat);
double alpha1= -sin(lon0);
double alpha2= cos(lon0);
double alpha3=0;
double beta1=-sin(lat0)*cos(lon0);
double beta2=-sin(lat0)*sin(lon0);
double beta3=cos(lat0);
uv.set(i, 0, u1*alpha1+u2*alpha2+u3*alpha3);
uv.set(i, 1, u1*beta1+u2*beta2+u3*beta3);
uv.set(i, 2, 0);
}
}
void latLonVector(dataCacheMap *, fullMatrix<double> &uv, fullMatrix<double> &latlon, fullMatrix<double> &u, fullMatrix<double> &v){
for (size_t i = 0; i< latlon.size1(); i++) {
double lat=latlon(i,0);
double lon=latlon(i,1);
double uLat=u(i,0);
double uLon=v(i,0);
double lon0=146.5892*M_PI/180;
double lat0=-17.147569*M_PI/180;
double u1= -uLat*sin(lat)*cos(lon)-uLon*sin(lon);
double u2= -uLat*sin(lat)*sin(lon)+uLon*cos(lon);
double u3= uLat*cos(lat);
double alpha1= -sin(lon0);
double alpha2= cos(lon0);
double alpha3=0;
double beta1=-sin(lat0)*cos(lon0);
double beta2=-sin(lat0)*sin(lon0);
double beta3=cos(lat0);
uv.set(i, 0, u1*alpha1+u2*alpha2+u3*alpha3);
uv.set(i, 1, u1*beta1+u2*beta2+u3*beta3);
uv.set(i, 2, 0);
}
}
void transport2velocity(dataCacheMap *,fullMatrix<double> &uv, fullMatrix<double> &UV, fullMatrix<double> &bath, fullMatrix<double> &solution){
for (size_t i = 0; i< solution.size1(); i++) {
uv.set(i,0,UV(i,0)/(bath(i,0)+solution(i,0)));
uv.set(i,1,UV(i,1)/(bath(i,0)+solution(i,0)));
uv.set(i,2,UV(i,2)/(bath(i,0)+solution(i,0)));
}
}
void computeEl(dataCacheMap *,fullMatrix<double> &val,fullMatrix<double> &solution) {
for (size_t i = 0; i< val.size1(); i++) {
val.set(i, 0, solution(i,0));
}
}
void computeVel(dataCacheMap *, fullMatrix<double> &val, fullMatrix<double> &solution){
for (size_t i = 0; i< val.size1(); i++) {
val.set(i, 0, solution(i,1));
val.set(i, 1, solution(i,2));
}
}
void computeErrorEl(dataCacheMap *,fullMatrix<double> &err,fullMatrix<double> &solution1, fullMatrix<double> &solution2) {
for (size_t i = 0; i< err.size1(); i++) {
err.set(i, 0, solution1(i,0)-solution2(i, 0));
}
}
void computeErrorVel(dataCacheMap *, fullMatrix<double> &err, fullMatrix<double> &solution1, fullMatrix<double> &solution2){
for (size_t i = 0; i< err.size1(); i++) {
err.set(i, 0, solution1(i,1)-solution2(i, 1));
err.set(i, 1, solution1(i,2)-solution2(i, 2));
}
}
void merge(dataCacheMap *,fullMatrix<double> &sol,fullMatrix<double> &eta, fullMatrix<double> &uv){
for (size_t i = 0; i< sol.size1(); i++) {
sol.set