diff --git a/2d/waves_vegetation/ProcessGauges.ipynb b/2d/waves_vegetation/ProcessGauges.ipynb new file mode 100644 index 00000000..230154a2 --- /dev/null +++ b/2d/waves_vegetation/ProcessGauges.ipynb @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:34173c7a91125f2f0d8f64a932453e1965a134a6c26cc8a84c1c9b866c5d41b5 +size 204330 diff --git a/2d/waves_vegetation/README.rst b/2d/waves_vegetation/README.rst new file mode 100644 index 00000000..f85367ca --- /dev/null +++ b/2d/waves_vegetation/README.rst @@ -0,0 +1,21 @@ +Flow Through Vegetation (Anderson and Smith, 2014) +================================================== + +The problem's domain is a 63.4m long, 1.5m wide flume of varying depth. A +wave is generated at the far edge of the flume, which is at a depth of +1.95m. +The deep section is 5.4m-long, followed by a 1:44 slope for 19.5m +that connects to a 12.2-m long, 1.5m-deep flat testing area. +Within the testing area there are rods simulating idealized vegetation placed in a diamond formation. +After the testing area there is a 11.3m, 1:20 slope with damping properties (to reduce reflected waves). + +.. figure:: ./image_here.bmp + :width: 100% + :align: center + +This study uses PROTEUS to model wave attenuation in shallow-water vegetated areas. + +References +---------- + +- Anderson, M.E. and Smith, J.M. (2013), Wave attenuation by flexible, idealized salt marsh vegetation. Coastal Engineering, Volume 83 p.82-92. diff --git a/2d/waves_vegetation/lonestar-runs/ls_consrv_n.py b/2d/waves_vegetation/lonestar-runs/ls_consrv_n.py new file mode 100644 index 00000000..498ecfd6 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/ls_consrv_n.py @@ -0,0 +1,48 @@ +from proteus import * +from tank import * +from ls_consrv_p import * + +timeIntegrator = ForwardIntegrator +timeIntegration = NoIntegration + +femSpaces = {0:basis} + +subgridError = None +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +shockCapturing = None + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'mcorr_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*mcorr_nl_atol_res +nl_atol_res = mcorr_nl_atol_res +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 diff --git a/2d/waves_vegetation/lonestar-runs/ls_consrv_p.py b/2d/waves_vegetation/lonestar-runs/ls_consrv_p.py new file mode 100644 index 00000000..84f85ac9 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/ls_consrv_p.py @@ -0,0 +1,26 @@ +from proteus import * +from proteus.default_p import * +from tank import * +from proteus.mprans import MCorr + +LevelModelType = MCorr.LevelModel + +coefficients = MCorr.Coefficients(LSModel_index=2,V_model=0,me_model=4,VOFModel_index=1, + applyCorrection=applyCorrection,nd=nd,checkMass=False,useMetrics=useMetrics, + epsFactHeaviside=epsFact_consrv_heaviside, + epsFactDirac=epsFact_consrv_dirac, + epsFactDiffusion=epsFact_consrv_diffusion) + +class zero_phi: + def __init__(self): + pass + def uOfX(self,X): + return 0.0 + def uOfXT(self,X,t): + return 0.0 + +initialConditions = {0:zero_phi()} + + + + diff --git a/2d/waves_vegetation/lonestar-runs/ls_n.py b/2d/waves_vegetation/lonestar-runs/ls_n.py new file mode 100644 index 00000000..4cfc1db3 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/ls_n.py @@ -0,0 +1,62 @@ +from proteus import * +from ls_p import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*ls_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +conservativeFlux = None +numericalFluxType = NCLS.NumericalFlux +subgridError = NCLS.SubgridError(coefficients,nd) +shockCapturing = NCLS.ShockCapturing(coefficients,nd,shockCapturingFactor=ls_shockCapturingFactor,lag=ls_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'ncls_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = ls_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*ls_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + diff --git a/2d/waves_vegetation/lonestar-runs/ls_p.py b/2d/waves_vegetation/lonestar-runs/ls_p.py new file mode 100644 index 00000000..fad3ff3a --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/ls_p.py @@ -0,0 +1,29 @@ +from proteus import * +from proteus.default_p import * +from tank import * +from proteus.mprans import NCLS + +LevelModelType = NCLS.LevelModel + +coefficients = NCLS.Coefficients(V_model=0,RD_model=3,ME_model=2, + checkMass=False, useMetrics=useMetrics, + epsFact=epsFact_consrv_heaviside,sc_uref=ls_sc_uref,sc_beta=ls_sc_beta,movingDomain=movingDomain) + +def getDBC_ls(x,flag): + if flag == boundaryTags['left']: + return wavePhi +# elif flag == boundaryTags['right']: +# return outflowPhi + else: + return None + +dirichletConditions = {0:getDBC_ls} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/2d/waves_vegetation/lonestar-runs/makeRuns.py b/2d/waves_vegetation/lonestar-runs/makeRuns.py new file mode 100644 index 00000000..bdac58af --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/makeRuns.py @@ -0,0 +1,22 @@ +import numpy as np +import os +import shutil + +params = np.loadtxt("tank_parameters.csv", delimiter=',', skiprows=1) +os.system("mkdir runs") +for i in range(params.shape[0]): + filepath = "runs/"+ "run" + `i` + os.system("mkdir " + filepath) + os.system("cp *.py " + filepath) + os.system("cp tank.stampede.slurm " + filepath) + f = open(filepath + "/context.options",'w') + f.write("parallel=True ") + f.write("wave_type='single-peaked' ") + f.write("gauges=True ") + f.write("depth=" + `params[i][0]` + " ") + f.write("wave_height=" + `params[i][1]` + " ") + f.write("peak_period=" + `params[i][2]` + " ") + f.write("peak_wavelength=" + `params[i][3]` + " ") + f.write("tank_height=" + `params[i][4]`) + f.close() + diff --git a/2d/waves_vegetation/lonestar-runs/petsc.options.asm b/2d/waves_vegetation/lonestar-runs/petsc.options.asm new file mode 100644 index 00000000..3c331794 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/petsc.options.asm @@ -0,0 +1,7 @@ +-rans2p_ksp_type gmres -rans2p_pc_type asm -rans2p_pc_asm_type basic -rans2p_ksp_max_it 2000 +-rans2p_ksp_gmres_modifiedgramschmidt -rans2p_ksp_gmres_restart 300 -rans2p_sub_ksp_type preonly -rans2p_sub_pc_factor_mat_solver_package superlu -rans2p_ksp_knoll -rans2p_sub_pc_type lu +-ncls_ksp_type gmres -ncls_pc_type hypre -ncls_pc_hypre_type boomeramg -ncls_ksp_gmres_restart 300 -ncls_ksp_knoll -ncls_ksp_max_it 2000 +-vof_ksp_type gmres -vof_pc_type hypre -vof_pc_hypre_type boomeramg -vof_ksp_gmres_restart 300 -vof_ksp_knoll -vof_ksp_max_it 2000 +-rdls_ksp_type gmres -rdls_pc_type asm -rdls_pc_asm_type basic -rdls_ksp_gmres_modifiedgramschmidt -rdls_ksp_gmres_restart 300 -rdls_ksp_knoll -rdls_sub_ksp_type preonly -rdls_sub_pc_factor_mat_solver_package superlu -rdls_sub_pc_type lu -rdls_ksp_max_it 2000 +-mcorr_ksp_type cg -mcorr_pc_type hypre -mcorr_pc_hypre_type boomeramg -mcorr_ksp_max_it 2000 +-log_summary diff --git a/2d/waves_vegetation/lonestar-runs/redist_n.py b/2d/waves_vegetation/lonestar-runs/redist_n.py new file mode 100644 index 00000000..237f5e17 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/redist_n.py @@ -0,0 +1,67 @@ +from proteus import * +from redist_p import * +from tank import * + +nl_atol_res = rd_nl_atol_res +tolFac = 0.0 +nl_atol_res = rd_nl_atol_res + +linTolFac = 0.01 +l_atol_res = 0.01*rd_nl_atol_res + +if redist_Newton: + timeIntegration = NoIntegration + stepController = Newton_controller + maxNonlinearIts = 50 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'r' + levelNonlinearSolverConvergenceTest = 'r' + linearSolverConvergenceTest = 'r-true' + useEisenstatWalker = False +else: + timeIntegration = BackwardEuler_cfl + stepController = RDLS.PsiTC + runCFL=2.0 + psitc['nStepsForce']=3 + psitc['nStepsMax']=50 + psitc['reduceRatio']=2.0 + psitc['startRatio']=1.0 + rtol_res[0] = 0.0 + atol_res[0] = rd_nl_atol_res + useEisenstatWalker = False + maxNonlinearIts = 1 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'rits' + levelNonlinearSolverConvergenceTest = 'rits' + linearSolverConvergenceTest = 'r-true' + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +subgridError = RDLS.SubgridError(coefficients,nd) +shockCapturing = RDLS.ShockCapturing(coefficients,nd,shockCapturingFactor=rd_shockCapturingFactor,lag=rd_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = NLGaussSeidel +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rdls_' + diff --git a/2d/waves_vegetation/lonestar-runs/redist_p.py b/2d/waves_vegetation/lonestar-runs/redist_p.py new file mode 100644 index 00000000..5700202e --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/redist_p.py @@ -0,0 +1,32 @@ +from proteus import * +from proteus.default_p import * +from math import * +from tank import * +from proteus.mprans import RDLS +""" +The redistancing equation in the sloshbox test problem. +""" + +LevelModelType = RDLS.LevelModel + +coefficients = RDLS.Coefficients(applyRedistancing=applyRedistancing, + epsFact=epsFact_redistance, + nModelId=2, + rdModelId=3, + useMetrics=useMetrics, + backgroundDiffusionFactor=backgroundDiffusionFactor) + +def getDBC_rd(x,flag): + pass + +dirichletConditions = {0:getDBC_rd} +weakDirichletConditions = {0:RDLS.setZeroLSweakDirichletBCsSimple} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/2d/waves_vegetation/lonestar-runs/tank.py b/2d/waves_vegetation/lonestar-runs/tank.py new file mode 100644 index 00000000..7ecfe840 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/tank.py @@ -0,0 +1,584 @@ +from math import * +import proteus.MeshTools +from proteus import Domain +from proteus.default_n import * +from proteus.Profiling import logEvent +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.ctransportCoefficients import smoothedHeaviside_integral +from proteus import Gauges +from proteus.Gauges import PointGauges,LineGauges,LineIntegralGauges +from proteus import WaveTools as WT + +from proteus import Context +opts=Context.Options([ + ("wave_type", 'linear', "type of waves generated: 'linear', 'Nonlinear', 'single-peaked', 'double-peaked'"), + ("depth", 0.457, "water depth at leading edge of vegetation (not at wave generator)[m]"), + ("wave_height", 0.192, "wave height at leading edget of vegetation [m]"), + ("peak_period", 2.0, "Peak period [s]"), + ("peak_period2", 1.5, "Second peak period (only used in double-peaked case)[s]"), + ("peak_wavelength",3.91,"Peak wavelength in [m]"), + ("parallel", False, "Run in parallel"), + ("gauges", True, "Enable gauges"), + ("tank_height", 1.0, "height of wave tank")]) + +#wave generator +windVelocity = (0.0,0.0) +veg_platform_height = 17.2/44.0 + 6.1/20.0 +depth = veg_platform_height + opts.depth +inflowHeightMean = depth +inflowVelocityMean = (0.0,0.0) +period = opts.peak_period +omega = 2.0*math.pi/period +waveheight = opts.wave_height +amplitude = waveheight/ 2.0 +wavelength = opts.peak_wavelength +k = 2.0*math.pi/wavelength +waveDir = numpy.array([1,0,0]) +g = numpy.array([0,-9.81,0]) +if opts.wave_type == 'linear': + waves = WT.MonochromaticWaves(period = period, # Peak period + waveHeight = waveheight, # Height + depth = depth, # Depth + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + waveType="Linear") +elif opts.wave_type == 'Nonlinear': + waves = WT.MonochromaticWaves(period = period, # Peak period + waveHeight = waveheight, # Height + wavelength = wavelength, + depth = depth, # Depth + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + waveType="Fenton", + Ycoeff = [0.04160592, #Surface elevation Fourier coefficients for non-dimensionalised solution + 0.00555874, + 0.00065892, + 0.00008144, + 0.00001078, + 0.00000151, + 0.00000023, + 0.00000007], + Bcoeff = [0.05395079, + 0.00357780, + 0.00020506, + 0.00000719, + -0.00000016, + -0.00000005, + 0.00000000, + 0.00000000]) +elif opts.wave_type == 'single-peaked': + waves = WT.RandomWaves( Tp = period, # Peak period + Hs = waveheight, # Height + mwl = inflowHeightMean, # Sea water level + depth = depth, # Depth + #fp = 1./period, #peak Frequency + bandFactor = 2.0, #fmin=fp/Bandfactor, fmax = Bandfactor * fp + N = 101, #No of frequencies for signal reconstruction + #mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, + spectName = "JONSWAP") # Gravity vector, defines the vertical + #gamma=3.3, + #spec_fun = JONSWAP) +elif opts.wave_type == 'double-peaked': + waves = WT.DoublePeakedRandomWaves( #Tp = period, # Peak period + Hs = waveheight, # Height + d = depth, # Depth + fp = 1./period, #peak Frequency + bandFactor = 2.0, #fmin=fp/Bandfactor, fmax = Bandfactor * fp + N = 101, #No of frequencies for signal reconstruction + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + gamma=10.0, + spec_fun = WT.JONSWAP, + Tp_2 = opts.peak_period2) + +gauges=opts.gauges +# Discretization -- input options +genMesh=True +movingDomain=False +applyRedistancing=True +useOldPETSc=False +useSuperlu=not opts.parallel +timeDiscretization='be'#'be','vbdf','flcbdf' +spaceOrder = 1 +useHex = False +useRBLES = 0.0 +useMetrics = 1.0 +applyCorrection=True +useVF = 1.0 +useOnlyVF = False +useRANS = 0 # 0 -- None + # 1 -- K-Epsilon + # 2 -- K-Omega +# Input checks +if spaceOrder not in [1,2]: + print "INVALID: spaceOrder" + spaceOrder + sys.exit() + +if useRBLES not in [0.0, 1.0]: + print "INVALID: useRBLES" + useRBLES + sys.exit() + +if useMetrics not in [0.0, 1.0]: + print "INVALID: useMetrics" + sys.exit() + +# Discretization +nd = 2 +if spaceOrder == 1: + hFactor=1.0 + if useHex: + basis=C0_AffineLinearOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,2) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,2) + else: + basis=C0_AffineLinearOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,3) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,3) +elif spaceOrder == 2: + hFactor=0.5 + if useHex: + basis=C0_AffineLagrangeOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,4) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,4) + else: + basis=C0_AffineQuadraticOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,4) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,4) + +# Domain and mesh + +#for debugging, make the tank short +L = (45.4,opts.tank_height) +he = 0.025 #float(wavelength)/130.0 #100.0 #50.0#0.0#100 + +GenerationZoneLength = wavelength +AbsorptionZoneLength= 45.4-37.9 +spongeLayer = True +xSponge = GenerationZoneLength +xRelaxCenter = xSponge/2.0 +epsFact_solid = xSponge/2.0 +#zone 2 +xSponge_2 = 37.9 +xRelaxCenter_2 = 0.5*(37.9+45.4) +epsFact_solid_2 = AbsorptionZoneLength/2.0 + +weak_bc_penalty_constant = 100.0 +nLevels = 1 +#parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.element +parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.node +nLayersOfOverlapForParallel = 0 +structured=False + + +gauge_dx=0.075 +PGL=[] +gauge_top = 19.5/44.0 + 1.5 +veg_gauge_bottom_y = 17.2/44.0 + 6.1/20.0 +LGL=[[(6.1, (6.1 - 5.4)/44.0, 0.0), (6.1,gauge_top,0.0)],#1 Goda + [(6.4, (6.4 - 5.4)/44.0, 0.0), (6.4,gauge_top,0.0)],#2 Goda + [(7.0, (7.0 - 5.4)/44.0, 0.0), (7.0,gauge_top,0.0)],#3 Goda + [(26.0, veg_gauge_bottom_y, 0.0), (26.0,gauge_top,0.0)],#4 veg + [(26.9, veg_gauge_bottom_y, 0.0), (26.9,gauge_top,0.0)],#5 + [(27.4, veg_gauge_bottom_y, 0.0), (27.4,gauge_top,0.0)],#6 + [(27.9, veg_gauge_bottom_y, 0.0), (27.9,gauge_top,0.0)],#7 + [(28.5, veg_gauge_bottom_y, 0.0), (28.5,gauge_top,0.0)],#8 + [(29.5, veg_gauge_bottom_y, 0.0), (29.5,gauge_top,0.0)],#9 + [(31.0, veg_gauge_bottom_y, 0.0), (31.0,gauge_top,0.0)],#10 + [(32.7, veg_gauge_bottom_y, 0.0), (32.7,gauge_top,0.0)],#11 + [(34.4, veg_gauge_bottom_y, 0.0), (34.4,gauge_top,0.0)],#12 + [(36.2, veg_gauge_bottom_y, 0.0), (36.2,gauge_top,0.0)]]#13 +#for i in range(0,int(L[0]/gauge_dx+1)): #+1 only if gauge_dx is an exact +# PGL.append([gauge_dx*i,0.5,0]) +# LGL.append([(gauge_dx*i,0.0,0),(gauge_dx*i,L[1],0)]) + + +gaugeLocations=tuple(map(tuple,PGL)) +columnLines=tuple(map(tuple,LGL)) + + +#pointGauges = PointGauges(gauges=((('u','v'), gaugeLocations), +# (('p',), gaugeLocations)), +# activeTime = (0, 1000.0), +# sampleRate = 0, +# fileName = 'combined_gauge_0_0.5_sample_all.txt') + +#print gaugeLocations +#print columnLines + +fields = ('vof',) + +columnGauge = LineIntegralGauges(gauges=((fields, columnLines),), + fileName='column_gauge.csv') + +#v_resolution = max(he,0.05) +#linePoints = int((gauge_top - veg_gauge_bottom_y)/v_resolution) +lineGauges = LineGauges(gauges=((('u','v'),#fields in gauge set + (#lines for these fields + ((26.0, veg_gauge_bottom_y, 0.0),(26.0, gauge_top, 0.0)), + ),#end lines + ),#end gauge set + ),#end gauges + fileName="vegZoneVelocity.csv") + +#lineGauges_phi = LineGauges_phi(lineGauges.endpoints,linePoints=20) + + +if useHex: + nnx=ceil(L[0]/he)+1 + nny=ceil(L[1]/he)+1 + hex=True + domain = Domain.RectangularDomain(L) +else: + boundaries=['left','right','bottom','top','front','back'] + boundaryTags=dict([(key,i+1) for (i,key) in enumerate(boundaries)]) + if structured: + nnx=ceil(L[0]/he)+1 + nny=ceil(L[1]/he)+1 + elif spongeLayer: + #tp = opts.tank_height + bp = 20.0*(opts.tank_height-0.2) + vertices=[[0.0, 0.0 ],#0 begin wave paddle bottom + [5.4, 0.0 ],#1 end wave paddle bottom, begin incline 1 + [5.4 + 17.2, 17.2/44.0 ],#2 end incline 1, begin incline 2 + [5.4 + 17.2 + 6.1, 17.2/44.0 + 6.1/20.0 ],#3 end incline 2, begin pre veg platform + [5.4 + 17.2 + 6.1 + 1.2, 17.2/44.0 + 6.1/20.0 ],#4 end pre veg platform, begin veg zone + [5.4 + 17.2 + 6.1 + 1.2 + 9.8, 17.2/44.0 + 6.1/20.0 ],#5 end veg zone, begin post veg platform + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2, 17.2/44.0 + 6.1/20.0 ],#6 -- sponge, end post veg platorm, begin slope bottom + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2 + bp, 17.2/44.0 + 6.1/20.0 + bp/20.0],#7 end slope bottom + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2 + bp, 17.2/44.0 + 6.1/20.0 + bp/20.0+0.2],#8 end slope top + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2, 19.5/44.0 + 1.5],#9 -- sponge begin sponge top + [0.0, 19.5/44.0 + 1.5]]#10 begin wave paddle top + + vertexFlags=[boundaryTags['bottom'],#0 + boundaryTags['bottom'],#1 + boundaryTags['bottom'],#2 + boundaryTags['bottom'],#3 + boundaryTags['bottom'],#4 + boundaryTags['bottom'],#5 + boundaryTags['bottom'],#6 + boundaryTags['bottom'],#7 + boundaryTags['top'],#8 + boundaryTags['top'],#9 + boundaryTags['top']]#10 + segments=[[0,1],#0 + [1,2],#1 + [2,3],#2 + [3,4],#3 + [4,5],#4 + [5,6],#5 + [6,7],#6 + [7,8],#7 + [8,9],#8 + [9,10],#8 + [10,0],#9 + [6,9]]#10 + + segmentFlags=[boundaryTags['bottom'],#0 + boundaryTags['bottom'],#1 + boundaryTags['bottom'],#2 + boundaryTags['bottom'],#3 + boundaryTags['bottom'],#4 + boundaryTags['bottom'],#5 + boundaryTags['bottom'],#6 + boundaryTags['right'],#7 + boundaryTags['top'],#8 + boundaryTags['top'],#9 + boundaryTags['left'],#10 + 0]#11 + + regions=[[0.5,0.5], + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2+1.0, 17.2/44.0 + 6.1/20.0 + 1.0]] + regionFlags=[1,2] + domain = Domain.PlanarStraightLineGraphDomain(vertices=vertices, + vertexFlags=vertexFlags, + segments=segments, + segmentFlags=segmentFlags, + regions=regions, + regionFlags=regionFlags) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="VApq30Dena%8.8f" % ((he**2)/2.0,) + print triangleOptions + logEvent("""Mesh generated using: triangle -%s %s""" % (triangleOptions,domain.polyfile+".poly")) + porosityTypes = numpy.array([1.0, + 1.0, + 1.0]) + dragAlphaTypes = numpy.array([0.0, + 0.0, + 0.5/1.004e-6]) + dragBetaTypes = numpy.array([0.0,0.0,0.0]) + + epsFact_solidTypes = np.array([0.0,0.0,epsFact_solid_2]) + + else: + vertices=[[0.0,0.0],#0 + [L[0],0.0],#1 + [L[0],L[1]],#2 + [0.0,L[1]]]#3 + + vertexFlags=[boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['top'], + boundaryTags['top']] + segments=[[0,1], + [1,2], + [2,3], + [3,0] + ] + segmentFlags=[boundaryTags['bottom'], + boundaryTags['right'], + boundaryTags['top'], + boundaryTags['left']] + + regions=[ [ 0.1*L[0] , 0.1*L[1] ], + [0.95*L[0] , 0.95*L[1] ] ] + regionFlags=[1,2] + domain = Domain.PlanarStraightLineGraphDomain(vertices=vertices, + vertexFlags=vertexFlags, + segments=segments, + segmentFlags=segmentFlags, + regions=regions, + regionFlags=regionFlags) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="VApq30Dena%8.8f" % ((he**2)/2.0,) + + logEvent("""Mesh generated using: triangle -%s %s""" % (triangleOptions,domain.polyfile+".poly")) +# Time stepping +T=200.0 #480.0 #480.0 #40*period +dt_fixed = period/30.0#2.0*0.5/20.0#T/2.0#period/21.0 +dt_init = min(0.001*dt_fixed,0.001) +runCFL=10.90 +nDTout = int(round(T/dt_fixed)) + +# Numerical parameters +ns_forceStrongDirichlet = False +backgroundDiffusionFactor=0.0 +if useMetrics: + ns_shockCapturingFactor = 0.25 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.25 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.25 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.25 + rd_lag_shockCapturing = False + epsFact_density = 3.0 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 1.5 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.1 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.1 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 +else: + ns_shockCapturingFactor = 0.9 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.9 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.9 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.9 + rd_lag_shockCapturing = False + epsFact_density = 1.5 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 0.33 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.9 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.9 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 + +ns_nl_atol_res = max(1.0e-10,0.001*he**2) +vof_nl_atol_res = max(1.0e-10,0.001*he**2) +ls_nl_atol_res = max(1.0e-10,0.001*he**2) +rd_nl_atol_res = max(1.0e-10,0.05*he) +mcorr_nl_atol_res = max(1.0e-10,0.001*he**2) +kappa_nl_atol_res = max(1.0e-10,0.001*he**2) +dissipation_nl_atol_res = max(1.0e-10,0.001*he**2) + +#turbulence +ns_closure=0 #1-classic smagorinsky, 2-dynamic smagorinsky, 3 -- k-epsilon, 4 -- k-omega +if useRANS == 1: + ns_closure = 3 +elif useRANS == 2: + ns_closure == 4 +# Water +rho_0 = 998.2 +nu_0 = 1.004e-6 + +# Air +rho_1 = 1.205 +nu_1 = 1.500e-5 + +# Surface tension +sigma_01 = 0.0 + +# Gravity +g = [0.0,-9.8] + +# Initial condition +waterLine_x = 2*L[0] +waterLine_z = inflowHeightMean + + +def signedDistance(x): + phi_x = x[0]-waterLine_x + phi_z = x[1]-waterLine_z + if phi_x < 0.0: + if phi_z < 0.0: + return max(phi_x,phi_z) + else: + return phi_z + else: + if phi_z < 0.0: + return phi_x + else: + return sqrt(phi_x**2 + phi_z**2) + + +def theta(x,t): + return k*x[0] - omega*t + pi/2.0 + +def z(x): + return x[1] - inflowHeightMean + +#sigma = omega - k*inflowVelocityMean[0] +h = inflowHeightMean # - transect[0][1] if lower left hand corner is not at z=0 + +#waveData + +def waveHeight(x,t): + return inflowHeightMean + waves.eta(x[0],x[1],x[2],t) +def waveVelocity_u(x,t): + return waves.u(x[0],x[1],x[2],t,"x") +def waveVelocity_v(x,t): + return waves.u(x[0],x[1],x[2],t,"y") +def waveVelocity_w(x,t): + return waves.u(x[0],x[1],x[2],t,"z") + +#solution variables + +def wavePhi(x,t): + return x[1] - waveHeight(x,t) + +def waveVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)) + +def twpflowVelocity_u(x,t): + waterspeed = waveVelocity_u(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + u = H*windVelocity[0] + (1.0-H)*waterspeed + return u + +def twpflowVelocity_v(x,t): + waterspeed = waveVelocity_v(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + return H*windVelocity[1]+(1.0-H)*waterspeed + +def twpflowFlux(x,t): + return -twpflowVelocity_u(x,t) + +outflowHeight=inflowHeightMean + +def outflowVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,x[1] - outflowHeight) + +def outflowPhi(x,t): + return x[1] - outflowHeight + +def outflowPressure(x,t): + if x[1]>inflowHeightMean: + return (L[1]-x[1])*rho_1*abs(g[1]) + else: + return (L[1]-inflowHeightMean)*rho_1*abs(g[1])+(inflowHeightMean-x[1])*rho_0*abs(g[1]) + + + #p_L = L[1]*rho_1*g[1] + #phi_L = L[1] - outflowHeight + #phi = x[1] - outflowHeight + #return p_L -g[1]*(rho_0*(phi_L - phi)+(rho_1 -rho_0)*(smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi_L) + # -smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi))) + +def twpflowVelocity_w(x,t): + return 0.0 + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + +RelaxationZone = namedtuple("RelaxationZone","center_x sign u v w") + +class RelaxationZoneWaveGenerator(AV_base): + """ Prescribe a velocity penalty scaling in a material zone via a Darcy-Forchheimer penalty + + :param zones: A dictionary mapping integer material types to Zones, where a Zone is a named tuple + specifying the x coordinate of the zone center and the velocity components + """ + def __init__(self,zones): + assert isinstance(zones,dict) + self.zones = zones + def calculate(self): + for l,m in enumerate(self.model.levelModelList): + for eN in range(m.coefficients.q_phi.shape[0]): + mType = m.mesh.elementMaterialTypes[eN] + if self.zones.has_key(mType): + for k in range(m.coefficients.q_phi.shape[1]): + t = m.timeIntegration.t + x = m.q['x'][eN,k] + m.coefficients.q_phi_solid[eN,k] = self.zones[mType].sign*(self.zones[mType].center_x - x[0]) + m.coefficients.q_velocity_solid[eN,k,0] = self.zones[mType].u(x,t) + m.coefficients.q_velocity_solid[eN,k,1] = self.zones[mType].v(x,t) + #m.coefficients.q_velocity_solid[eN,k,2] = self.zones[mType].w(x,t) + m.q['phi_solid'] = m.coefficients.q_phi_solid + m.q['velocity_solid'] = m.coefficients.q_velocity_solid + +rzWaveGenerator = RelaxationZoneWaveGenerator(zones={ + # 1:RelaxationZone(xRelaxCenter, + # 1.0, + # twpflowVelocity_u, + # twpflowVelocity_v, + # twpflowVelocity_w), + 2:RelaxationZone(xRelaxCenter_2, + -1.0, #currently Hs=1-exp_function + zeroVel, + zeroVel, + zeroVel)}) diff --git a/2d/waves_vegetation/lonestar-runs/tank.stampede.slurm b/2d/waves_vegetation/lonestar-runs/tank.stampede.slurm new file mode 100644 index 00000000..e90581f5 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/tank.stampede.slurm @@ -0,0 +1,17 @@ +#!/bin/bash +#SBATCH -J wavetank #Job Name +#SBATCH -o oe.wavetank.asm.o%j #Output and Error File +#SBATCH -n 12 #Number of mpi asks +#SBATCH -p normal #Queue +#SBATCH -t 48:00:00 #Run Time +#SBATCH --mail-user=steve.a.mattis@gmail.com #email +#SBATCH --mail-type=begin #when to email +#SBATCH --mail-type=end #when to email +#SBATCH -A ADCIRC #account +set -x +#source ${WORK}/src/proteus-hashstack8/envConfig +#source /scratch/01082/smattis/src/proteus-hashdist-11-2015/envConfig +source /home1/01082/smattis/src/proteus/envConfig +mkdir $SLURM_JOB_NAME.$SLURM_JOB_ID +#/usr/local/bin/ibrun /scratch/01082/smattis/src/proteus-hashdist-11-2015/proteus/stampede.gnu/bin/python /scratch/01082/smattis/src/proteus-hashdist-11-2015/proteus/stampede.gnu/bin/parun tank_so.py -l 3 -v -D $SLURM_JOB_NAME.$SLURM_JOB_ID -O ../../../inputTemplates/petsc.options.asm -o context.options -p +ibrun parun tank_so.py -l 3 -v -D $SLURM_JOB_NAME.$SLURM_JOB_ID -O ../../petsc.options.asm -o context.options #-p \ No newline at end of file diff --git a/2d/waves_vegetation/lonestar-runs/tank_batch.py b/2d/waves_vegetation/lonestar-runs/tank_batch.py new file mode 100644 index 00000000..38fabd1d --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/tank_batch.py @@ -0,0 +1,4 @@ +simFlagsList[0]['storeQuantities']= ["q:'phi_solid'","q:'velocity_solid'"] +#simFlagsList[0]['storeQuantities']= ["q:velocity_solid"] +start +quit diff --git a/2d/waves_vegetation/lonestar-runs/tank_parameters.csv b/2d/waves_vegetation/lonestar-runs/tank_parameters.csv new file mode 100644 index 00000000..657f2e9f --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/tank_parameters.csv @@ -0,0 +1,16 @@ +depth,wave_height,peak_period,peak_wavelength,tank_height +0.533,0.111,1.5,2.89,0.9 +0.533,0.11,1.75,3.53,0.9 +0.533,0.112,2,4.16,0.9 +0.457,0.081,1.5,2.74,0.8 +0.457,0.109,1.5,2.74,0.8 +0.457,0.139,1.5,2.74,0.8 +0.457,0.05,2,3.91,0.8 +0.457,0.107,2,3.91,0.8 +0.457,0.153,2,3.91,0.8 +0.457,0.192,2,3.91,0.8 +0.305,0.113,1.25,1.88,0.8 +0.305,0.11,1.5,2.36,0.8 +0.305,0.112,1.75,2.82,0.8 +0.305,0.111,2,3.28,0.8 +0.305,0.112,2.25,3.73,0.8 diff --git a/2d/waves_vegetation/lonestar-runs/tank_so.py b/2d/waves_vegetation/lonestar-runs/tank_so.py new file mode 100644 index 00000000..d06f6ffa --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/tank_so.py @@ -0,0 +1,35 @@ +from proteus.default_so import * +import tank +from proteus import Context +Context.setFromModule(tank) +ctx = Context.get() + +if tank.useOnlyVF: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n")] +else: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n"), + ("ls_p", "ls_n"), + ("redist_p", "redist_n"), + ("ls_consrv_p", "ls_consrv_n")] + + +if tank.useRANS > 0: + pnList.append(("kappa_p", + "kappa_n")) + pnList.append(("dissipation_p", + "dissipation_n")) +name = "tank_p" + +if tank.timeDiscretization == 'flcbdf': + systemStepControllerType = Sequential_MinFLCBDFModelStep + systemStepControllerType = Sequential_MinAdaptiveModelStep +else: + systemStepControllerType = Sequential_MinAdaptiveModelStep + +needEBQ_GLOBAL = False +needEBQ = False + +tnList = [0.0,tank.dt_init]+[i*tank.dt_fixed for i in range(1,tank.nDTout+1)] +archiveFlag = ArchiveFlags.EVERY_USER_STEP diff --git a/2d/waves_vegetation/lonestar-runs/twp_navier_stokes_n.py b/2d/waves_vegetation/lonestar-runs/twp_navier_stokes_n.py new file mode 100644 index 00000000..49d4d062 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/twp_navier_stokes_n.py @@ -0,0 +1,68 @@ +from proteus import * +from twp_navier_stokes_p import * +from tank import * +from proteus import Context +ctx = Context.get() + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller_sys + stepController = Min_dt_cfl_controller + time_tol = 10.0*ns_nl_atol_res + atol_u = {1:time_tol,2:time_tol} + rtol_u = {1:time_tol,2:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis, + 1:basis, + 2:basis} + +massLumping = False +numericalFluxType = None +conservativeFlux = None + +numericalFluxType = RANS2P.NumericalFlux +subgridError = RANS2P.SubgridError(coefficients,nd,lag=ns_lag_subgridError,hFactor=hFactor,nStepsToDelay=1) +shockCapturing = RANS2P.ShockCapturing(coefficients,nd,ns_shockCapturingFactor,lag=ns_lag_shockCapturing,nStepsToDelay=1) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None + +#linearSmoother = SimpleNavierStokes2D +linearSmoother = None + +matrix = SparseMatrix + +multilevelLinearSolver = KSP_petsc4py +levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rans2p_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*ns_nl_atol_res +nl_atol_res = ns_nl_atol_res +useEisenstatWalker = False +maxNonlinearIts = 100 +maxLineSearches = 0 +conservativeFlux = {0:'pwl-bdm-opt'} + +if ctx.gauges: + auxiliaryVariables=[lineGauges] +#auxiliaryVariables=[pointGauges,rzWaveGenerator] diff --git a/2d/waves_vegetation/lonestar-runs/twp_navier_stokes_p.py b/2d/waves_vegetation/lonestar-runs/twp_navier_stokes_p.py new file mode 100644 index 00000000..10a8a9dc --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/twp_navier_stokes_p.py @@ -0,0 +1,141 @@ +from proteus import * +from proteus.default_p import * +from tank import * +from proteus.mprans import RANS2P + +LevelModelType = RANS2P.LevelModel +if useOnlyVF: + LS_model = None +else: + LS_model = 2 +if useRANS >= 1: + Closure_0_model = 5; Closure_1_model=6 + if useOnlyVF: + Closure_0_model=2; Closure_1_model=3 + if movingDomain: + Closure_0_model += 1; Closure_1_model += 1 +else: + Closure_0_model = None + Closure_1_model = None + +if spongeLayer or levee or slopingSpongeLayer: + coefficients = RANS2P.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain, + porosityTypes=porosityTypes, + dragAlphaTypes=dragAlphaTypes, + dragBetaTypes=dragBetaTypes, + epsFact_solid = epsFact_solidTypes) +else: + coefficients = RANS2P.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain) + +def getDBC_p(x,flag): + if flag == boundaryTags['top']: + return lambda x,t: 0.0 +# elif flag == boundaryTags['right']: +# return outflowPressure + +def getDBC_u(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_u +# elif flag == boundaryTags['right']: +# return lambda x,t: 0.0 + +def getDBC_v(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_v +# elif flag == boundaryTags['right']: +# return lambda x,t: 0.0 + +dirichletConditions = {0:getDBC_p, + 1:getDBC_u, + 2:getDBC_v} + +def getAFBC_p(x,flag): + if flag == boundaryTags['left']: + return twpflowFlux + elif flag == boundaryTags['bottom'] or flag == boundaryTags['right'] or flag == 0: + return lambda x,t: 0.0 + +def getAFBC_u(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['right'] or flag == 0: + return lambda x,t: 0.0 + +def getAFBC_v(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['right'] or flag == 0: + return lambda x,t: 0.0 + +def getDFBC_u(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + +def getDFBC_v(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + +advectiveFluxBoundaryConditions = {0:getAFBC_p, + 1:getAFBC_u, + 2:getAFBC_v} + +diffusiveFluxBoundaryConditions = {0:{}, + 1:{1:getDFBC_u}, + 2:{2:getDFBC_v}} + +class PerturbedSurface_p: + def __init__(self,waterLevel): + self.waterLevel=waterLevel + def uOfXT(self,x,t): + if signedDistance(x) < 0: + return -(L[1] - self.waterLevel)*rho_1*g[1] - (self.waterLevel - x[1])*rho_0*g[1] + else: + return -(L[1] - self.waterLevel)*rho_1*g[1] + +class AtRest: + def __init__(self): + pass + def uOfXT(self,x,t): + return 0.0 + +initialConditions = {0:PerturbedSurface_p(waterLine_z), + 1:AtRest(), + 2:AtRest()} diff --git a/2d/waves_vegetation/lonestar-runs/vof_n.py b/2d/waves_vegetation/lonestar-runs/vof_n.py new file mode 100644 index 00000000..5a15d474 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/vof_n.py @@ -0,0 +1,68 @@ +from proteus import * +from tank import * +from vof_p import * +from proteus import Context +ctx = Context.get() + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*vof_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = VOF.NumericalFlux +conservativeFlux = None +subgridError = VOF.SubgridError(coefficients=coefficients,nd=nd) +shockCapturing = VOF.ShockCapturing(coefficients,nd,shockCapturingFactor=vof_shockCapturingFactor,lag=vof_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'vof_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = vof_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*vof_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + + +if ctx.gauges: + auxiliaryVariables = [columnGauge] diff --git a/2d/waves_vegetation/lonestar-runs/vof_p.py b/2d/waves_vegetation/lonestar-runs/vof_p.py new file mode 100644 index 00000000..94965ec8 --- /dev/null +++ b/2d/waves_vegetation/lonestar-runs/vof_p.py @@ -0,0 +1,46 @@ +from proteus import * +from proteus.default_p import * +from proteus.ctransportCoefficients import smoothedHeaviside +from tank import * +from proteus.mprans import VOF + +LevelModelType = VOF.LevelModel +if useOnlyVF: + RD_model = None + LS_model = None +else: + RD_model = 3 + LS_model = 2 + +coefficients = VOF.Coefficients(LS_model=LS_model,V_model=0,RD_model=RD_model,ME_model=1, + checkMass=False,useMetrics=useMetrics, + epsFact=epsFact_vof,sc_uref=vof_sc_uref,sc_beta=vof_sc_beta,movingDomain=movingDomain) + +def getDBC_vof(x,flag): + if flag == boundaryTags['left']: + return waveVF + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return lambda x,t: 1.0 +# elif flag == boundaryTags['right']: +# return outflowVF + +dirichletConditions = {0:getDBC_vof} + +def getAFBC_vof(x,flag): + if flag == boundaryTags['left']: + return None + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return None + # elif flag == boundaryTags['right']: + # return None + else: + return lambda x,t: 0.0 + +advectiveFluxBoundaryConditions = {0:getAFBC_vof} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_H: + def uOfXT(self,x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,signedDistance(x)) + +initialConditions = {0:PerturbedSurface_H()} diff --git a/2d/waves_vegetation/ls_consrv_n.py b/2d/waves_vegetation/ls_consrv_n.py new file mode 100644 index 00000000..498ecfd6 --- /dev/null +++ b/2d/waves_vegetation/ls_consrv_n.py @@ -0,0 +1,48 @@ +from proteus import * +from tank import * +from ls_consrv_p import * + +timeIntegrator = ForwardIntegrator +timeIntegration = NoIntegration + +femSpaces = {0:basis} + +subgridError = None +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +shockCapturing = None + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'mcorr_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*mcorr_nl_atol_res +nl_atol_res = mcorr_nl_atol_res +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 diff --git a/2d/waves_vegetation/ls_consrv_p.py b/2d/waves_vegetation/ls_consrv_p.py new file mode 100644 index 00000000..84f85ac9 --- /dev/null +++ b/2d/waves_vegetation/ls_consrv_p.py @@ -0,0 +1,26 @@ +from proteus import * +from proteus.default_p import * +from tank import * +from proteus.mprans import MCorr + +LevelModelType = MCorr.LevelModel + +coefficients = MCorr.Coefficients(LSModel_index=2,V_model=0,me_model=4,VOFModel_index=1, + applyCorrection=applyCorrection,nd=nd,checkMass=False,useMetrics=useMetrics, + epsFactHeaviside=epsFact_consrv_heaviside, + epsFactDirac=epsFact_consrv_dirac, + epsFactDiffusion=epsFact_consrv_diffusion) + +class zero_phi: + def __init__(self): + pass + def uOfX(self,X): + return 0.0 + def uOfXT(self,X,t): + return 0.0 + +initialConditions = {0:zero_phi()} + + + + diff --git a/2d/waves_vegetation/ls_n.py b/2d/waves_vegetation/ls_n.py new file mode 100644 index 00000000..4cfc1db3 --- /dev/null +++ b/2d/waves_vegetation/ls_n.py @@ -0,0 +1,62 @@ +from proteus import * +from ls_p import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*ls_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +conservativeFlux = None +numericalFluxType = NCLS.NumericalFlux +subgridError = NCLS.SubgridError(coefficients,nd) +shockCapturing = NCLS.ShockCapturing(coefficients,nd,shockCapturingFactor=ls_shockCapturingFactor,lag=ls_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'ncls_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = ls_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*ls_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + diff --git a/2d/waves_vegetation/ls_p.py b/2d/waves_vegetation/ls_p.py new file mode 100644 index 00000000..fad3ff3a --- /dev/null +++ b/2d/waves_vegetation/ls_p.py @@ -0,0 +1,29 @@ +from proteus import * +from proteus.default_p import * +from tank import * +from proteus.mprans import NCLS + +LevelModelType = NCLS.LevelModel + +coefficients = NCLS.Coefficients(V_model=0,RD_model=3,ME_model=2, + checkMass=False, useMetrics=useMetrics, + epsFact=epsFact_consrv_heaviside,sc_uref=ls_sc_uref,sc_beta=ls_sc_beta,movingDomain=movingDomain) + +def getDBC_ls(x,flag): + if flag == boundaryTags['left']: + return wavePhi +# elif flag == boundaryTags['right']: +# return outflowPhi + else: + return None + +dirichletConditions = {0:getDBC_ls} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/2d/waves_vegetation/redist_n.py b/2d/waves_vegetation/redist_n.py new file mode 100644 index 00000000..237f5e17 --- /dev/null +++ b/2d/waves_vegetation/redist_n.py @@ -0,0 +1,67 @@ +from proteus import * +from redist_p import * +from tank import * + +nl_atol_res = rd_nl_atol_res +tolFac = 0.0 +nl_atol_res = rd_nl_atol_res + +linTolFac = 0.01 +l_atol_res = 0.01*rd_nl_atol_res + +if redist_Newton: + timeIntegration = NoIntegration + stepController = Newton_controller + maxNonlinearIts = 50 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'r' + levelNonlinearSolverConvergenceTest = 'r' + linearSolverConvergenceTest = 'r-true' + useEisenstatWalker = False +else: + timeIntegration = BackwardEuler_cfl + stepController = RDLS.PsiTC + runCFL=2.0 + psitc['nStepsForce']=3 + psitc['nStepsMax']=50 + psitc['reduceRatio']=2.0 + psitc['startRatio']=1.0 + rtol_res[0] = 0.0 + atol_res[0] = rd_nl_atol_res + useEisenstatWalker = False + maxNonlinearIts = 1 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'rits' + levelNonlinearSolverConvergenceTest = 'rits' + linearSolverConvergenceTest = 'r-true' + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +subgridError = RDLS.SubgridError(coefficients,nd) +shockCapturing = RDLS.ShockCapturing(coefficients,nd,shockCapturingFactor=rd_shockCapturingFactor,lag=rd_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = NLGaussSeidel +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rdls_' + diff --git a/2d/waves_vegetation/redist_p.py b/2d/waves_vegetation/redist_p.py new file mode 100644 index 00000000..5700202e --- /dev/null +++ b/2d/waves_vegetation/redist_p.py @@ -0,0 +1,32 @@ +from proteus import * +from proteus.default_p import * +from math import * +from tank import * +from proteus.mprans import RDLS +""" +The redistancing equation in the sloshbox test problem. +""" + +LevelModelType = RDLS.LevelModel + +coefficients = RDLS.Coefficients(applyRedistancing=applyRedistancing, + epsFact=epsFact_redistance, + nModelId=2, + rdModelId=3, + useMetrics=useMetrics, + backgroundDiffusionFactor=backgroundDiffusionFactor) + +def getDBC_rd(x,flag): + pass + +dirichletConditions = {0:getDBC_rd} +weakDirichletConditions = {0:RDLS.setZeroLSweakDirichletBCsSimple} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/2d/waves_vegetation/tank.py b/2d/waves_vegetation/tank.py new file mode 100644 index 00000000..5558e2b6 --- /dev/null +++ b/2d/waves_vegetation/tank.py @@ -0,0 +1,579 @@ +from math import * +import proteus.MeshTools +from proteus import Domain +from proteus.default_n import * +from proteus.Profiling import logEvent +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.ctransportCoefficients import smoothedHeaviside_integral +from proteus import Gauges +from proteus.Gauges import PointGauges,LineGauges,LineIntegralGauges +from proteus import WaveTools as WT + +from proteus import Context +opts=Context.Options([ + ("wave_type", 'linear', "type of waves generated: 'linear', 'Nonlinear', 'single-peaked', 'double-peaked'"), + ("depth", 0.457, "water depth at leading edge of vegetation (not at wave generator)[m]"), + ("wave_height", 0.192, "wave height at leading edget of vegetation [m]"), + ("peak_period", 2.0, "Peak period [s]"), + ("peak_period2", 1.5, "Second peak period (only used in double-peaked case)[s]"), + ("peak_wavelength",3.91,"Peak wavelength in [m]"), + ("parallel", False, "Run in parallel"), + ("gauges", True, "Enable gauges")]) + +#wave generator +windVelocity = (0.0,0.0) +veg_platform_height = 17.2/44.0 + 6.1/20.0 +depth = veg_platform_height + opts.depth +inflowHeightMean = depth +inflowVelocityMean = (0.0,0.0) +period = opts.peak_period +omega = 2.0*math.pi/period +waveheight = opts.wave_height +amplitude = waveheight/ 2.0 +wavelength = opts.peak_wavelength +k = 2.0*math.pi/wavelength +waveDir = numpy.array([1,0,0]) +g = numpy.array([0,-9.81,0]) +if opts.wave_type == 'linear': + waves = WT.MonochromaticWaves(period = period, # Peak period + waveHeight = waveheight, # Height + depth = depth, # Depth + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + waveType="Linear") +elif opts.wave_type == 'Nonlinear': + waves = WT.MonochromaticWaves(period = period, # Peak period + waveHeight = waveheight, # Height + wavelength = wavelength, + depth = depth, # Depth + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + waveType="Fenton", + Ycoeff = [0.04160592, #Surface elevation Fourier coefficients for non-dimensionalised solution + 0.00555874, + 0.00065892, + 0.00008144, + 0.00001078, + 0.00000151, + 0.00000023, + 0.00000007], + Bcoeff = [0.05395079, + 0.00357780, + 0.00020506, + 0.00000719, + -0.00000016, + -0.00000005, + 0.00000000, + 0.00000000]) +elif opts.wave_type == 'single-peaked': + waves = WT.RandomWaves( Tp = period, # Peak period + Hs = waveheight, # Height + d = depth, # Depth + fp = 1./period, #peak Frequency + bandFactor = 2.0, #fmin=fp/Bandfactor, fmax = Bandfactor * fp + N = 101, #No of frequencies for signal reconstruction + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + gamma=3.3, + spec_fun = WT.JONSWAP) +elif opts.wave_type == 'double-peaked': + waves = WT.DoublePeakedRandomWaves( Tp = period, # Peak period + Hs = waveheight, # Height + d = depth, # Depth + fp = 1./period, #peak Frequency + bandFactor = 2.0, #fmin=fp/Bandfactor, fmax = Bandfactor * fp + N = 101, #No of frequencies for signal reconstruction + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + gamma=10.0, + spec_fun = WT.JONSWAP, + Tp_2 = opts.peak_period2) + +gauges=opts.gauges +# Discretization -- input options +genMesh=True +movingDomain=False +applyRedistancing=True +useOldPETSc=False +useSuperlu=not opts.parallel +timeDiscretization='be'#'be','vbdf','flcbdf' +spaceOrder = 1 +useHex = False +useRBLES = 0.0 +useMetrics = 1.0 +applyCorrection=True +useVF = 1.0 +useOnlyVF = False +useRANS = 0 # 0 -- None + # 1 -- K-Epsilon + # 2 -- K-Omega +# Input checks +if spaceOrder not in [1,2]: + print "INVALID: spaceOrder" + spaceOrder + sys.exit() + +if useRBLES not in [0.0, 1.0]: + print "INVALID: useRBLES" + useRBLES + sys.exit() + +if useMetrics not in [0.0, 1.0]: + print "INVALID: useMetrics" + sys.exit() + +# Discretization +nd = 2 +if spaceOrder == 1: + hFactor=1.0 + if useHex: + basis=C0_AffineLinearOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,2) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,2) + else: + basis=C0_AffineLinearOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,3) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,3) +elif spaceOrder == 2: + hFactor=0.5 + if useHex: + basis=C0_AffineLagrangeOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,4) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,4) + else: + basis=C0_AffineQuadraticOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,4) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,4) + +# Domain and mesh + +#for debugging, make the tank short +L = (45.4,1.0) +he = float(wavelength)/50.0#0.0#100 + +GenerationZoneLength = wavelength +AbsorptionZoneLength= 45.4-37.9 +spongeLayer = True +xSponge = GenerationZoneLength +xRelaxCenter = xSponge/2.0 +epsFact_solid = xSponge/2.0 +#zone 2 +xSponge_2 = 37.9 +xRelaxCenter_2 = 0.5*(37.9+45.4) +epsFact_solid_2 = AbsorptionZoneLength/2.0 + +weak_bc_penalty_constant = 100.0 +nLevels = 1 +#parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.element +parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.node +nLayersOfOverlapForParallel = 0 +structured=False + + +gauge_dx=0.075 +PGL=[] +gauge_top = 19.5/44.0 + 1.5 +veg_gauge_bottom_y = 17.2/44.0 + 6.1/20.0 +LGL=[[(6.1, (6.1 - 5.4)/44.0, 0.0), (6.1,gauge_top,0.0)],#1 Goda + [(6.4, (6.4 - 5.4)/44.0, 0.0), (6.4,gauge_top,0.0)],#2 Goda + [(7.0, (7.0 - 5.4)/44.0, 0.0), (7.0,gauge_top,0.0)],#3 Goda + [(26.0, veg_gauge_bottom_y, 0.0), (26.0,gauge_top,0.0)],#4 veg + [(26.9, veg_gauge_bottom_y, 0.0), (26.9,gauge_top,0.0)],#5 + [(27.4, veg_gauge_bottom_y, 0.0), (27.4,gauge_top,0.0)],#6 + [(27.9, veg_gauge_bottom_y, 0.0), (27.9,gauge_top,0.0)],#7 + [(28.5, veg_gauge_bottom_y, 0.0), (28.5,gauge_top,0.0)],#8 + [(29.5, veg_gauge_bottom_y, 0.0), (29.5,gauge_top,0.0)],#9 + [(31.0, veg_gauge_bottom_y, 0.0), (31.0,gauge_top,0.0)],#10 + [(32.7, veg_gauge_bottom_y, 0.0), (32.7,gauge_top,0.0)],#11 + [(34.4, veg_gauge_bottom_y, 0.0), (34.4,gauge_top,0.0)],#12 + [(36.2, veg_gauge_bottom_y, 0.0), (36.2,gauge_top,0.0)]]#13 +#for i in range(0,int(L[0]/gauge_dx+1)): #+1 only if gauge_dx is an exact +# PGL.append([gauge_dx*i,0.5,0]) +# LGL.append([(gauge_dx*i,0.0,0),(gauge_dx*i,L[1],0)]) + + +gaugeLocations=tuple(map(tuple,PGL)) +columnLines=tuple(map(tuple,LGL)) + + +#pointGauges = PointGauges(gauges=((('u','v'), gaugeLocations), +# (('p',), gaugeLocations)), +# activeTime = (0, 1000.0), +# sampleRate = 0, +# fileName = 'combined_gauge_0_0.5_sample_all.txt') + +#print gaugeLocations +#print columnLines + +fields = ('vof',) + +columnGauge = LineIntegralGauges(gauges=((fields, columnLines),), + fileName='column_gauge.csv') + +#v_resolution = max(he,0.05) +#linePoints = int((gauge_top - veg_gauge_bottom_y)/v_resolution) +lineGauges = LineGauges(gauges=((('u','v'),#fields in gauge set + (#lines for these fields + ((26.0, veg_gauge_bottom_y, 0.0),(26.0, gauge_top, 0.0)), + ),#end lines + ),#end gauge set + ),#end gauges + fileName="vegZoneVelocity.csv") + +#lineGauges_phi = LineGauges_phi(lineGauges.endpoints,linePoints=20) + + +if useHex: + nnx=ceil(L[0]/he)+1 + nny=ceil(L[1]/he)+1 + hex=True + domain = Domain.RectangularDomain(L) +else: + boundaries=['left','right','bottom','top','front','back'] + boundaryTags=dict([(key,i+1) for (i,key) in enumerate(boundaries)]) + if structured: + nnx=ceil(L[0]/he)+1 + nny=ceil(L[1]/he)+1 + elif spongeLayer: + vertices=[[0.0, 0.0 ],#0 + [5.4, 0.0 ],#1 + [5.4 + 17.2, 17.2/44.0 ],#2 + [5.4 + 17.2 + 6.1, 17.2/44.0 + 6.1/20.0 ],#3 + [5.4 + 17.2 + 6.1 + 1.2, 17.2/44.0 + 6.1/20.0 ],#4 + [5.4 + 17.2 + 6.1 + 1.2 + 9.8, 17.2/44.0 + 6.1/20.0 ],#5 + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2, 17.2/44.0 + 6.1/20.0 ],#6 -- sponge + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2 + 20.95, 17.2/44.0 + 6.1/20.0 + 20.95/20.0],#7 + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2 + 20.95, 17.2/44.0 + 6.1/20.0 + 20.95/20.0+0.2],#8 + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2, 19.5/44.0 + 1.5],#9 -- sponge + [0.0, 19.5/44.0 + 1.5]]#10 + + vertexFlags=[boundaryTags['bottom'],#0 + boundaryTags['bottom'],#1 + boundaryTags['bottom'],#2 + boundaryTags['bottom'],#3 + boundaryTags['bottom'],#4 + boundaryTags['bottom'],#5 + boundaryTags['bottom'],#6 + boundaryTags['bottom'],#7 + boundaryTags['top'],#8 + boundaryTags['top'],#9 + boundaryTags['top']]#10 + segments=[[0,1],#0 + [1,2],#1 + [2,3],#2 + [3,4],#3 + [4,5],#4 + [5,6],#5 + [6,7],#6 + [7,8],#7 + [8,9],#8 + [9,10],#8 + [10,0],#9 + [6,9]]#10 + + segmentFlags=[boundaryTags['bottom'],#0 + boundaryTags['bottom'],#1 + boundaryTags['bottom'],#2 + boundaryTags['bottom'],#3 + boundaryTags['bottom'],#4 + boundaryTags['bottom'],#5 + boundaryTags['bottom'],#6 + boundaryTags['right'],#7 + boundaryTags['top'],#8 + boundaryTags['top'],#9 + boundaryTags['left'],#10 + 0]#11 + + regions=[[0.5,0.5], + [5.4 + 17.2 + 6.1 + 1.2 + 9.8 + 1.2+1.0, 17.2/44.0 + 6.1/20.0 + 1.0]] + regionFlags=[1,2] + domain = Domain.PlanarStraightLineGraphDomain(vertices=vertices, + vertexFlags=vertexFlags, + segments=segments, + segmentFlags=segmentFlags, + regions=regions, + regionFlags=regionFlags) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="VApq30Dena%8.8f" % ((he**2)/2.0,) + + logEvent("""Mesh generated using: triangle -%s %s""" % (triangleOptions,domain.polyfile+".poly")) + porosityTypes = numpy.array([1.0, + 1.0, + 1.0]) + dragAlphaTypes = numpy.array([0.0, + 0.0, + 0.5/1.004e-6]) + dragBetaTypes = numpy.array([0.0,0.0,0.0]) + + epsFact_solidTypes = np.array([0.0,0.0,epsFact_solid_2]) + + else: + vertices=[[0.0,0.0],#0 + [L[0],0.0],#1 + [L[0],L[1]],#2 + [0.0,L[1]]]#3 + + vertexFlags=[boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['top'], + boundaryTags['top']] + segments=[[0,1], + [1,2], + [2,3], + [3,0] + ] + segmentFlags=[boundaryTags['bottom'], + boundaryTags['right'], + boundaryTags['top'], + boundaryTags['left']] + + regions=[ [ 0.1*L[0] , 0.1*L[1] ], + [0.95*L[0] , 0.95*L[1] ] ] + regionFlags=[1,2] + domain = Domain.PlanarStraightLineGraphDomain(vertices=vertices, + vertexFlags=vertexFlags, + segments=segments, + segmentFlags=segmentFlags, + regions=regions, + regionFlags=regionFlags) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="VApq30Dena%8.8f" % ((he**2)/2.0,) + + logEvent("""Mesh generated using: triangle -%s %s""" % (triangleOptions,domain.polyfile+".poly")) +# Time stepping +T=40*period +dt_fixed = period/11.0#2.0*0.5/20.0#T/2.0#period/21.0 +dt_init = min(0.001*dt_fixed,0.001) +runCFL=0.90 +nDTout = int(round(T/dt_fixed)) + +# Numerical parameters +ns_forceStrongDirichlet = False +backgroundDiffusionFactor=0.0 +if useMetrics: + ns_shockCapturingFactor = 0.25 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.25 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.25 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.25 + rd_lag_shockCapturing = False + epsFact_density = 3.0 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 1.5 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.1 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.1 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 +else: + ns_shockCapturingFactor = 0.9 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.9 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.9 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.9 + rd_lag_shockCapturing = False + epsFact_density = 1.5 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 0.33 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.9 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.9 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 + +ns_nl_atol_res = max(1.0e-10,0.001*he**2) +vof_nl_atol_res = max(1.0e-10,0.001*he**2) +ls_nl_atol_res = max(1.0e-10,0.001*he**2) +rd_nl_atol_res = max(1.0e-10,0.05*he) +mcorr_nl_atol_res = max(1.0e-10,0.001*he**2) +kappa_nl_atol_res = max(1.0e-10,0.001*he**2) +dissipation_nl_atol_res = max(1.0e-10,0.001*he**2) + +#turbulence +ns_closure=0 #1-classic smagorinsky, 2-dynamic smagorinsky, 3 -- k-epsilon, 4 -- k-omega +if useRANS == 1: + ns_closure = 3 +elif useRANS == 2: + ns_closure == 4 +# Water +rho_0 = 998.2 +nu_0 = 1.004e-6 + +# Air +rho_1 = 1.205 +nu_1 = 1.500e-5 + +# Surface tension +sigma_01 = 0.0 + +# Gravity +g = [0.0,-9.8] + +# Initial condition +waterLine_x = 2*L[0] +waterLine_z = inflowHeightMean + + +def signedDistance(x): + phi_x = x[0]-waterLine_x + phi_z = x[1]-waterLine_z + if phi_x < 0.0: + if phi_z < 0.0: + return max(phi_x,phi_z) + else: + return phi_z + else: + if phi_z < 0.0: + return phi_x + else: + return sqrt(phi_x**2 + phi_z**2) + + +def theta(x,t): + return k*x[0] - omega*t + pi/2.0 + +def z(x): + return x[1] - inflowHeightMean + +#sigma = omega - k*inflowVelocityMean[0] +h = inflowHeightMean # - transect[0][1] if lower left hand corner is not at z=0 + +#waveData + +def waveHeight(x,t): + return inflowHeightMean + waves.eta(x[0],x[1],x[2],t) +def waveVelocity_u(x,t): + return waves.u(x[0],x[1],x[2],t,"x") +def waveVelocity_v(x,t): + return waves.u(x[0],x[1],x[2],t,"y") +def waveVelocity_w(x,t): + return waves.u(x[0],x[1],x[2],t,"z") + +#solution variables + +def wavePhi(x,t): + return x[1] - waveHeight(x,t) + +def waveVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)) + +def twpflowVelocity_u(x,t): + waterspeed = waveVelocity_u(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + u = H*windVelocity[0] + (1.0-H)*waterspeed + return u + +def twpflowVelocity_v(x,t): + waterspeed = waveVelocity_v(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + return H*windVelocity[1]+(1.0-H)*waterspeed + +def twpflowFlux(x,t): + return -twpflowVelocity_u(x,t) + +outflowHeight=inflowHeightMean + +def outflowVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,x[1] - outflowHeight) + +def outflowPhi(x,t): + return x[1] - outflowHeight + +def outflowPressure(x,t): + if x[1]>inflowHeightMean: + return (L[1]-x[1])*rho_1*abs(g[1]) + else: + return (L[1]-inflowHeightMean)*rho_1*abs(g[1])+(inflowHeightMean-x[1])*rho_0*abs(g[1]) + + + #p_L = L[1]*rho_1*g[1] + #phi_L = L[1] - outflowHeight + #phi = x[1] - outflowHeight + #return p_L -g[1]*(rho_0*(phi_L - phi)+(rho_1 -rho_0)*(smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi_L) + # -smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi))) + +def twpflowVelocity_w(x,t): + return 0.0 + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + +RelaxationZone = namedtuple("RelaxationZone","center_x sign u v w") + +class RelaxationZoneWaveGenerator(AV_base): + """ Prescribe a velocity penalty scaling in a material zone via a Darcy-Forchheimer penalty + + :param zones: A dictionary mapping integer material types to Zones, where a Zone is a named tuple + specifying the x coordinate of the zone center and the velocity components + """ + def __init__(self,zones): + assert isinstance(zones,dict) + self.zones = zones + def calculate(self): + for l,m in enumerate(self.model.levelModelList): + for eN in range(m.coefficients.q_phi.shape[0]): + mType = m.mesh.elementMaterialTypes[eN] + if self.zones.has_key(mType): + for k in range(m.coefficients.q_phi.shape[1]): + t = m.timeIntegration.t + x = m.q['x'][eN,k] + m.coefficients.q_phi_solid[eN,k] = self.zones[mType].sign*(self.zones[mType].center_x - x[0]) + m.coefficients.q_velocity_solid[eN,k,0] = self.zones[mType].u(x,t) + m.coefficients.q_velocity_solid[eN,k,1] = self.zones[mType].v(x,t) + #m.coefficients.q_velocity_solid[eN,k,2] = self.zones[mType].w(x,t) + m.q['phi_solid'] = m.coefficients.q_phi_solid + m.q['velocity_solid'] = m.coefficients.q_velocity_solid + +rzWaveGenerator = RelaxationZoneWaveGenerator(zones={ + # 1:RelaxationZone(xRelaxCenter, + # 1.0, + # twpflowVelocity_u, + # twpflowVelocity_v, + # twpflowVelocity_w), + 2:RelaxationZone(xRelaxCenter_2, + -1.0, #currently Hs=1-exp_function + zeroVel, + zeroVel, + zeroVel)}) diff --git a/2d/waves_vegetation/tank_batch.py b/2d/waves_vegetation/tank_batch.py new file mode 100644 index 00000000..38fabd1d --- /dev/null +++ b/2d/waves_vegetation/tank_batch.py @@ -0,0 +1,4 @@ +simFlagsList[0]['storeQuantities']= ["q:'phi_solid'","q:'velocity_solid'"] +#simFlagsList[0]['storeQuantities']= ["q:velocity_solid"] +start +quit diff --git a/2d/waves_vegetation/tank_so.py b/2d/waves_vegetation/tank_so.py new file mode 100644 index 00000000..d06f6ffa --- /dev/null +++ b/2d/waves_vegetation/tank_so.py @@ -0,0 +1,35 @@ +from proteus.default_so import * +import tank +from proteus import Context +Context.setFromModule(tank) +ctx = Context.get() + +if tank.useOnlyVF: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n")] +else: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n"), + ("ls_p", "ls_n"), + ("redist_p", "redist_n"), + ("ls_consrv_p", "ls_consrv_n")] + + +if tank.useRANS > 0: + pnList.append(("kappa_p", + "kappa_n")) + pnList.append(("dissipation_p", + "dissipation_n")) +name = "tank_p" + +if tank.timeDiscretization == 'flcbdf': + systemStepControllerType = Sequential_MinFLCBDFModelStep + systemStepControllerType = Sequential_MinAdaptiveModelStep +else: + systemStepControllerType = Sequential_MinAdaptiveModelStep + +needEBQ_GLOBAL = False +needEBQ = False + +tnList = [0.0,tank.dt_init]+[i*tank.dt_fixed for i in range(1,tank.nDTout+1)] +archiveFlag = ArchiveFlags.EVERY_USER_STEP diff --git a/2d/waves_vegetation/twp_navier_stokes_n.py b/2d/waves_vegetation/twp_navier_stokes_n.py new file mode 100644 index 00000000..fe79dd36 --- /dev/null +++ b/2d/waves_vegetation/twp_navier_stokes_n.py @@ -0,0 +1,67 @@ +from proteus import * +from twp_navier_stokes_p import * +from tank import * +from proteus import Context +ctx = Context.get() + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller_sys + stepController = Min_dt_cfl_controller + time_tol = 10.0*ns_nl_atol_res + atol_u = {1:time_tol,2:time_tol} + rtol_u = {1:time_tol,2:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis, + 1:basis, + 2:basis} + +massLumping = False +numericalFluxType = None +conservativeFlux = None + +numericalFluxType = RANS2P.NumericalFlux +subgridError = RANS2P.SubgridError(coefficients,nd,lag=ns_lag_subgridError,hFactor=hFactor,nStepsToDelay=1) +shockCapturing = RANS2P.ShockCapturing(coefficients,nd,ns_shockCapturingFactor,lag=ns_lag_shockCapturing,nStepsToDelay=1) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None + +#linearSmoother = SimpleNavierStokes2D + +matrix = SparseMatrix + +multilevelLinearSolver = KSP_petsc4py +levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rans2p_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*ns_nl_atol_res +nl_atol_res = ns_nl_atol_res +useEisenstatWalker = False +maxNonlinearIts = 100 +maxLineSearches = 0 +conservativeFlux = {0:'pwl-bdm-opt'} + +if ctx.gauges: + auxiliaryVariables=[lineGauges] +#auxiliaryVariables=[pointGauges,rzWaveGenerator] diff --git a/2d/waves_vegetation/twp_navier_stokes_p.py b/2d/waves_vegetation/twp_navier_stokes_p.py new file mode 100644 index 00000000..10a8a9dc --- /dev/null +++ b/2d/waves_vegetation/twp_navier_stokes_p.py @@ -0,0 +1,141 @@ +from proteus import * +from proteus.default_p import * +from tank import * +from proteus.mprans import RANS2P + +LevelModelType = RANS2P.LevelModel +if useOnlyVF: + LS_model = None +else: + LS_model = 2 +if useRANS >= 1: + Closure_0_model = 5; Closure_1_model=6 + if useOnlyVF: + Closure_0_model=2; Closure_1_model=3 + if movingDomain: + Closure_0_model += 1; Closure_1_model += 1 +else: + Closure_0_model = None + Closure_1_model = None + +if spongeLayer or levee or slopingSpongeLayer: + coefficients = RANS2P.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain, + porosityTypes=porosityTypes, + dragAlphaTypes=dragAlphaTypes, + dragBetaTypes=dragBetaTypes, + epsFact_solid = epsFact_solidTypes) +else: + coefficients = RANS2P.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain) + +def getDBC_p(x,flag): + if flag == boundaryTags['top']: + return lambda x,t: 0.0 +# elif flag == boundaryTags['right']: +# return outflowPressure + +def getDBC_u(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_u +# elif flag == boundaryTags['right']: +# return lambda x,t: 0.0 + +def getDBC_v(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_v +# elif flag == boundaryTags['right']: +# return lambda x,t: 0.0 + +dirichletConditions = {0:getDBC_p, + 1:getDBC_u, + 2:getDBC_v} + +def getAFBC_p(x,flag): + if flag == boundaryTags['left']: + return twpflowFlux + elif flag == boundaryTags['bottom'] or flag == boundaryTags['right'] or flag == 0: + return lambda x,t: 0.0 + +def getAFBC_u(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['right'] or flag == 0: + return lambda x,t: 0.0 + +def getAFBC_v(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['right'] or flag == 0: + return lambda x,t: 0.0 + +def getDFBC_u(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + +def getDFBC_v(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + +advectiveFluxBoundaryConditions = {0:getAFBC_p, + 1:getAFBC_u, + 2:getAFBC_v} + +diffusiveFluxBoundaryConditions = {0:{}, + 1:{1:getDFBC_u}, + 2:{2:getDFBC_v}} + +class PerturbedSurface_p: + def __init__(self,waterLevel): + self.waterLevel=waterLevel + def uOfXT(self,x,t): + if signedDistance(x) < 0: + return -(L[1] - self.waterLevel)*rho_1*g[1] - (self.waterLevel - x[1])*rho_0*g[1] + else: + return -(L[1] - self.waterLevel)*rho_1*g[1] + +class AtRest: + def __init__(self): + pass + def uOfXT(self,x,t): + return 0.0 + +initialConditions = {0:PerturbedSurface_p(waterLine_z), + 1:AtRest(), + 2:AtRest()} diff --git a/2d/waves_vegetation/vof_n.py b/2d/waves_vegetation/vof_n.py new file mode 100644 index 00000000..5a15d474 --- /dev/null +++ b/2d/waves_vegetation/vof_n.py @@ -0,0 +1,68 @@ +from proteus import * +from tank import * +from vof_p import * +from proteus import Context +ctx = Context.get() + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*vof_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = VOF.NumericalFlux +conservativeFlux = None +subgridError = VOF.SubgridError(coefficients=coefficients,nd=nd) +shockCapturing = VOF.ShockCapturing(coefficients,nd,shockCapturingFactor=vof_shockCapturingFactor,lag=vof_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'vof_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = vof_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*vof_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + + +if ctx.gauges: + auxiliaryVariables = [columnGauge] diff --git a/2d/waves_vegetation/vof_p.py b/2d/waves_vegetation/vof_p.py new file mode 100644 index 00000000..94965ec8 --- /dev/null +++ b/2d/waves_vegetation/vof_p.py @@ -0,0 +1,46 @@ +from proteus import * +from proteus.default_p import * +from proteus.ctransportCoefficients import smoothedHeaviside +from tank import * +from proteus.mprans import VOF + +LevelModelType = VOF.LevelModel +if useOnlyVF: + RD_model = None + LS_model = None +else: + RD_model = 3 + LS_model = 2 + +coefficients = VOF.Coefficients(LS_model=LS_model,V_model=0,RD_model=RD_model,ME_model=1, + checkMass=False,useMetrics=useMetrics, + epsFact=epsFact_vof,sc_uref=vof_sc_uref,sc_beta=vof_sc_beta,movingDomain=movingDomain) + +def getDBC_vof(x,flag): + if flag == boundaryTags['left']: + return waveVF + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return lambda x,t: 1.0 +# elif flag == boundaryTags['right']: +# return outflowVF + +dirichletConditions = {0:getDBC_vof} + +def getAFBC_vof(x,flag): + if flag == boundaryTags['left']: + return None + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return None + # elif flag == boundaryTags['right']: + # return None + else: + return lambda x,t: 0.0 + +advectiveFluxBoundaryConditions = {0:getAFBC_vof} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_H: + def uOfXT(self,x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,signedDistance(x)) + +initialConditions = {0:PerturbedSurface_H()} diff --git a/3d/waves_vegetation/generate_wave_series/generate_monochromatic_linear.py b/3d/waves_vegetation/generate_wave_series/generate_monochromatic_linear.py new file mode 100644 index 00000000..8d3c4ce3 --- /dev/null +++ b/3d/waves_vegetation/generate_wave_series/generate_monochromatic_linear.py @@ -0,0 +1,30 @@ +from proteus import WaveTools +import numpy as np +import math + +mwl = 0.5 +mw = WaveTools.MonochromaticWaves(period = 1.0, + waveHeight = 0.1, + mwl = mwl, + depth = 1.0, + g = np.array([0.0, 0.0, -9.81]), + waveDir = np.array([1.0, 0.0, 0.0]), + wavelength=None, + waveType="Linear", + Ycoeff = None, + Bcoeff =None, + meanVelocity = 0.0, + phi0 = math.pi/2.0) + +f = open('monochromatic_linear_time_series.txt', 'w') +num_points = 10000 +tList = np.linspace(0.0,100.0,num_points) +etaList = np.zeros(tList.shape) +for i in range(num_points): + etaList[i] = mw.eta(0.0,0.0,0.0,tList[i]) + mwl + +combined = np.vstack((tList,etaList)).transpose() +np.savetxt("monochromatic_linear_time_series.txt", combined) + + + diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_consrv_n.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_consrv_n.py new file mode 100644 index 00000000..d6fdb0a9 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_consrv_n.py @@ -0,0 +1,48 @@ +from proteus import * +from tank3D import * +from ls_consrv_p import * + +timeIntegrator = ForwardIntegrator +timeIntegration = NoIntegration + +femSpaces = {0:basis} + +subgridError = None +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +shockCapturing = None + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'mcorr_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*mcorr_nl_atol_res +nl_atol_res = mcorr_nl_atol_res +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_consrv_p.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_consrv_p.py new file mode 100644 index 00000000..ceb09427 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_consrv_p.py @@ -0,0 +1,26 @@ +from proteus import * +from proteus.default_p import * +from tank3D import * +from proteus.mprans import MCorr + +LevelModelType = MCorr.LevelModel + +coefficients = MCorr.Coefficients(LSModel_index=2,V_model=0,me_model=4,VOFModel_index=1, + applyCorrection=applyCorrection,nd=nd,checkMass=False,useMetrics=useMetrics, + epsFactHeaviside=epsFact_consrv_heaviside, + epsFactDirac=epsFact_consrv_dirac, + epsFactDiffusion=epsFact_consrv_diffusion) + +class zero_phi: + def __init__(self): + pass + def uOfX(self,X): + return 0.0 + def uOfXT(self,X,t): + return 0.0 + +initialConditions = {0:zero_phi()} + + + + diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_n.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_n.py new file mode 100644 index 00000000..fde02e06 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_n.py @@ -0,0 +1,63 @@ +from proteus import * +from ls_p import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*ls_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +conservativeFlux = None +numericalFluxType = NCLS.NumericalFlux +subgridError = NCLS.SubgridError(coefficients,nd) +shockCapturing = NCLS.ShockCapturing(coefficients,nd,shockCapturingFactor=ls_shockCapturingFactor,lag=ls_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'ncls_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = ls_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*ls_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + + diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_p.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_p.py new file mode 100644 index 00000000..9d6237c1 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/ls_p.py @@ -0,0 +1,29 @@ +from proteus import * +from proteus.default_p import * +from tank3D import * +from proteus.mprans import NCLS + +LevelModelType = NCLS.LevelModel + +coefficients = NCLS.Coefficients(V_model=0,RD_model=3,ME_model=2, + checkMass=False, useMetrics=useMetrics, + epsFact=epsFact_consrv_heaviside,sc_uref=ls_sc_uref,sc_beta=ls_sc_beta,movingDomain=movingDomain) + +def getDBC_ls(x,flag): + if flag == boundaryTags['left']: + return wavePhi +# elif flag == boundaryTags['right']: +# return outflowPhi + else: + return None + +dirichletConditions = {0:getDBC_ls} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/redist_n.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/redist_n.py new file mode 100644 index 00000000..4f85688b --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/redist_n.py @@ -0,0 +1,66 @@ +from proteus import * +from redist_p import * +from tank3D import * + +nl_atol_res = rd_nl_atol_res +tolFac = 0.0 +nl_atol_res = rd_nl_atol_res + +linTolFac = 0.01 +l_atol_res = 0.01*rd_nl_atol_res + +if redist_Newton: + timeIntegration = NoIntegration + stepController = Newton_controller + maxNonlinearIts = 50 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'r' + levelNonlinearSolverConvergenceTest = 'r' + linearSolverConvergenceTest = 'r-true' + useEisenstatWalker = False +else: + timeIntegration = BackwardEuler_cfl + stepController = RDLS.PsiTC + runCFL=2.0 + psitc['nStepsForce']=3 + psitc['nStepsMax']=50 + psitc['reduceRatio']=2.0 + psitc['startRatio']=1.0 + rtol_res[0] = 0.0 + atol_res[0] = rd_nl_atol_res + useEisenstatWalker = False + maxNonlinearIts = 1 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'rits' + levelNonlinearSolverConvergenceTest = 'rits' + linearSolverConvergenceTest = 'r-true' + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +subgridError = RDLS.SubgridError(coefficients,nd) +shockCapturing = RDLS.ShockCapturing(coefficients,nd,shockCapturingFactor=rd_shockCapturingFactor,lag=rd_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = NLGaussSeidel +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rdls_' diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/redist_p.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/redist_p.py new file mode 100644 index 00000000..ea199426 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/redist_p.py @@ -0,0 +1,32 @@ +from proteus import * +from proteus.default_p import * +from math import * +from tank3D import * +from proteus.mprans import RDLS +""" +The redistancing equation in the sloshbox test problem. +""" + +LevelModelType = RDLS.LevelModel + +coefficients = RDLS.Coefficients(applyRedistancing=applyRedistancing, + epsFact=epsFact_redistance, + nModelId=2, + rdModelId=3, + useMetrics=useMetrics, + backgroundDiffusionFactor=backgroundDiffusionFactor) + +def getDBC_rd(x,flag): + pass + +dirichletConditions = {0:getDBC_rd} +weakDirichletConditions = {0:RDLS.setZeroLSweakDirichletBCsSimple} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/tank3D.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/tank3D.py new file mode 100644 index 00000000..bd12a064 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/tank3D.py @@ -0,0 +1,586 @@ +from math import * +import proteus.MeshTools +from proteus import Domain +from proteus.default_n import * +from proteus.Profiling import logEvent +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.ctransportCoefficients import smoothedHeaviside_integral +from proteus import Gauges +from proteus.Gauges import PointGauges,LineGauges,LineIntegralGauges +import vegZoneVelocityInterp as vegZoneInterp + + +from proteus import Context +opts=Context.Options([ + ("wave_type", 'linear', "type of waves generated: 'linear', 'Nonlinear', 'single-peaked', 'double-peaked'"), + ("depth", 0.457, "water depth at leading edge of vegetation (not at wave generator)[m]"), + ("wave_height", 0.192, "wave height at leading edget of vegetation [m]"), + ("peak_period", 2.0, "Peak period [s]"), + ("peak_period2", 1.5, "Second peak period (only used in double-peaked case)[s]"), + ("peak_wavelength",3.91,"Peak wavelength in [m]"), + ("parallel", False, "Run in parallel")]) + +#wave generator +windVelocity = (0.0,0.0,0.0) +#veg_platform_height = 17.2/44.0 + 6.1/20.0 +depth = opts.depth +inflowHeightMean = depth +inflowVelocityMean = (0.0,0.0,0,0) +period = opts.peak_period +omega = 2.0*math.pi/period +waveheight = opts.wave_height +amplitude = waveheight/ 2.0 +wavelength = opts.peak_wavelength +k = 2.0*math.pi/wavelength +waveDir = numpy.array([1,0,0]) +g = numpy.array([0.0,0.0,-9.81]) + +# Discretization -- input options +genMesh=True +movingDomain=False +applyRedistancing=True +useOldPETSc=False +useSuperlu=not opts.parallel +timeDiscretization='be'#'be','vbdf','flcbdf' +spaceOrder = 1 +useHex = False +useRBLES = 0.0 +useMetrics = 1.0 +applyCorrection=True +useVF = 1.0 +useOnlyVF = False +useRANS = 0 # 0 -- None + # 1 -- K-Epsilon + # 2 -- K-Omega +# Input checks +if spaceOrder not in [1,2]: + print "INVALID: spaceOrder" + spaceOrder + sys.exit() + +if useRBLES not in [0.0, 1.0]: + print "INVALID: useRBLES" + useRBLES + sys.exit() + +if useMetrics not in [0.0, 1.0]: + print "INVALID: useMetrics" + sys.exit() + +# Discretization +nd = 3 +if spaceOrder == 1: + hFactor=1.0 + if useHex: + basis=C0_AffineLinearOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,2) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,2) + else: + basis=C0_AffineLinearOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,3) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,3) +elif spaceOrder == 2: + hFactor=0.5 + if useHex: + basis=C0_AffineLagrangeOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,4) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,4) + else: + basis=C0_AffineQuadraticOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,4) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,4) + +# Domain and mesh + +#for debugging, make the tank short + +he = float(wavelength)/30.0#0.0#100 +L = (15.0, 1.5, 1.0) + +GenerationZoneLength = 1.2 +AbsorptionZoneLength= 2.8 +spongeLayer = True +xSponge = GenerationZoneLength +xRelaxCenter = xSponge/2.0 +epsFact_solid = xSponge/2.0 +#zone 2 +xSponge_2 = L[0] - AbsorptionZoneLength +xRelaxCenter_2 = 0.5*(xSponge_2+L[0]) +epsFact_solid_2 = AbsorptionZoneLength/2.0 + +weak_bc_penalty_constant = 100.0 +nLevels = 1 +#parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.element +parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.node +nLayersOfOverlapForParallel = 0 +structured=False + +if useHex: + nnx=4*Refinement+1 + nny=2*Refinement+1 + hex=True + domain = Domain.RectangularDomain(L) +else: + boundaries=['empty','left','right','bottom','top','front','back'] + boundaryTags=dict([(key,i+1) for (i,key) in enumerate(boundaries)]) + if structured: + nnx=4*Refinement + nny=2*Refinement + domain = Domain.RectangularDomain(L) + elif spongeLayer: + vertices=[[0.0,0.0,0.0],#0 + [xSponge,0.0,0.0],#1 + [xSponge_2,0.0,0.0],#2 + [L[0],0.0,0.0],#3 + [L[0],L[1],0.0],#4 + [xSponge_2,L[1],0.0],#5 + [xSponge,L[1],0.0],#6 + [0.0,L[1],0.0]]#7 + + + vertexFlags=[boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom']] + + + for v,vf in zip(vertices,vertexFlags): + vertices.append([v[0],v[1],L[2]]) + vertexFlags.append(boundaryTags['top']) + + print vertices + print vertexFlags + + segments=[[0,1], + [1,2], + [2,3], + [3,4], + [4,5], + [5,6], + [6,7], + [7,0], + [1,6], + [2,5]] + + segmentFlags=[boundaryTags['front'], + boundaryTags['front'], + boundaryTags['front'], + boundaryTags['right'], + boundaryTags['back'], + boundaryTags['back'], + boundaryTags['back'], + boundaryTags['left'], + boundaryTags['empty'], + boundaryTags['empty'] ] + + + facets=[] + facetFlags=[] + + for s,sF in zip(segments,segmentFlags): + facets.append([[s[0],s[1],s[1]+8,s[0]+8]]) + facetFlags.append(sF) + + bf=[[0,1,6,7],[1,2,5,6],[2,3,4,5]] + tf=[] + for i in range(0,3): + facets.append([bf[i]]) + tf=[ss + 8 for ss in bf[i]] + facets.append([tf]) + + for i in range(0,3): + facetFlags.append(boundaryTags['bottom']) + facetFlags.append(boundaryTags['top']) + + print facets + print facetFlags + + regions=[[xRelaxCenter, 0.5*L[1],0.0], + [xRelaxCenter_2, 0.5*L[1], 0.0], + [0.5*L[0],0.1*L[1], 0.0]] + regionFlags=[1,2,3] + + domain = Domain.PiecewiseLinearComplexDomain(vertices=vertices, + vertexFlags=vertexFlags, + facets=facets, + facetFlags=facetFlags, + regions=regions, + regionFlags=regionFlags, + ) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="KVApq1.4q12feena%21.16e" % ((he**3)/6.0,) + + + logEvent("""Mesh generated using: tetgen -%s %s""" % (triangleOptions,domain.polyfile+".poly")) + + porosityTypes = numpy.array([1.0, + 1.0, + 1.0, + 1.0]) + dragAlphaTypes = numpy.array([0.0, + 0.0, + 0.0, + 0.5/1.004e-6]) + + dragBetaTypes = numpy.array([0.0,0.0,0.0,0.0]) + + epsFact_solidTypes = np.array([0.0,0.0,0.0,0.0]) + + else: + vertices=[[0.0,0.0,0.0],#0 + [L[0],0.0,0.0],#1 + [L[0],L[1],0.0],#2 + [0.0,L[1],0.0]]#3 + + + vertexFlags=[boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom']] + + + for v,vf in zip(vertices,vertexFlags): + vertices.append([v[0],v[1],L[2]]) + vertexFlags.append(boundaryTags['top']) + + segments=[[0,1], + [1,2], + [2,3], + [3,0]] + + segmentFlags=[boundaryTags['front'], + boundaryTags['right'], + boundaryTags['back'], + boundaryTags['left']] + + facets=[] + facetFlags=[] + + for s,sF in zip(segments,segmentFlags): + facets.append([[s[0],s[1],s[1]+4,s[0]+4]]) + facetFlags.append(sF) + + bf=[[0,1,2,3]] + tf=[] + for i in range(0,1): + facets.append([bf[i]]) + tf=[ss + 4 for ss in bf[i]] + facets.append([tf]) + + for i in range(0,1): + facetFlags.append(boundaryTags['bottom']) + facetFlags.append(boundaryTags['top']) + + for s,sF in zip(segments,segmentFlags): + segments.append([s[1]+4,s[0]+4]) + segmentFlags.append(sF) + + + regions=[[0.5*L[0],0.5*L[1], 0.0]] + regionFlags=[1] + + domain = Domain.PiecewiseLinearComplexDomain(vertices=vertices, + vertexFlags=vertexFlags, + facets=facets, + facetFlags=facetFlags, + regions=regions, + regionFlags=regionFlags) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="KVApq1.4q12feena%21.16e" % ((he**3)/6.0,) + + + logEvent("""Mesh generated using: tetgen -%s %s""" % (triangleOptions,domain.polyfile+".poly")) + + + +# Time stepping +T=40*period +dt_fixed = period/11.0#2.0*0.5/20.0#T/2.0#period/21.0 +dt_init = min(0.001*dt_fixed,0.001) +runCFL=0.90 +nDTout = int(round(T/dt_fixed)) + +# Numerical parameters +ns_forceStrongDirichlet = False#True +backgroundDiffusionFactor=0.0 +if useMetrics: + ns_shockCapturingFactor = 0.25 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.25 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.25 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.25 + rd_lag_shockCapturing = False + epsFact_density = 3.0 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 1.5 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.1 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.1 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 +else: + ns_shockCapturingFactor = 0.9 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.9 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.9 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.9 + rd_lag_shockCapturing = False + epsFact_density = 1.5 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 0.33 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.9 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.9 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 + +ns_nl_atol_res = max(1.0e-10,0.0001*he**2) +vof_nl_atol_res = max(1.0e-10,0.0001*he**2) +ls_nl_atol_res = max(1.0e-10,0.0001*he**2) +rd_nl_atol_res = max(1.0e-10,0.005*he) +mcorr_nl_atol_res = max(1.0e-10,0.0001*he**2) +kappa_nl_atol_res = max(1.0e-10,0.001*he**2) +dissipation_nl_atol_res = max(1.0e-10,0.001*he**2) + +#turbulence +ns_closure=0 #1-classic smagorinsky, 2-dynamic smagorinsky, 3 -- k-epsilon, 4 -- k-omega +if useRANS == 1: + ns_closure = 3 +elif useRANS == 2: + ns_closure == 4 +# Water +rho_0 = 998.2 +nu_0 = 1.004e-6 + +# Air +rho_1 = 1.205 +nu_1 = 1.500e-5 + +# Surface tension +sigma_01 = 0.0 + +# Gravity +g = [0.0,0.0,-9.8] + +# Initial condition +waterLine_x = 2*L[0] +waterLine_z = inflowHeightMean + + +def signedDistance(x): + phi_x = x[0]-waterLine_x + phi_z = x[2]-waterLine_z + if phi_x < 0.0: + if phi_z < 0.0: + return max(phi_x,phi_z) + else: + return phi_z + else: + if phi_z < 0.0: + return phi_x + else: + return sqrt(phi_x**2 + phi_z**2) + + +def theta(x,t): + return k*x[0] - omega*t + pi/2.0 + +def z(x): + return x[2] - inflowHeightMean + +#sigma = omega - k*inflowVelocityMean[0] +h = inflowHeightMean # - transect[0][1] if lower left hand corner is not at z=0 + +#waveData + +def waveHeight(x,t): + return float(vegZoneInterp.interp_phi.__call__(t%vegZoneInterp.time[-1])) +def waveVelocity_u(x,t): + return vegZoneInterp.interpU.__call__(t%vegZoneInterp.time[-1],x[2])[0][0] +def waveVelocity_v(x,t): + return 0.0 +def waveVelocity_w(x,t): + return vegZoneInterp.interpW.__call__(t%vegZoneInterp.time[-1],x[2])[0][0] + +#solution variables + +def wavePhi(x,t): + return x[2] - waveHeight(x,t) + +def waveVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)) + +def twpflowVelocity_u(x,t): + waterspeed = waveVelocity_u(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + u = H*windVelocity[0] + (1.0-H)*waterspeed + return u + +def twpflowVelocity_v(x,t): + waterspeed = waveVelocity_v(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + return H*windVelocity[1]+(1.0-H)*waterspeed + +def twpflowFlux(x,t): + return -twpflowVelocity_u(x,t) + +outflowHeight=inflowHeightMean + +def outflowVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,x[2] - outflowHeight) + +def outflowPhi(x,t): + return x[2] - outflowHeight + +def outflowPressure(x,t): + if x[2]>inflowHeightMean: + return (L[2]-x[2])*rho_1*abs(g[2]) + else: + return (L[2]-inflowHeightMean)*rho_1*abs(g[2])+(inflowHeightMean-x[2])*rho_0*abs(g[2]) + + + #p_L = L[1]*rho_1*g[1] + #phi_L = L[1] - outflowHeight + #phi = x[1] - outflowHeight + #return p_L -g[1]*(rho_0*(phi_L - phi)+(rho_1 -rho_0)*(smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi_L) + # -smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi))) + +def twpflowVelocity_w(x,t): + return 0.0 + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + +RelaxationZone = namedtuple("RelaxationZone","center_x sign u v w") + +class RelaxationZoneWaveGenerator(AV_base): + """ Prescribe a velocity penalty scaling in a material zone via a Darcy-Forchheimer penalty + + :param zones: A dictionary mapping integer material types to Zones, where a Zone is a named tuple + specifying the x coordinate of the zone center and the velocity components + """ + def __init__(self,zones): + assert isinstance(zones,dict) + self.zones = zones + def calculate(self): + for l,m in enumerate(self.model.levelModelList): + for eN in range(m.coefficients.q_phi.shape[0]): + mType = m.mesh.elementMaterialTypes[eN] + if self.zones.has_key(mType): + for k in range(m.coefficients.q_phi.shape[1]): + t = m.timeIntegration.t + x = m.q['x'][eN,k] + m.coefficients.q_phi_solid[eN,k] = self.zones[mType].sign*(self.zones[mType].center_x - x[0]) + m.coefficients.q_velocity_solid[eN,k,0] = self.zones[mType].u(x,t) + m.coefficients.q_velocity_solid[eN,k,1] = self.zones[mType].v(x,t) + #m.coefficients.q_velocity_solid[eN,k,2] = self.zones[mType].w(x,t) + m.q['phi_solid'] = m.coefficients.q_phi_solid + m.q['velocity_solid'] = m.coefficients.q_velocity_solid + +rzWaveGenerator = RelaxationZoneWaveGenerator(zones={ + # 1:RelaxationZone(xRelaxCenter, + # 1.0, + # twpflowVelocity_u, + # twpflowVelocity_v, + # twpflowVelocity_w), + 1:RelaxationZone(xRelaxCenter_2, + -1.0, #currently Hs=1-exp_function + zeroVel, + zeroVel, + zeroVel)}) + +beam_quadOrder=3 +beam_useSparse=False +beamFilename="wavetankBeams" +#nBeamElements=max(nBeamElements,3) + +#beam info +beamLocation=[] +beamLength=[] +beamRadius=[] +EI=[] +GJ=[] +lam = 0.25 #0.05 #3.0*2.54/100.0 #57.4e-3 +lamx = 3.0**0.5*lam +xs = 1.2 +ys = 0.0 +xList=[] +yList = [] +while xs <= 11.0: + xList.append(xs) + xs += lam +while ys<= L[1]: + yList.append(ys) + ys+=lamx +for i in xList: + for j in yList: + beamLocation.append((i,j)) + beamLength.append(0.415) + beamRadius.append(0.0032) + EI.append(3.0e-4) # needs to be fixed + GJ.append(1.5e-4) # needs to be fixed + +xs = 1.2+0.5*lam +ys = 0.5*lamx +xList=[] +yList = [] +while xs <= 11.0: + xList.append(xs) + xs += lam + +while ys<= L[1]: + yList.append(ys) + ys+=lamx + +for i in xList: + for j in yList: + beamLocation.append((i,j)) + beamLength.append(0.415) + beamRadius.append(0.0032) + EI.append(3.0e-4) # needs to be fixed + GJ.append(1.5e-4) # needs to be fixed +nBeamElements = int(beamLength[0]/he*0.5) +nBeamElements=max(nBeamElements,3) +print nBeamElements diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/tank3D_so.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/tank3D_so.py new file mode 100644 index 00000000..f8c11c9b --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/tank3D_so.py @@ -0,0 +1,41 @@ +from proteus.default_so import * +import tank3D + +if tank3D.useOnlyVF: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n")] +else: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n"), + ("ls_p", "ls_n"), + ("redist_p", "redist_n"), + ("ls_consrv_p", "ls_consrv_n")] + + +if tank3D.useRANS > 0: + pnList.append(("kappa_p", + "kappa_n")) + pnList.append(("dissipation_p", + "dissipation_n")) +name = "tank3D_p" + +if tank3D.timeDiscretization == 'flcbdf': + systemStepControllerType = Sequential_MinFLCBDFModelStep + systemStepControllerType = Sequential_MinAdaptiveModelStep +else: + systemStepControllerType = Sequential_MinAdaptiveModelStep + +needEBQ_GLOBAL = False +needEBQ = False + +tnList = [0.0,tank3D.dt_init]+[i*tank3D.dt_fixed for i in range(1,tank3D.nDTout+1)] + +info = open("TimeList.txt","w") + + +for time in tnList: + info.write(str(time)+"\n") +info.close() + + +archiveFlag = ArchiveFlags.EVERY_SEQUENCE_STEP diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/twp_navier_stokes_n.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/twp_navier_stokes_n.py new file mode 100644 index 00000000..4e5ff8d1 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/twp_navier_stokes_n.py @@ -0,0 +1,69 @@ +from proteus import * +from twp_navier_stokes_p import * +from tank3D import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller_sys + stepController = Min_dt_cfl_controller + time_tol = 10.0*ns_nl_atol_res + atol_u = {1:time_tol,2:time_tol,3:time_tol} + rtol_u = {1:time_tol,2:time_tol,3:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis, + 1:basis, + 2:basis, + 3:basis} + +massLumping = False +numericalFluxType = None +conservativeFlux = None + +numericalFluxType = RANS2P.NumericalFlux +subgridError = RANS2P.SubgridError(coefficients,nd,lag=ns_lag_subgridError,hFactor=hFactor) +shockCapturing = RANS2P.ShockCapturing(coefficients,nd,ns_shockCapturingFactor,lag=ns_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None + +linearSmoother = SimpleNavierStokes3D + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rans2p_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*ns_nl_atol_res +nl_atol_res = ns_nl_atol_res +useEisenstatWalker = False +maxNonlinearIts = 50 +maxLineSearches = 0 +conservativeFlux = {0:'pwl-bdm-opt'} + + +#auxiliaryVariables=[pointGauges,rzWaveGenerator] diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/twp_navier_stokes_p.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/twp_navier_stokes_p.py new file mode 100644 index 00000000..9f5c7798 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/twp_navier_stokes_p.py @@ -0,0 +1,218 @@ +from proteus import * +from proteus.default_p import * +from tank3D import * +from proteus.mprans import RANS2P +from proteus.mprans import RANS2P_IB + +LevelModelType = RANS2P_IB.LevelModel +if useOnlyVF: + LS_model = None +else: + LS_model = 2 +if useRANS >= 1: + Closure_0_model = 5; Closure_1_model=6 + if useOnlyVF: + Closure_0_model=2; Closure_1_model=3 + if movingDomain: + Closure_0_model += 1; Closure_1_model += 1 +else: + Closure_0_model = None + Closure_1_model = None + +if spongeLayer or levee or slopingSpongeLayer: + coefficients = RANS2P_IB.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain, + porosityTypes=porosityTypes, + dragAlphaTypes=dragAlphaTypes, + dragBetaTypes=dragBetaTypes, + epsFact_solid = epsFact_solidTypes, + beamLocation=beamLocation, + beamLength=beamLength, + beamRadius=beamRadius, + EI=EI, + GJ=GJ, + nBeamElements=nBeamElements, + beam_quadOrder=beam_quadOrder, + beamFilename=beamFilename, + beam_Cd=1.2, + beam_nlTol=1.0e-5, + beamRigid=True) +else: + coefficients = RANS2P_IB.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain, + beamLocation=beamLocation, + beamLength=beamLength, + beamRadius=beamRadius, + EI=EI, + GJ=GJ, + nBeamElements=nBeamElements, + beam_quadOrder=beam_quadOrder, + beamFilename=beamFilename, + beam_Cd=1.2, + beam_nlTol=1.0e-5, + beamRigid=True) + + +def getDBC_p(x,flag): + if flag == boundaryTags['top']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return outflowPressure + +def getDBC_u(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_u + + +def getDBC_v(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_v + +def getDBC_w(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_w + +dirichletConditions = {0:getDBC_p, + 1:getDBC_u, + 2:getDBC_v, + 3:getDBC_w} + +def getAFBC_p(x,flag): + if flag == boundaryTags['left']: + return lambda x,t: -twpflowVelocity_u(x,t) + elif flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getAFBC_u(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getAFBC_v(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getAFBC_w(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getDFBC_u(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + + +def getDFBC_v(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + + +def getDFBC_w(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 +# if flag == boundaryTags['bottom']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['back']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['front']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['right']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['top']: +# return lambda x,t: 0.0 + + +advectiveFluxBoundaryConditions = {0:getAFBC_p, + 1:getAFBC_u, + 2:getAFBC_v, + 3:getAFBC_w} + +diffusiveFluxBoundaryConditions = {0:{}, + 1:{1:getDFBC_u}, + 2:{2:getDFBC_v}, + 3:{3:getDFBC_w}} + +class PerturbedSurface_p: + def __init__(self,waterLevel): + self.waterLevel=waterLevel + def uOfXT(self,x,t): + # if signedDistance(x) < 0: + # return -(L[2] - self.waterLevel)*rho_1*g[2] - (self.waterLevel - x[2])*rho_0*g[2] + # else: + # return -(L[2] - self.waterLevel)*rho_1*g[2] + if x[2]>inflowHeightMean: + return (L[2]-x[2])*rho_1*abs(g[2]) + else: + return (L[2]-inflowHeightMean)*rho_1*abs(g[2])+(inflowHeightMean-x[2])*rho_0*abs(g[2]) + +class AtRest: + def __init__(self): + pass + def uOfXT(self,x,t): + return 0.0 + +class WaterVelocity_u: + def __init__(self): + pass + def uOfXT(self,x,t): + + return waterVelocity(x,0) + +initialConditions = {0:PerturbedSurface_p(waterLine_z), + 1:AtRest(), + 2:AtRest(), + 3:AtRest()} diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/vegZoneVelocityInterp.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/vegZoneVelocityInterp.py new file mode 100644 index 00000000..6880b430 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/vegZoneVelocityInterp.py @@ -0,0 +1,43 @@ +import numpy as np +import scipy.interpolate as interp +import matplotlib.pyplot as plt + +vals = np.loadtxt("vegZoneVelocity.csv", + skiprows = 1, + delimiter=',') + +f = open("vegZoneVelocity.csv", 'r') +first = f.readline() +f.close() +g = first.split(',') +zListU=[] +zListW=[] +for i in range(1,len(g)): + h = g[i].split() + if h[0] == 'u': + zListU.append(float(h[3])) + elif h[0] == 'v': + zListW.append(float(h[3])) + +time = vals[:,0] +zU = np.array(zListU)-zListU[0] +zW = np.array(zListW)-zListW[0] + +Uarray = vals[:,1:len(zListU)+1] +Warray = vals[:,len(zListU)+1::] + +interpU = interp.RectBivariateSpline(time,zU,Uarray, kx=1,ky=3) +interpW = interp.RectBivariateSpline(time,zW, Warray, kx=1, ky=3) +#plt.contourf(tt,zzU,Uarray) +#plt.imshow(Uarray) +#plt.show() +vals = np.loadtxt("column_gauge.csv", + skiprows=1, + delimiter=',') + + +time = vals[:,0] +water_height = vals[:,4] + +interp_phi=interp.interp1d(time,water_height) + diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/vof_n.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/vof_n.py new file mode 100644 index 00000000..db0de251 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/vof_n.py @@ -0,0 +1,65 @@ +from proteus import * +from tank3D import * +from vof_p import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*vof_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = VOF.NumericalFlux +conservativeFlux = None +subgridError = VOF.SubgridError(coefficients=coefficients,nd=nd) +shockCapturing = VOF.ShockCapturing(coefficients,nd,shockCapturingFactor=vof_shockCapturingFactor,lag=vof_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'vof_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = vof_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*vof_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + +#auxiliaryVariables = [columnGauge] + diff --git a/3d/waves_vegetation/wavetank_vegetation_from_2d/vof_p.py b/3d/waves_vegetation/wavetank_vegetation_from_2d/vof_p.py new file mode 100644 index 00000000..ea33ff3b --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_from_2d/vof_p.py @@ -0,0 +1,47 @@ +from proteus import * +from proteus.default_p import * +from proteus.ctransportCoefficients import smoothedHeaviside +from tank3D import * +from proteus.mprans import VOF + +LevelModelType = VOF.LevelModel +if useOnlyVF: + RD_model = None + LS_model = None +else: + RD_model = 3 + LS_model = 2 + +coefficients = VOF.Coefficients(LS_model=LS_model,V_model=0,RD_model=RD_model,ME_model=1, + checkMass=False,useMetrics=useMetrics, + epsFact=epsFact_vof,sc_uref=vof_sc_uref,sc_beta=vof_sc_beta,movingDomain=movingDomain) + +def getDBC_vof(x,flag): + if flag == boundaryTags['left']: + return waveVF + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return lambda x,t: 1.0 +# elif flag == boundaryTags['right']: +# return outflowVF + + +dirichletConditions = {0:getDBC_vof} + +def getAFBC_vof(x,flag): + if flag == boundaryTags['left']: + return None + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return None +# elif flag == boundaryTags['right']: +# return None + else: + return lambda x,t: 0.0 + +advectiveFluxBoundaryConditions = {0:getAFBC_vof} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_H: + def uOfXT(self,x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,signedDistance(x)) + +initialConditions = {0:PerturbedSurface_H()} diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_consrv_n.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_consrv_n.py new file mode 100644 index 00000000..d6fdb0a9 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_consrv_n.py @@ -0,0 +1,48 @@ +from proteus import * +from tank3D import * +from ls_consrv_p import * + +timeIntegrator = ForwardIntegrator +timeIntegration = NoIntegration + +femSpaces = {0:basis} + +subgridError = None +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +shockCapturing = None + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'mcorr_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*mcorr_nl_atol_res +nl_atol_res = mcorr_nl_atol_res +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_consrv_p.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_consrv_p.py new file mode 100644 index 00000000..ceb09427 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_consrv_p.py @@ -0,0 +1,26 @@ +from proteus import * +from proteus.default_p import * +from tank3D import * +from proteus.mprans import MCorr + +LevelModelType = MCorr.LevelModel + +coefficients = MCorr.Coefficients(LSModel_index=2,V_model=0,me_model=4,VOFModel_index=1, + applyCorrection=applyCorrection,nd=nd,checkMass=False,useMetrics=useMetrics, + epsFactHeaviside=epsFact_consrv_heaviside, + epsFactDirac=epsFact_consrv_dirac, + epsFactDiffusion=epsFact_consrv_diffusion) + +class zero_phi: + def __init__(self): + pass + def uOfX(self,X): + return 0.0 + def uOfXT(self,X,t): + return 0.0 + +initialConditions = {0:zero_phi()} + + + + diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_n.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_n.py new file mode 100644 index 00000000..fde02e06 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_n.py @@ -0,0 +1,63 @@ +from proteus import * +from ls_p import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*ls_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +conservativeFlux = None +numericalFluxType = NCLS.NumericalFlux +subgridError = NCLS.SubgridError(coefficients,nd) +shockCapturing = NCLS.ShockCapturing(coefficients,nd,shockCapturingFactor=ls_shockCapturingFactor,lag=ls_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'ncls_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = ls_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*ls_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + + diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_p.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_p.py new file mode 100644 index 00000000..9d6237c1 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/ls_p.py @@ -0,0 +1,29 @@ +from proteus import * +from proteus.default_p import * +from tank3D import * +from proteus.mprans import NCLS + +LevelModelType = NCLS.LevelModel + +coefficients = NCLS.Coefficients(V_model=0,RD_model=3,ME_model=2, + checkMass=False, useMetrics=useMetrics, + epsFact=epsFact_consrv_heaviside,sc_uref=ls_sc_uref,sc_beta=ls_sc_beta,movingDomain=movingDomain) + +def getDBC_ls(x,flag): + if flag == boundaryTags['left']: + return wavePhi +# elif flag == boundaryTags['right']: +# return outflowPhi + else: + return None + +dirichletConditions = {0:getDBC_ls} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/redist_n.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/redist_n.py new file mode 100644 index 00000000..4f85688b --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/redist_n.py @@ -0,0 +1,66 @@ +from proteus import * +from redist_p import * +from tank3D import * + +nl_atol_res = rd_nl_atol_res +tolFac = 0.0 +nl_atol_res = rd_nl_atol_res + +linTolFac = 0.01 +l_atol_res = 0.01*rd_nl_atol_res + +if redist_Newton: + timeIntegration = NoIntegration + stepController = Newton_controller + maxNonlinearIts = 50 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'r' + levelNonlinearSolverConvergenceTest = 'r' + linearSolverConvergenceTest = 'r-true' + useEisenstatWalker = False +else: + timeIntegration = BackwardEuler_cfl + stepController = RDLS.PsiTC + runCFL=2.0 + psitc['nStepsForce']=3 + psitc['nStepsMax']=50 + psitc['reduceRatio']=2.0 + psitc['startRatio']=1.0 + rtol_res[0] = 0.0 + atol_res[0] = rd_nl_atol_res + useEisenstatWalker = False + maxNonlinearIts = 1 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'rits' + levelNonlinearSolverConvergenceTest = 'rits' + linearSolverConvergenceTest = 'r-true' + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +subgridError = RDLS.SubgridError(coefficients,nd) +shockCapturing = RDLS.ShockCapturing(coefficients,nd,shockCapturingFactor=rd_shockCapturingFactor,lag=rd_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = NLGaussSeidel +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rdls_' diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/redist_p.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/redist_p.py new file mode 100644 index 00000000..ea199426 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/redist_p.py @@ -0,0 +1,32 @@ +from proteus import * +from proteus.default_p import * +from math import * +from tank3D import * +from proteus.mprans import RDLS +""" +The redistancing equation in the sloshbox test problem. +""" + +LevelModelType = RDLS.LevelModel + +coefficients = RDLS.Coefficients(applyRedistancing=applyRedistancing, + epsFact=epsFact_redistance, + nModelId=2, + rdModelId=3, + useMetrics=useMetrics, + backgroundDiffusionFactor=backgroundDiffusionFactor) + +def getDBC_rd(x,flag): + pass + +dirichletConditions = {0:getDBC_rd} +weakDirichletConditions = {0:RDLS.setZeroLSweakDirichletBCsSimple} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/tank3D.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/tank3D.py new file mode 100644 index 00000000..8fe42d32 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/tank3D.py @@ -0,0 +1,652 @@ +from math import * +import proteus.MeshTools +from proteus import Domain +from proteus.default_n import * +from proteus.Profiling import logEvent +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.ctransportCoefficients import smoothedHeaviside_integral +from proteus import Gauges +from proteus.Gauges import PointGauges,LineGauges,LineIntegralGauges +from proteus import WaveTools as WT + +from proteus import Context +opts=Context.Options([ + ("wave_type", 'linear', "type of waves generated: 'linear', 'Nonlinear', 'single-peaked', 'double-peaked'"), + ("depth", 0.457, "water depth at leading edge of vegetation (not at wave generator)[m]"), + ("wave_height", 0.192, "wave height at leading edget of vegetation [m]"), + ("peak_period", 2.0, "Peak period [s]"), + ("peak_period2", 1.5, "Second peak period (only used in double-peaked case)[s]"), + ("peak_wavelength",3.91,"Peak wavelength in [m]"), + ("parallel", False, "Run in parallel")]) + +#wave generator +windVelocity = (0.0,0.0,0.0) +#veg_platform_height = 17.2/44.0 + 6.1/20.0 +depth = opts.depth +inflowHeightMean = depth +inflowVelocityMean = (0.0,0.0,0,0) +period = opts.peak_period +omega = 2.0*math.pi/period +waveheight = opts.wave_height +amplitude = waveheight/ 2.0 +wavelength = opts.peak_wavelength +k = 2.0*math.pi/wavelength +waveDir = numpy.array([1,0,0]) +g = numpy.array([0.0,0.0,-9.81]) +if opts.wave_type == 'linear': + waves = WT.MonochromaticWaves(period = period, # Peak period + waveHeight = waveheight, # Height + depth = depth, # Depth + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + waveType="Linear") +elif opts.wave_type == 'Nonlinear': + waves = WT.MonochromaticWaves(period = period, # Peak period + waveHeight = waveheight, # Height + wavelength = wavelength, + depth = depth, # Depth + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + waveType="Fenton", + Ycoeff = [0.04160592, #Surface elevation Fourier coefficients for non-dimensionalised solution + 0.00555874, + 0.00065892, + 0.00008144, + 0.00001078, + 0.00000151, + 0.00000023, + 0.00000007], + Bcoeff = [0.05395079, + 0.00357780, + 0.00020506, + 0.00000719, + -0.00000016, + -0.00000005, + 0.00000000, + 0.00000000]) +elif opts.wave_type == 'single-peaked': + waves = WT.RandomWaves( Tp = period, # Peak period + Hs = waveheight, # Height + d = depth, # Depth + fp = 1./period, #peak Frequency + bandFactor = 2.0, #fmin=fp/Bandfactor, fmax = Bandfactor * fp + N = 101, #No of frequencies for signal reconstruction + mwl = inflowHeightMean, # Sea water level + waveDir = waveDir, # waveDirection + g = g, # Gravity vector, defines the vertical + gamma=3.3, + spec_fun = WT.JONSWAP) +# elif opts.wave_type == 'double-peaked': +# waves = WT.DoublePeakedRandomWaves( Tp = period, # Peak period +# Hs = waveheight, # Height +# d = depth, # Depth +# fp = 1./period, #peak Frequency +# bandFactor = 2.0, #fmin=fp/Bandfactor, fmax = Bandfactor * fp +# N = 101, #No of frequencies for signal reconstruction +# Hs = 2.0, #m significant wave height +# d = 2.0, #m depth +# fp = 1.0/5.0, #peak frequency +# bandFactor = 2.0, #controls width of band around fp +# N = 101, #number of frequency bins +# mwl = 0.0, #mean water level +# waveDir = np.array([1,0,0]), +# g = np.array([0, -9.81, 0]), #accelerationof gravity +# spec_fun = JONSWAP, +# gamma=3.3, mwl = inflowHeightMean, # Sea water level +# waveDir = waveDir, # waveDirection +# g = g, # Gravity vector, defines the vertical +# gamma=10.0, +# spec_fun = WT.JONSWAP, +# Tp_2 = opts.peak_period2) + +# Discretization -- input options +genMesh=True +movingDomain=False +applyRedistancing=True +useOldPETSc=False +useSuperlu=not opts.parallel +timeDiscretization='be'#'be','vbdf','flcbdf' +spaceOrder = 1 +useHex = False +useRBLES = 0.0 +useMetrics = 1.0 +applyCorrection=True +useVF = 1.0 +useOnlyVF = False +useRANS = 0 # 0 -- None + # 1 -- K-Epsilon + # 2 -- K-Omega +# Input checks +if spaceOrder not in [1,2]: + print "INVALID: spaceOrder" + spaceOrder + sys.exit() + +if useRBLES not in [0.0, 1.0]: + print "INVALID: useRBLES" + useRBLES + sys.exit() + +if useMetrics not in [0.0, 1.0]: + print "INVALID: useMetrics" + sys.exit() + +# Discretization +nd = 3 +if spaceOrder == 1: + hFactor=1.0 + if useHex: + basis=C0_AffineLinearOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,2) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,2) + else: + basis=C0_AffineLinearOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,3) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,3) +elif spaceOrder == 2: + hFactor=0.5 + if useHex: + basis=C0_AffineLagrangeOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,4) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,4) + else: + basis=C0_AffineQuadraticOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,4) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,4) + +# Domain and mesh + +#for debugging, make the tank short + +he = float(wavelength)/30.0#0.0#100 +L = (15.0, 1.5, 1.0) + +GenerationZoneLength = 1.2 +AbsorptionZoneLength= 2.8 +spongeLayer = True +xSponge = GenerationZoneLength +xRelaxCenter = xSponge/2.0 +epsFact_solid = xSponge/2.0 +#zone 2 +xSponge_2 = L[0] - AbsorptionZoneLength +xRelaxCenter_2 = 0.5*(xSponge_2+L[0]) +epsFact_solid_2 = AbsorptionZoneLength/2.0 + +weak_bc_penalty_constant = 100.0 +nLevels = 1 +#parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.element +parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.node +nLayersOfOverlapForParallel = 0 +structured=False + +if useHex: + nnx=4*Refinement+1 + nny=2*Refinement+1 + hex=True + domain = Domain.RectangularDomain(L) +else: + boundaries=['empty','left','right','bottom','top','front','back'] + boundaryTags=dict([(key,i+1) for (i,key) in enumerate(boundaries)]) + if structured: + nnx=4*Refinement + nny=2*Refinement + domain = Domain.RectangularDomain(L) + elif spongeLayer: + vertices=[[0.0,0.0,0.0],#0 + [xSponge,0.0,0.0],#1 + [xSponge_2,0.0,0.0],#2 + [L[0],0.0,0.0],#3 + [L[0],L[1],0.0],#4 + [xSponge_2,L[1],0.0],#5 + [xSponge,L[1],0.0],#6 + [0.0,L[1],0.0]]#7 + + + vertexFlags=[boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom']] + + + for v,vf in zip(vertices,vertexFlags): + vertices.append([v[0],v[1],L[2]]) + vertexFlags.append(boundaryTags['top']) + + print vertices + print vertexFlags + + segments=[[0,1], + [1,2], + [2,3], + [3,4], + [4,5], + [5,6], + [6,7], + [7,0], + [1,6], + [2,5]] + + segmentFlags=[boundaryTags['front'], + boundaryTags['front'], + boundaryTags['front'], + boundaryTags['right'], + boundaryTags['back'], + boundaryTags['back'], + boundaryTags['back'], + boundaryTags['left'], + boundaryTags['empty'], + boundaryTags['empty'] ] + + + facets=[] + facetFlags=[] + + for s,sF in zip(segments,segmentFlags): + facets.append([[s[0],s[1],s[1]+8,s[0]+8]]) + facetFlags.append(sF) + + bf=[[0,1,6,7],[1,2,5,6],[2,3,4,5]] + tf=[] + for i in range(0,3): + facets.append([bf[i]]) + tf=[ss + 8 for ss in bf[i]] + facets.append([tf]) + + for i in range(0,3): + facetFlags.append(boundaryTags['bottom']) + facetFlags.append(boundaryTags['top']) + + print facets + print facetFlags + + regions=[[xRelaxCenter, 0.5*L[1],0.0], + [xRelaxCenter_2, 0.5*L[1], 0.0], + [0.5*L[0],0.1*L[1], 0.0]] + regionFlags=[1,2,3] + + domain = Domain.PiecewiseLinearComplexDomain(vertices=vertices, + vertexFlags=vertexFlags, + facets=facets, + facetFlags=facetFlags, + regions=regions, + regionFlags=regionFlags, + ) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="KVApq1.4q12feena%21.16e" % ((he**3)/6.0,) + + + logEvent("""Mesh generated using: tetgen -%s %s""" % (triangleOptions,domain.polyfile+".poly")) + + porosityTypes = numpy.array([1.0, + 1.0, + 1.0, + 1.0]) + dragAlphaTypes = numpy.array([0.0, + 0.0, + 0.0, + 0.5/1.004e-6]) + + dragBetaTypes = numpy.array([0.0,0.0,0.0,0.0]) + + epsFact_solidTypes = np.array([0.0,0.0,0.0,0.0]) + + else: + vertices=[[0.0,0.0,0.0],#0 + [L[0],0.0,0.0],#1 + [L[0],L[1],0.0],#2 + [0.0,L[1],0.0]]#3 + + + vertexFlags=[boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom']] + + + for v,vf in zip(vertices,vertexFlags): + vertices.append([v[0],v[1],L[2]]) + vertexFlags.append(boundaryTags['top']) + + segments=[[0,1], + [1,2], + [2,3], + [3,0]] + + segmentFlags=[boundaryTags['front'], + boundaryTags['right'], + boundaryTags['back'], + boundaryTags['left']] + + facets=[] + facetFlags=[] + + for s,sF in zip(segments,segmentFlags): + facets.append([[s[0],s[1],s[1]+4,s[0]+4]]) + facetFlags.append(sF) + + bf=[[0,1,2,3]] + tf=[] + for i in range(0,1): + facets.append([bf[i]]) + tf=[ss + 4 for ss in bf[i]] + facets.append([tf]) + + for i in range(0,1): + facetFlags.append(boundaryTags['bottom']) + facetFlags.append(boundaryTags['top']) + + for s,sF in zip(segments,segmentFlags): + segments.append([s[1]+4,s[0]+4]) + segmentFlags.append(sF) + + + regions=[[0.5*L[0],0.5*L[1], 0.0]] + regionFlags=[1] + + domain = Domain.PiecewiseLinearComplexDomain(vertices=vertices, + vertexFlags=vertexFlags, + facets=facets, + facetFlags=facetFlags, + regions=regions, + regionFlags=regionFlags) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="KVApq1.4q12feena%21.16e" % ((he**3)/6.0,) + + + logEvent("""Mesh generated using: tetgen -%s %s""" % (triangleOptions,domain.polyfile+".poly")) + + + +# Time stepping +T=40*period +dt_fixed = period/11.0#2.0*0.5/20.0#T/2.0#period/21.0 +dt_init = min(0.001*dt_fixed,0.001) +runCFL=0.90 +nDTout = int(round(T/dt_fixed)) + +# Numerical parameters +ns_forceStrongDirichlet = False#True +backgroundDiffusionFactor=0.0 +if useMetrics: + ns_shockCapturingFactor = 0.25 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.25 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.25 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.25 + rd_lag_shockCapturing = False + epsFact_density = 3.0 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 1.5 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.1 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.1 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 +else: + ns_shockCapturingFactor = 0.9 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.9 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.9 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.9 + rd_lag_shockCapturing = False + epsFact_density = 1.5 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 0.33 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.9 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.9 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 + +ns_nl_atol_res = max(1.0e-10,0.0001*he**2) +vof_nl_atol_res = max(1.0e-10,0.0001*he**2) +ls_nl_atol_res = max(1.0e-10,0.0001*he**2) +rd_nl_atol_res = max(1.0e-10,0.005*he) +mcorr_nl_atol_res = max(1.0e-10,0.0001*he**2) +kappa_nl_atol_res = max(1.0e-10,0.001*he**2) +dissipation_nl_atol_res = max(1.0e-10,0.001*he**2) + +#turbulence +ns_closure=0 #1-classic smagorinsky, 2-dynamic smagorinsky, 3 -- k-epsilon, 4 -- k-omega +if useRANS == 1: + ns_closure = 3 +elif useRANS == 2: + ns_closure == 4 +# Water +rho_0 = 998.2 +nu_0 = 1.004e-6 + +# Air +rho_1 = 1.205 +nu_1 = 1.500e-5 + +# Surface tension +sigma_01 = 0.0 + +# Gravity +g = [0.0,0.0,-9.8] + +# Initial condition +waterLine_x = 2*L[0] +waterLine_z = inflowHeightMean + + +def signedDistance(x): + phi_x = x[0]-waterLine_x + phi_z = x[2]-waterLine_z + if phi_x < 0.0: + if phi_z < 0.0: + return max(phi_x,phi_z) + else: + return phi_z + else: + if phi_z < 0.0: + return phi_x + else: + return sqrt(phi_x**2 + phi_z**2) + + +def theta(x,t): + return k*x[0] - omega*t + pi/2.0 + +def z(x): + return x[2] - inflowHeightMean + +#sigma = omega - k*inflowVelocityMean[0] +h = inflowHeightMean # - transect[0][1] if lower left hand corner is not at z=0 + +#waveData + +def waveHeight(x,t): + return inflowHeightMean + waves.eta(x[0],x[1],x[2],t) +def waveVelocity_u(x,t): + return waves.u(x[0],x[1],x[2],t,"x") +def waveVelocity_v(x,t): + return waves.u(x[0],x[1],x[2],t,"y") +def waveVelocity_w(x,t): + return waves.u(x[0],x[1],x[2],t,"z") + +#solution variables + +def wavePhi(x,t): + return x[2] - waveHeight(x,t) + +def waveVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)) + +def twpflowVelocity_u(x,t): + waterspeed = waveVelocity_u(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + u = H*windVelocity[0] + (1.0-H)*waterspeed + return u + +def twpflowVelocity_v(x,t): + waterspeed = waveVelocity_v(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + return H*windVelocity[1]+(1.0-H)*waterspeed + +def twpflowFlux(x,t): + return -twpflowVelocity_u(x,t) + +outflowHeight=inflowHeightMean + +def outflowVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,x[2] - outflowHeight) + +def outflowPhi(x,t): + return x[2] - outflowHeight + +def outflowPressure(x,t): + if x[2]>inflowHeightMean: + return (L[2]-x[2])*rho_1*abs(g[2]) + else: + return (L[2]-inflowHeightMean)*rho_1*abs(g[2])+(inflowHeightMean-x[2])*rho_0*abs(g[2]) + + + #p_L = L[1]*rho_1*g[1] + #phi_L = L[1] - outflowHeight + #phi = x[1] - outflowHeight + #return p_L -g[1]*(rho_0*(phi_L - phi)+(rho_1 -rho_0)*(smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi_L) + # -smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi))) + +def twpflowVelocity_w(x,t): + return 0.0 + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + +RelaxationZone = namedtuple("RelaxationZone","center_x sign u v w") + +class RelaxationZoneWaveGenerator(AV_base): + """ Prescribe a velocity penalty scaling in a material zone via a Darcy-Forchheimer penalty + + :param zones: A dictionary mapping integer material types to Zones, where a Zone is a named tuple + specifying the x coordinate of the zone center and the velocity components + """ + def __init__(self,zones): + assert isinstance(zones,dict) + self.zones = zones + def calculate(self): + for l,m in enumerate(self.model.levelModelList): + for eN in range(m.coefficients.q_phi.shape[0]): + mType = m.mesh.elementMaterialTypes[eN] + if self.zones.has_key(mType): + for k in range(m.coefficients.q_phi.shape[1]): + t = m.timeIntegration.t + x = m.q['x'][eN,k] + m.coefficients.q_phi_solid[eN,k] = self.zones[mType].sign*(self.zones[mType].center_x - x[0]) + m.coefficients.q_velocity_solid[eN,k,0] = self.zones[mType].u(x,t) + m.coefficients.q_velocity_solid[eN,k,1] = self.zones[mType].v(x,t) + #m.coefficients.q_velocity_solid[eN,k,2] = self.zones[mType].w(x,t) + m.q['phi_solid'] = m.coefficients.q_phi_solid + m.q['velocity_solid'] = m.coefficients.q_velocity_solid + +rzWaveGenerator = RelaxationZoneWaveGenerator(zones={ + # 1:RelaxationZone(xRelaxCenter, + # 1.0, + # twpflowVelocity_u, + # twpflowVelocity_v, + # twpflowVelocity_w), + 1:RelaxationZone(xRelaxCenter_2, + -1.0, #currently Hs=1-exp_function + zeroVel, + zeroVel, + zeroVel)}) + +beam_quadOrder=3 +beam_useSparse=False +beamFilename="wavetankBeams" +#nBeamElements=max(nBeamElements,3) + +#beam info +beamLocation=[] +beamLength=[] +beamRadius=[] +EI=[] +GJ=[] +lam = 0.25 #0.05 #3.0*2.54/100.0 #57.4e-3 +lamx = 3.0**0.5*lam +xs = 1.2 +ys = 0.0 +xList=[] +yList = [] +while xs <= 11.0: + xList.append(xs) + xs += lam +while ys<= L[1]: + yList.append(ys) + ys+=lamx +for i in xList: + for j in yList: + beamLocation.append((i,j)) + beamLength.append(0.415) + beamRadius.append(0.0032) + EI.append(3.0e-4) # needs to be fixed + GJ.append(1.5e-4) # needs to be fixed + +xs = 1.2+0.5*lam +ys = 0.5*lamx +xList=[] +yList = [] +while xs <= 11.0: + xList.append(xs) + xs += lam + +while ys<= L[1]: + yList.append(ys) + ys+=lamx + +for i in xList: + for j in yList: + beamLocation.append((i,j)) + beamLength.append(0.415) + beamRadius.append(0.0032) + EI.append(3.0e-4) # needs to be fixed + GJ.append(1.5e-4) # needs to be fixed +nBeamElements = int(beamLength[0]/he*0.5) +nBeamElements=max(nBeamElements,3) +print nBeamElements diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/tank3D_so.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/tank3D_so.py new file mode 100644 index 00000000..f8c11c9b --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/tank3D_so.py @@ -0,0 +1,41 @@ +from proteus.default_so import * +import tank3D + +if tank3D.useOnlyVF: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n")] +else: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n"), + ("ls_p", "ls_n"), + ("redist_p", "redist_n"), + ("ls_consrv_p", "ls_consrv_n")] + + +if tank3D.useRANS > 0: + pnList.append(("kappa_p", + "kappa_n")) + pnList.append(("dissipation_p", + "dissipation_n")) +name = "tank3D_p" + +if tank3D.timeDiscretization == 'flcbdf': + systemStepControllerType = Sequential_MinFLCBDFModelStep + systemStepControllerType = Sequential_MinAdaptiveModelStep +else: + systemStepControllerType = Sequential_MinAdaptiveModelStep + +needEBQ_GLOBAL = False +needEBQ = False + +tnList = [0.0,tank3D.dt_init]+[i*tank3D.dt_fixed for i in range(1,tank3D.nDTout+1)] + +info = open("TimeList.txt","w") + + +for time in tnList: + info.write(str(time)+"\n") +info.close() + + +archiveFlag = ArchiveFlags.EVERY_SEQUENCE_STEP diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/twp_navier_stokes_n.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/twp_navier_stokes_n.py new file mode 100644 index 00000000..4e5ff8d1 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/twp_navier_stokes_n.py @@ -0,0 +1,69 @@ +from proteus import * +from twp_navier_stokes_p import * +from tank3D import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller_sys + stepController = Min_dt_cfl_controller + time_tol = 10.0*ns_nl_atol_res + atol_u = {1:time_tol,2:time_tol,3:time_tol} + rtol_u = {1:time_tol,2:time_tol,3:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis, + 1:basis, + 2:basis, + 3:basis} + +massLumping = False +numericalFluxType = None +conservativeFlux = None + +numericalFluxType = RANS2P.NumericalFlux +subgridError = RANS2P.SubgridError(coefficients,nd,lag=ns_lag_subgridError,hFactor=hFactor) +shockCapturing = RANS2P.ShockCapturing(coefficients,nd,ns_shockCapturingFactor,lag=ns_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None + +linearSmoother = SimpleNavierStokes3D + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rans2p_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*ns_nl_atol_res +nl_atol_res = ns_nl_atol_res +useEisenstatWalker = False +maxNonlinearIts = 50 +maxLineSearches = 0 +conservativeFlux = {0:'pwl-bdm-opt'} + + +#auxiliaryVariables=[pointGauges,rzWaveGenerator] diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/twp_navier_stokes_p.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/twp_navier_stokes_p.py new file mode 100644 index 00000000..9f5c7798 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/twp_navier_stokes_p.py @@ -0,0 +1,218 @@ +from proteus import * +from proteus.default_p import * +from tank3D import * +from proteus.mprans import RANS2P +from proteus.mprans import RANS2P_IB + +LevelModelType = RANS2P_IB.LevelModel +if useOnlyVF: + LS_model = None +else: + LS_model = 2 +if useRANS >= 1: + Closure_0_model = 5; Closure_1_model=6 + if useOnlyVF: + Closure_0_model=2; Closure_1_model=3 + if movingDomain: + Closure_0_model += 1; Closure_1_model += 1 +else: + Closure_0_model = None + Closure_1_model = None + +if spongeLayer or levee or slopingSpongeLayer: + coefficients = RANS2P_IB.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain, + porosityTypes=porosityTypes, + dragAlphaTypes=dragAlphaTypes, + dragBetaTypes=dragBetaTypes, + epsFact_solid = epsFact_solidTypes, + beamLocation=beamLocation, + beamLength=beamLength, + beamRadius=beamRadius, + EI=EI, + GJ=GJ, + nBeamElements=nBeamElements, + beam_quadOrder=beam_quadOrder, + beamFilename=beamFilename, + beam_Cd=1.2, + beam_nlTol=1.0e-5, + beamRigid=True) +else: + coefficients = RANS2P_IB.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain, + beamLocation=beamLocation, + beamLength=beamLength, + beamRadius=beamRadius, + EI=EI, + GJ=GJ, + nBeamElements=nBeamElements, + beam_quadOrder=beam_quadOrder, + beamFilename=beamFilename, + beam_Cd=1.2, + beam_nlTol=1.0e-5, + beamRigid=True) + + +def getDBC_p(x,flag): + if flag == boundaryTags['top']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return outflowPressure + +def getDBC_u(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_u + + +def getDBC_v(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_v + +def getDBC_w(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_w + +dirichletConditions = {0:getDBC_p, + 1:getDBC_u, + 2:getDBC_v, + 3:getDBC_w} + +def getAFBC_p(x,flag): + if flag == boundaryTags['left']: + return lambda x,t: -twpflowVelocity_u(x,t) + elif flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getAFBC_u(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getAFBC_v(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getAFBC_w(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getDFBC_u(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + + +def getDFBC_v(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + + +def getDFBC_w(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 +# if flag == boundaryTags['bottom']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['back']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['front']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['right']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['top']: +# return lambda x,t: 0.0 + + +advectiveFluxBoundaryConditions = {0:getAFBC_p, + 1:getAFBC_u, + 2:getAFBC_v, + 3:getAFBC_w} + +diffusiveFluxBoundaryConditions = {0:{}, + 1:{1:getDFBC_u}, + 2:{2:getDFBC_v}, + 3:{3:getDFBC_w}} + +class PerturbedSurface_p: + def __init__(self,waterLevel): + self.waterLevel=waterLevel + def uOfXT(self,x,t): + # if signedDistance(x) < 0: + # return -(L[2] - self.waterLevel)*rho_1*g[2] - (self.waterLevel - x[2])*rho_0*g[2] + # else: + # return -(L[2] - self.waterLevel)*rho_1*g[2] + if x[2]>inflowHeightMean: + return (L[2]-x[2])*rho_1*abs(g[2]) + else: + return (L[2]-inflowHeightMean)*rho_1*abs(g[2])+(inflowHeightMean-x[2])*rho_0*abs(g[2]) + +class AtRest: + def __init__(self): + pass + def uOfXT(self,x,t): + return 0.0 + +class WaterVelocity_u: + def __init__(self): + pass + def uOfXT(self,x,t): + + return waterVelocity(x,0) + +initialConditions = {0:PerturbedSurface_p(waterLine_z), + 1:AtRest(), + 2:AtRest(), + 3:AtRest()} diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/vof_n.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/vof_n.py new file mode 100644 index 00000000..db0de251 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/vof_n.py @@ -0,0 +1,65 @@ +from proteus import * +from tank3D import * +from vof_p import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*vof_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = VOF.NumericalFlux +conservativeFlux = None +subgridError = VOF.SubgridError(coefficients=coefficients,nd=nd) +shockCapturing = VOF.ShockCapturing(coefficients,nd,shockCapturingFactor=vof_shockCapturingFactor,lag=vof_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'vof_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = vof_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*vof_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + +#auxiliaryVariables = [columnGauge] + diff --git a/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/vof_p.py b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/vof_p.py new file mode 100644 index 00000000..ea33ff3b --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_generate_from_wavetools/vof_p.py @@ -0,0 +1,47 @@ +from proteus import * +from proteus.default_p import * +from proteus.ctransportCoefficients import smoothedHeaviside +from tank3D import * +from proteus.mprans import VOF + +LevelModelType = VOF.LevelModel +if useOnlyVF: + RD_model = None + LS_model = None +else: + RD_model = 3 + LS_model = 2 + +coefficients = VOF.Coefficients(LS_model=LS_model,V_model=0,RD_model=RD_model,ME_model=1, + checkMass=False,useMetrics=useMetrics, + epsFact=epsFact_vof,sc_uref=vof_sc_uref,sc_beta=vof_sc_beta,movingDomain=movingDomain) + +def getDBC_vof(x,flag): + if flag == boundaryTags['left']: + return waveVF + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return lambda x,t: 1.0 +# elif flag == boundaryTags['right']: +# return outflowVF + + +dirichletConditions = {0:getDBC_vof} + +def getAFBC_vof(x,flag): + if flag == boundaryTags['left']: + return None + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return None +# elif flag == boundaryTags['right']: +# return None + else: + return lambda x,t: 0.0 + +advectiveFluxBoundaryConditions = {0:getAFBC_vof} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_H: + def uOfXT(self,x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,signedDistance(x)) + +initialConditions = {0:PerturbedSurface_H()} diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/ls_consrv_n.py b/3d/waves_vegetation/wavetank_vegetation_time_series/ls_consrv_n.py new file mode 100644 index 00000000..d6fdb0a9 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/ls_consrv_n.py @@ -0,0 +1,48 @@ +from proteus import * +from tank3D import * +from ls_consrv_p import * + +timeIntegrator = ForwardIntegrator +timeIntegration = NoIntegration + +femSpaces = {0:basis} + +subgridError = None +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +shockCapturing = None + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'mcorr_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*mcorr_nl_atol_res +nl_atol_res = mcorr_nl_atol_res +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/ls_consrv_p.py b/3d/waves_vegetation/wavetank_vegetation_time_series/ls_consrv_p.py new file mode 100644 index 00000000..ceb09427 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/ls_consrv_p.py @@ -0,0 +1,26 @@ +from proteus import * +from proteus.default_p import * +from tank3D import * +from proteus.mprans import MCorr + +LevelModelType = MCorr.LevelModel + +coefficients = MCorr.Coefficients(LSModel_index=2,V_model=0,me_model=4,VOFModel_index=1, + applyCorrection=applyCorrection,nd=nd,checkMass=False,useMetrics=useMetrics, + epsFactHeaviside=epsFact_consrv_heaviside, + epsFactDirac=epsFact_consrv_dirac, + epsFactDiffusion=epsFact_consrv_diffusion) + +class zero_phi: + def __init__(self): + pass + def uOfX(self,X): + return 0.0 + def uOfXT(self,X,t): + return 0.0 + +initialConditions = {0:zero_phi()} + + + + diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/ls_n.py b/3d/waves_vegetation/wavetank_vegetation_time_series/ls_n.py new file mode 100644 index 00000000..fde02e06 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/ls_n.py @@ -0,0 +1,63 @@ +from proteus import * +from ls_p import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*ls_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +conservativeFlux = None +numericalFluxType = NCLS.NumericalFlux +subgridError = NCLS.SubgridError(coefficients,nd) +shockCapturing = NCLS.ShockCapturing(coefficients,nd,shockCapturingFactor=ls_shockCapturingFactor,lag=ls_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'ncls_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = ls_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*ls_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + + diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/ls_p.py b/3d/waves_vegetation/wavetank_vegetation_time_series/ls_p.py new file mode 100644 index 00000000..9d6237c1 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/ls_p.py @@ -0,0 +1,29 @@ +from proteus import * +from proteus.default_p import * +from tank3D import * +from proteus.mprans import NCLS + +LevelModelType = NCLS.LevelModel + +coefficients = NCLS.Coefficients(V_model=0,RD_model=3,ME_model=2, + checkMass=False, useMetrics=useMetrics, + epsFact=epsFact_consrv_heaviside,sc_uref=ls_sc_uref,sc_beta=ls_sc_beta,movingDomain=movingDomain) + +def getDBC_ls(x,flag): + if flag == boundaryTags['left']: + return wavePhi +# elif flag == boundaryTags['right']: +# return outflowPhi + else: + return None + +dirichletConditions = {0:getDBC_ls} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/redist_n.py b/3d/waves_vegetation/wavetank_vegetation_time_series/redist_n.py new file mode 100644 index 00000000..4f85688b --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/redist_n.py @@ -0,0 +1,66 @@ +from proteus import * +from redist_p import * +from tank3D import * + +nl_atol_res = rd_nl_atol_res +tolFac = 0.0 +nl_atol_res = rd_nl_atol_res + +linTolFac = 0.01 +l_atol_res = 0.01*rd_nl_atol_res + +if redist_Newton: + timeIntegration = NoIntegration + stepController = Newton_controller + maxNonlinearIts = 50 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'r' + levelNonlinearSolverConvergenceTest = 'r' + linearSolverConvergenceTest = 'r-true' + useEisenstatWalker = False +else: + timeIntegration = BackwardEuler_cfl + stepController = RDLS.PsiTC + runCFL=2.0 + psitc['nStepsForce']=3 + psitc['nStepsMax']=50 + psitc['reduceRatio']=2.0 + psitc['startRatio']=1.0 + rtol_res[0] = 0.0 + atol_res[0] = rd_nl_atol_res + useEisenstatWalker = False + maxNonlinearIts = 1 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'rits' + levelNonlinearSolverConvergenceTest = 'rits' + linearSolverConvergenceTest = 'r-true' + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = DoNothing +conservativeFlux = None +subgridError = RDLS.SubgridError(coefficients,nd) +shockCapturing = RDLS.ShockCapturing(coefficients,nd,shockCapturingFactor=rd_shockCapturingFactor,lag=rd_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = NLGaussSeidel +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rdls_' diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/redist_p.py b/3d/waves_vegetation/wavetank_vegetation_time_series/redist_p.py new file mode 100644 index 00000000..ea199426 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/redist_p.py @@ -0,0 +1,32 @@ +from proteus import * +from proteus.default_p import * +from math import * +from tank3D import * +from proteus.mprans import RDLS +""" +The redistancing equation in the sloshbox test problem. +""" + +LevelModelType = RDLS.LevelModel + +coefficients = RDLS.Coefficients(applyRedistancing=applyRedistancing, + epsFact=epsFact_redistance, + nModelId=2, + rdModelId=3, + useMetrics=useMetrics, + backgroundDiffusionFactor=backgroundDiffusionFactor) + +def getDBC_rd(x,flag): + pass + +dirichletConditions = {0:getDBC_rd} +weakDirichletConditions = {0:RDLS.setZeroLSweakDirichletBCsSimple} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_phi: + def uOfXT(self,x,t): + return signedDistance(x) + +initialConditions = {0:PerturbedSurface_phi()} diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/tank3D.py b/3d/waves_vegetation/wavetank_vegetation_time_series/tank3D.py new file mode 100644 index 00000000..c744e92f --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/tank3D.py @@ -0,0 +1,600 @@ +from math import * +import proteus.MeshTools +from proteus import Domain +from proteus.default_n import * +from proteus.Profiling import logEvent +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.ctransportCoefficients import smoothedHeaviside_integral +from proteus import Gauges +from proteus.Gauges import PointGauges,LineGauges,LineIntegralGauges +from proteus.WaveTools import timeSeries + + + +from proteus import Context +opts=Context.Options([ + ("wave_type", 'linear', "type of waves generated: 'linear', 'Nonlinear', 'single-peaked', 'double-peaked'"), + ("depth", 0.457, "water depth at leading edge of vegetation (not at wave generator)[m]"), + ("wave_height", 0.192, "wave height at leading edget of vegetation [m]"), + ("peak_period", 2.0, "Peak period [s]"), + ("peak_period2", 1.5, "Second peak period (only used in double-peaked case)[s]"), + ("peak_wavelength",3.91,"Peak wavelength in [m]"), + ("parallel", False, "Run in parallel")]) + +#wave generator +windVelocity = (0.0,0.0,0.0) +#veg_platform_height = 17.2/44.0 + 6.1/20.0 +depth = opts.depth +inflowHeightMean = depth +inflowVelocityMean = (0.0,0.0,0,0) +period = opts.peak_period +omega = 2.0*math.pi/period +waveheight = opts.wave_height +amplitude = waveheight/ 2.0 +wavelength = opts.peak_wavelength +k = 2.0*math.pi/wavelength +waveDir = numpy.array([1,0,0]) +g = numpy.array([0.0,0.0,-9.81]) + +tS = timeSeries(timeSeriesFile = "monochromatic_linear_time_series.txt", + d = depth, + Npeaks = 1, #m depth + bandFactor = [2.0], #controls width of band around fp + peakFrequencies = [1.0], + N = 32, #number of frequency bins + Nwaves = 20, + mwl = inflowHeightMean, #mean water level + waveDir = np.array([1,0,0]), + g = np.array(g)) + + + +# Discretization -- input options +genMesh=True +movingDomain=False +applyRedistancing=True +useOldPETSc=False +useSuperlu=not opts.parallel +timeDiscretization='be'#'be','vbdf','flcbdf' +spaceOrder = 1 +useHex = False +useRBLES = 0.0 +useMetrics = 1.0 +applyCorrection=True +useVF = 1.0 +useOnlyVF = False +useRANS = 0 # 0 -- None + # 1 -- K-Epsilon + # 2 -- K-Omega +# Input checks +if spaceOrder not in [1,2]: + print "INVALID: spaceOrder" + spaceOrder + sys.exit() + +if useRBLES not in [0.0, 1.0]: + print "INVALID: useRBLES" + useRBLES + sys.exit() + +if useMetrics not in [0.0, 1.0]: + print "INVALID: useMetrics" + sys.exit() + +# Discretization +nd = 3 +if spaceOrder == 1: + hFactor=1.0 + if useHex: + basis=C0_AffineLinearOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,2) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,2) + else: + basis=C0_AffineLinearOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,3) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,3) +elif spaceOrder == 2: + hFactor=0.5 + if useHex: + basis=C0_AffineLagrangeOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,4) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,4) + else: + basis=C0_AffineQuadraticOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,4) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,4) + +# Domain and mesh + +#for debugging, make the tank short + +he = float(wavelength)/30.0#0.0#100 +L = (15.0, 1.5, 1.0) + +GenerationZoneLength = 1.2 +AbsorptionZoneLength= 2.8 +spongeLayer = True +xSponge = GenerationZoneLength +xRelaxCenter = xSponge/2.0 +epsFact_solid = xSponge/2.0 +#zone 2 +xSponge_2 = L[0] - AbsorptionZoneLength +xRelaxCenter_2 = 0.5*(xSponge_2+L[0]) +epsFact_solid_2 = AbsorptionZoneLength/2.0 + +weak_bc_penalty_constant = 100.0 +nLevels = 1 +#parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.element +parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.node +nLayersOfOverlapForParallel = 0 +structured=False + +if useHex: + nnx=4*Refinement+1 + nny=2*Refinement+1 + hex=True + domain = Domain.RectangularDomain(L) +else: + boundaries=['empty','left','right','bottom','top','front','back'] + boundaryTags=dict([(key,i+1) for (i,key) in enumerate(boundaries)]) + if structured: + nnx=4*Refinement + nny=2*Refinement + domain = Domain.RectangularDomain(L) + elif spongeLayer: + vertices=[[0.0,0.0,0.0],#0 + [xSponge,0.0,0.0],#1 + [xSponge_2,0.0,0.0],#2 + [L[0],0.0,0.0],#3 + [L[0],L[1],0.0],#4 + [xSponge_2,L[1],0.0],#5 + [xSponge,L[1],0.0],#6 + [0.0,L[1],0.0]]#7 + + + vertexFlags=[boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom']] + + + for v,vf in zip(vertices,vertexFlags): + vertices.append([v[0],v[1],L[2]]) + vertexFlags.append(boundaryTags['top']) + + print vertices + print vertexFlags + + segments=[[0,1], + [1,2], + [2,3], + [3,4], + [4,5], + [5,6], + [6,7], + [7,0], + [1,6], + [2,5]] + + segmentFlags=[boundaryTags['front'], + boundaryTags['front'], + boundaryTags['front'], + boundaryTags['right'], + boundaryTags['back'], + boundaryTags['back'], + boundaryTags['back'], + boundaryTags['left'], + boundaryTags['empty'], + boundaryTags['empty'] ] + + + facets=[] + facetFlags=[] + + for s,sF in zip(segments,segmentFlags): + facets.append([[s[0],s[1],s[1]+8,s[0]+8]]) + facetFlags.append(sF) + + bf=[[0,1,6,7],[1,2,5,6],[2,3,4,5]] + tf=[] + for i in range(0,3): + facets.append([bf[i]]) + tf=[ss + 8 for ss in bf[i]] + facets.append([tf]) + + for i in range(0,3): + facetFlags.append(boundaryTags['bottom']) + facetFlags.append(boundaryTags['top']) + + print facets + print facetFlags + + regions=[[xRelaxCenter, 0.5*L[1],0.0], + [xRelaxCenter_2, 0.5*L[1], 0.0], + [0.5*L[0],0.1*L[1], 0.0]] + regionFlags=[1,2,3] + + domain = Domain.PiecewiseLinearComplexDomain(vertices=vertices, + vertexFlags=vertexFlags, + facets=facets, + facetFlags=facetFlags, + regions=regions, + regionFlags=regionFlags, + ) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="KVApq1.4q12feena%21.16e" % ((he**3)/6.0,) + + + logEvent("""Mesh generated using: tetgen -%s %s""" % (triangleOptions,domain.polyfile+".poly")) + + porosityTypes = numpy.array([1.0, + 1.0, + 1.0, + 1.0]) + dragAlphaTypes = numpy.array([0.0, + 0.0, + 0.0, + 0.5/1.004e-6]) + + dragBetaTypes = numpy.array([0.0,0.0,0.0,0.0]) + + epsFact_solidTypes = np.array([0.0,0.0,0.0,0.0]) + + else: + vertices=[[0.0,0.0,0.0],#0 + [L[0],0.0,0.0],#1 + [L[0],L[1],0.0],#2 + [0.0,L[1],0.0]]#3 + + + vertexFlags=[boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['bottom']] + + + for v,vf in zip(vertices,vertexFlags): + vertices.append([v[0],v[1],L[2]]) + vertexFlags.append(boundaryTags['top']) + + segments=[[0,1], + [1,2], + [2,3], + [3,0]] + + segmentFlags=[boundaryTags['front'], + boundaryTags['right'], + boundaryTags['back'], + boundaryTags['left']] + + facets=[] + facetFlags=[] + + for s,sF in zip(segments,segmentFlags): + facets.append([[s[0],s[1],s[1]+4,s[0]+4]]) + facetFlags.append(sF) + + bf=[[0,1,2,3]] + tf=[] + for i in range(0,1): + facets.append([bf[i]]) + tf=[ss + 4 for ss in bf[i]] + facets.append([tf]) + + for i in range(0,1): + facetFlags.append(boundaryTags['bottom']) + facetFlags.append(boundaryTags['top']) + + for s,sF in zip(segments,segmentFlags): + segments.append([s[1]+4,s[0]+4]) + segmentFlags.append(sF) + + + regions=[[0.5*L[0],0.5*L[1], 0.0]] + regionFlags=[1] + + domain = Domain.PiecewiseLinearComplexDomain(vertices=vertices, + vertexFlags=vertexFlags, + facets=facets, + facetFlags=facetFlags, + regions=regions, + regionFlags=regionFlags) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + triangleOptions="KVApq1.4q12feena%21.16e" % ((he**3)/6.0,) + + + logEvent("""Mesh generated using: tetgen -%s %s""" % (triangleOptions,domain.polyfile+".poly")) + + + +# Time stepping +T=40*period +dt_fixed = period/11.0#2.0*0.5/20.0#T/2.0#period/21.0 +dt_init = min(0.001*dt_fixed,0.001) +runCFL=0.90 +nDTout = int(round(T/dt_fixed)) + +# Numerical parameters +ns_forceStrongDirichlet = False#True +backgroundDiffusionFactor=0.0 +if useMetrics: + ns_shockCapturingFactor = 0.25 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.25 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.25 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.25 + rd_lag_shockCapturing = False + epsFact_density = 3.0 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 1.5 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.1 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.1 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 +else: + ns_shockCapturingFactor = 0.9 + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = 0.9 + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = 1.0 + vof_shockCapturingFactor = 0.9 + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = 1.0 + rd_shockCapturingFactor = 0.9 + rd_lag_shockCapturing = False + epsFact_density = 1.5 + epsFact_viscosity = epsFact_curvature = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density + epsFact_redistance = 0.33 + epsFact_consrv_diffusion = 10.0 + redist_Newton = False + kappa_shockCapturingFactor = 0.9 + kappa_lag_shockCapturing = True#False + kappa_sc_uref = 1.0 + kappa_sc_beta = 1.0 + dissipation_shockCapturingFactor = 0.9 + dissipation_lag_shockCapturing = True#False + dissipation_sc_uref = 1.0 + dissipation_sc_beta = 1.0 + +ns_nl_atol_res = max(1.0e-10,0.0001*he**2) +vof_nl_atol_res = max(1.0e-10,0.0001*he**2) +ls_nl_atol_res = max(1.0e-10,0.0001*he**2) +rd_nl_atol_res = max(1.0e-10,0.005*he) +mcorr_nl_atol_res = max(1.0e-10,0.0001*he**2) +kappa_nl_atol_res = max(1.0e-10,0.001*he**2) +dissipation_nl_atol_res = max(1.0e-10,0.001*he**2) + +#turbulence +ns_closure=0 #1-classic smagorinsky, 2-dynamic smagorinsky, 3 -- k-epsilon, 4 -- k-omega +if useRANS == 1: + ns_closure = 3 +elif useRANS == 2: + ns_closure == 4 +# Water +rho_0 = 998.2 +nu_0 = 1.004e-6 + +# Air +rho_1 = 1.205 +nu_1 = 1.500e-5 + +# Surface tension +sigma_01 = 0.0 + +# Gravity +g = [0.0,0.0,-9.8] + +# Initial condition +waterLine_x = 2*L[0] +waterLine_z = inflowHeightMean + + +def signedDistance(x): + phi_x = x[0]-waterLine_x + phi_z = x[2]-waterLine_z + if phi_x < 0.0: + if phi_z < 0.0: + return max(phi_x,phi_z) + else: + return phi_z + else: + if phi_z < 0.0: + return phi_x + else: + return sqrt(phi_x**2 + phi_z**2) + + +def theta(x,t): + return k*x[0] - omega*t + pi/2.0 + +def z(x): + return x[2] - inflowHeightMean + +#sigma = omega - k*inflowVelocityMean[0] +h = inflowHeightMean # - transect[0][1] if lower left hand corner is not at z=0 + +#waveData + +def waveHeight(x,t): + return float(vegZoneInterp.interp_phi.__call__(t%vegZoneInterp.time[-1])) +def waveVelocity_u(x,t): + return vegZoneInterp.interpU.__call__(t%vegZoneInterp.time[-1],x[2])[0][0] +def waveVelocity_v(x,t): + return 0.0 +def waveVelocity_w(x,t): + return vegZoneInterp.interpW.__call__(t%vegZoneInterp.time[-1],x[2])[0][0] + +#solution variables + +def wavePhi(x,t): + return x[2] - waveHeight(x,t) + +def waveVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)) + +def twpflowVelocity_u(x,t): + waterspeed = waveVelocity_u(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + u = H*windVelocity[0] + (1.0-H)*waterspeed + return u + +def twpflowVelocity_v(x,t): + waterspeed = waveVelocity_v(x,t) + H = smoothedHeaviside(epsFact_consrv_heaviside*he,wavePhi(x,t)-epsFact_consrv_heaviside*he) + return H*windVelocity[1]+(1.0-H)*waterspeed + +def twpflowFlux(x,t): + return -twpflowVelocity_u(x,t) + +outflowHeight=inflowHeightMean + +def outflowVF(x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,x[2] - outflowHeight) + +def outflowPhi(x,t): + return x[2] - outflowHeight + +def outflowPressure(x,t): + if x[2]>inflowHeightMean: + return (L[2]-x[2])*rho_1*abs(g[2]) + else: + return (L[2]-inflowHeightMean)*rho_1*abs(g[2])+(inflowHeightMean-x[2])*rho_0*abs(g[2]) + + + #p_L = L[1]*rho_1*g[1] + #phi_L = L[1] - outflowHeight + #phi = x[1] - outflowHeight + #return p_L -g[1]*(rho_0*(phi_L - phi)+(rho_1 -rho_0)*(smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi_L) + # -smoothedHeaviside_integral(epsFact_consrv_heaviside*he,phi))) + +def twpflowVelocity_w(x,t): + return 0.0 + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + + +def zeroVel(x,t): + return 0.0 + +from collections import namedtuple + +RelaxationZone = namedtuple("RelaxationZone","center_x sign u v w") + +class RelaxationZoneWaveGenerator(AV_base): + """ Prescribe a velocity penalty scaling in a material zone via a Darcy-Forchheimer penalty + + :param zones: A dictionary mapping integer material types to Zones, where a Zone is a named tuple + specifying the x coordinate of the zone center and the velocity components + """ + def __init__(self,zones): + assert isinstance(zones,dict) + self.zones = zones + def calculate(self): + for l,m in enumerate(self.model.levelModelList): + for eN in range(m.coefficients.q_phi.shape[0]): + mType = m.mesh.elementMaterialTypes[eN] + if self.zones.has_key(mType): + for k in range(m.coefficients.q_phi.shape[1]): + t = m.timeIntegration.t + x = m.q['x'][eN,k] + m.coefficients.q_phi_solid[eN,k] = self.zones[mType].sign*(self.zones[mType].center_x - x[0]) + m.coefficients.q_velocity_solid[eN,k,0] = self.zones[mType].u(x,t) + m.coefficients.q_velocity_solid[eN,k,1] = self.zones[mType].v(x,t) + #m.coefficients.q_velocity_solid[eN,k,2] = self.zones[mType].w(x,t) + m.q['phi_solid'] = m.coefficients.q_phi_solid + m.q['velocity_solid'] = m.coefficients.q_velocity_solid + +rzWaveGenerator = RelaxationZoneWaveGenerator(zones={ + # 1:RelaxationZone(xRelaxCenter, + # 1.0, + # twpflowVelocity_u, + # twpflowVelocity_v, + # twpflowVelocity_w), + 1:RelaxationZone(xRelaxCenter_2, + -1.0, #currently Hs=1-exp_function + zeroVel, + zeroVel, + zeroVel)}) + +beam_quadOrder=3 +beam_useSparse=False +beamFilename="wavetankBeams" +#nBeamElements=max(nBeamElements,3) + +#beam info +beamLocation=[] +beamLength=[] +beamRadius=[] +EI=[] +GJ=[] +lam = 0.25 #0.05 #3.0*2.54/100.0 #57.4e-3 +lamx = 3.0**0.5*lam +xs = 1.2 +ys = 0.0 +xList=[] +yList = [] +while xs <= 11.0: + xList.append(xs) + xs += lam +while ys<= L[1]: + yList.append(ys) + ys+=lamx +for i in xList: + for j in yList: + beamLocation.append((i,j)) + beamLength.append(0.415) + beamRadius.append(0.0032) + EI.append(3.0e-4) # needs to be fixed + GJ.append(1.5e-4) # needs to be fixed + +xs = 1.2+0.5*lam +ys = 0.5*lamx +xList=[] +yList = [] +while xs <= 11.0: + xList.append(xs) + xs += lam + +while ys<= L[1]: + yList.append(ys) + ys+=lamx + +for i in xList: + for j in yList: + beamLocation.append((i,j)) + beamLength.append(0.415) + beamRadius.append(0.0032) + EI.append(3.0e-4) # needs to be fixed + GJ.append(1.5e-4) # needs to be fixed +nBeamElements = int(beamLength[0]/he*0.5) +nBeamElements=max(nBeamElements,3) +print nBeamElements diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/tank3D_so.py b/3d/waves_vegetation/wavetank_vegetation_time_series/tank3D_so.py new file mode 100644 index 00000000..f8c11c9b --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/tank3D_so.py @@ -0,0 +1,41 @@ +from proteus.default_so import * +import tank3D + +if tank3D.useOnlyVF: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n")] +else: + pnList = [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n"), + ("ls_p", "ls_n"), + ("redist_p", "redist_n"), + ("ls_consrv_p", "ls_consrv_n")] + + +if tank3D.useRANS > 0: + pnList.append(("kappa_p", + "kappa_n")) + pnList.append(("dissipation_p", + "dissipation_n")) +name = "tank3D_p" + +if tank3D.timeDiscretization == 'flcbdf': + systemStepControllerType = Sequential_MinFLCBDFModelStep + systemStepControllerType = Sequential_MinAdaptiveModelStep +else: + systemStepControllerType = Sequential_MinAdaptiveModelStep + +needEBQ_GLOBAL = False +needEBQ = False + +tnList = [0.0,tank3D.dt_init]+[i*tank3D.dt_fixed for i in range(1,tank3D.nDTout+1)] + +info = open("TimeList.txt","w") + + +for time in tnList: + info.write(str(time)+"\n") +info.close() + + +archiveFlag = ArchiveFlags.EVERY_SEQUENCE_STEP diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/twp_navier_stokes_n.py b/3d/waves_vegetation/wavetank_vegetation_time_series/twp_navier_stokes_n.py new file mode 100644 index 00000000..4e5ff8d1 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/twp_navier_stokes_n.py @@ -0,0 +1,69 @@ +from proteus import * +from twp_navier_stokes_p import * +from tank3D import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller_sys + stepController = Min_dt_cfl_controller + time_tol = 10.0*ns_nl_atol_res + atol_u = {1:time_tol,2:time_tol,3:time_tol} + rtol_u = {1:time_tol,2:time_tol,3:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis, + 1:basis, + 2:basis, + 3:basis} + +massLumping = False +numericalFluxType = None +conservativeFlux = None + +numericalFluxType = RANS2P.NumericalFlux +subgridError = RANS2P.SubgridError(coefficients,nd,lag=ns_lag_subgridError,hFactor=hFactor) +shockCapturing = RANS2P.ShockCapturing(coefficients,nd,ns_shockCapturingFactor,lag=ns_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None + +linearSmoother = SimpleNavierStokes3D + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'rans2p_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.01 +l_atol_res = 0.01*ns_nl_atol_res +nl_atol_res = ns_nl_atol_res +useEisenstatWalker = False +maxNonlinearIts = 50 +maxLineSearches = 0 +conservativeFlux = {0:'pwl-bdm-opt'} + + +#auxiliaryVariables=[pointGauges,rzWaveGenerator] diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/twp_navier_stokes_p.py b/3d/waves_vegetation/wavetank_vegetation_time_series/twp_navier_stokes_p.py new file mode 100644 index 00000000..9f5c7798 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/twp_navier_stokes_p.py @@ -0,0 +1,218 @@ +from proteus import * +from proteus.default_p import * +from tank3D import * +from proteus.mprans import RANS2P +from proteus.mprans import RANS2P_IB + +LevelModelType = RANS2P_IB.LevelModel +if useOnlyVF: + LS_model = None +else: + LS_model = 2 +if useRANS >= 1: + Closure_0_model = 5; Closure_1_model=6 + if useOnlyVF: + Closure_0_model=2; Closure_1_model=3 + if movingDomain: + Closure_0_model += 1; Closure_1_model += 1 +else: + Closure_0_model = None + Closure_1_model = None + +if spongeLayer or levee or slopingSpongeLayer: + coefficients = RANS2P_IB.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain, + porosityTypes=porosityTypes, + dragAlphaTypes=dragAlphaTypes, + dragBetaTypes=dragBetaTypes, + epsFact_solid = epsFact_solidTypes, + beamLocation=beamLocation, + beamLength=beamLength, + beamRadius=beamRadius, + EI=EI, + GJ=GJ, + nBeamElements=nBeamElements, + beam_quadOrder=beam_quadOrder, + beamFilename=beamFilename, + beam_Cd=1.2, + beam_nlTol=1.0e-5, + beamRigid=True) +else: + coefficients = RANS2P_IB.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + VF_model=1, + LS_model=LS_model, + Closure_0_model=Closure_0_model, + Closure_1_model=Closure_1_model, + epsFact_density=epsFact_density, + stokes=False, + useVF=useVF, + useRBLES=useRBLES, + useMetrics=useMetrics, + eb_adjoint_sigma=1.0, + eb_penalty_constant=weak_bc_penalty_constant, + forceStrongDirichlet=ns_forceStrongDirichlet, + turbulenceClosureModel=ns_closure, + movingDomain=movingDomain, + beamLocation=beamLocation, + beamLength=beamLength, + beamRadius=beamRadius, + EI=EI, + GJ=GJ, + nBeamElements=nBeamElements, + beam_quadOrder=beam_quadOrder, + beamFilename=beamFilename, + beam_Cd=1.2, + beam_nlTol=1.0e-5, + beamRigid=True) + + +def getDBC_p(x,flag): + if flag == boundaryTags['top']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return outflowPressure + +def getDBC_u(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_u + + +def getDBC_v(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_v + +def getDBC_w(x,flag): + if flag == boundaryTags['left']: + return twpflowVelocity_w + +dirichletConditions = {0:getDBC_p, + 1:getDBC_u, + 2:getDBC_v, + 3:getDBC_w} + +def getAFBC_p(x,flag): + if flag == boundaryTags['left']: + return lambda x,t: -twpflowVelocity_u(x,t) + elif flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getAFBC_u(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getAFBC_v(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getAFBC_w(x,flag): + if flag == boundaryTags['bottom'] or flag == boundaryTags['back'] or flag == boundaryTags['front']: + return lambda x,t: 0.0 + elif flag == boundaryTags['right']: + return lambda x,t: 0.0 + else: + return None + +def getDFBC_u(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + + +def getDFBC_v(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 + + +def getDFBC_w(x,flag): + if flag != boundaryTags['left']: + return lambda x,t: 0.0 +# if flag == boundaryTags['bottom']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['back']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['front']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['right']: +# return lambda x,t: 0.0 +# elif flag == boundaryTags['top']: +# return lambda x,t: 0.0 + + +advectiveFluxBoundaryConditions = {0:getAFBC_p, + 1:getAFBC_u, + 2:getAFBC_v, + 3:getAFBC_w} + +diffusiveFluxBoundaryConditions = {0:{}, + 1:{1:getDFBC_u}, + 2:{2:getDFBC_v}, + 3:{3:getDFBC_w}} + +class PerturbedSurface_p: + def __init__(self,waterLevel): + self.waterLevel=waterLevel + def uOfXT(self,x,t): + # if signedDistance(x) < 0: + # return -(L[2] - self.waterLevel)*rho_1*g[2] - (self.waterLevel - x[2])*rho_0*g[2] + # else: + # return -(L[2] - self.waterLevel)*rho_1*g[2] + if x[2]>inflowHeightMean: + return (L[2]-x[2])*rho_1*abs(g[2]) + else: + return (L[2]-inflowHeightMean)*rho_1*abs(g[2])+(inflowHeightMean-x[2])*rho_0*abs(g[2]) + +class AtRest: + def __init__(self): + pass + def uOfXT(self,x,t): + return 0.0 + +class WaterVelocity_u: + def __init__(self): + pass + def uOfXT(self,x,t): + + return waterVelocity(x,0) + +initialConditions = {0:PerturbedSurface_p(waterLine_z), + 1:AtRest(), + 2:AtRest(), + 3:AtRest()} diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/vof_n.py b/3d/waves_vegetation/wavetank_vegetation_time_series/vof_n.py new file mode 100644 index 00000000..de2fdd09 --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/vof_n.py @@ -0,0 +1,65 @@ +from proteus import * +from tank3D import * +from vof_p import * + +if timeDiscretization=='vbdf': + timeIntegration = VBDF + timeOrder=2 + stepController = Min_dt_cfl_controller +elif timeDiscretization=='flcbdf': + timeIntegration = FLCBDF + #stepController = FLCBDF_controller + stepController = Min_dt_cfl_controller + time_tol = 10.0*vof_nl_atol_res + atol_u = {0:time_tol} + rtol_u = {0:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:basis} + +massLumping = False +numericalFluxType = VOF.NumericalFlux +conservativeFlux = None +subgridError = VOF.SubgridError(coefficients=coefficients,nd=nd) +shockCapturing = VOF.ShockCapturing(coefficients,nd,shockCapturingFactor=vof_shockCapturingFactor,lag=vof_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = SparseMatrix + +if useOldPETSc: + multilevelLinearSolver = PETSc + levelLinearSolver = PETSc +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + +if useSuperlu: + multilevelLinearSolver = LU + levelLinearSolver = LU + +linear_solver_options_prefix = 'vof_' +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +nl_atol_res = vof_nl_atol_res + +linTolFac = 0.0 +l_atol_res = 0.1*vof_nl_atol_res + +useEisenstatWalker = False + +maxNonlinearIts = 50 +maxLineSearches = 0 + +auxiliaryVariables = [columnGauge] + diff --git a/3d/waves_vegetation/wavetank_vegetation_time_series/vof_p.py b/3d/waves_vegetation/wavetank_vegetation_time_series/vof_p.py new file mode 100644 index 00000000..ea33ff3b --- /dev/null +++ b/3d/waves_vegetation/wavetank_vegetation_time_series/vof_p.py @@ -0,0 +1,47 @@ +from proteus import * +from proteus.default_p import * +from proteus.ctransportCoefficients import smoothedHeaviside +from tank3D import * +from proteus.mprans import VOF + +LevelModelType = VOF.LevelModel +if useOnlyVF: + RD_model = None + LS_model = None +else: + RD_model = 3 + LS_model = 2 + +coefficients = VOF.Coefficients(LS_model=LS_model,V_model=0,RD_model=RD_model,ME_model=1, + checkMass=False,useMetrics=useMetrics, + epsFact=epsFact_vof,sc_uref=vof_sc_uref,sc_beta=vof_sc_beta,movingDomain=movingDomain) + +def getDBC_vof(x,flag): + if flag == boundaryTags['left']: + return waveVF + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return lambda x,t: 1.0 +# elif flag == boundaryTags['right']: +# return outflowVF + + +dirichletConditions = {0:getDBC_vof} + +def getAFBC_vof(x,flag): + if flag == boundaryTags['left']: + return None + elif flag == boundaryTags['top']:# or x[1] >= L[1] - 1.0e-12: + return None +# elif flag == boundaryTags['right']: +# return None + else: + return lambda x,t: 0.0 + +advectiveFluxBoundaryConditions = {0:getAFBC_vof} +diffusiveFluxBoundaryConditions = {0:{}} + +class PerturbedSurface_H: + def uOfXT(self,x,t): + return smoothedHeaviside(epsFact_consrv_heaviside*he,signedDistance(x)) + +initialConditions = {0:PerturbedSurface_H()} diff --git a/OpenFOAM_Benchmarks/README.rst b/OpenFOAM_Benchmarks/README.rst index 39b88474..a9ca2965 100644 --- a/OpenFOAM_Benchmarks/README.rst +++ b/OpenFOAM_Benchmarks/README.rst @@ -2,18 +2,18 @@ OpenFOAM Benchmarks ===================================================== -https://github.com/erdc-cm/air-water-vv/OpenFOAM_Benchmarks/ +https://github.com/erdc/air-water-vv/OpenFOAM_Benchmarks/ Case descritpion ---------------------------- -- dambreakColagrossi, see https://github.com/erdc-cm/air-water-vv/2d/dambreak_Colagrossi/README.rst +- dambreakColagrossi, see https://github.com/erdc/air-water-vv/2d/dambreak_Colagrossi/README.rst -- dambreakWithObstacle, see http://foam.sourceforge.net/docs/Guides-a4/UserGuide.pdf (section 2.3) and https://github.com/erdc-cm/air-water-vv/2d/dambreak_Ubbink/README.rst +- dambreakWithObstacle, see http://foam.sourceforge.net/docs/Guides-a4/UserGuide.pdf (section 2.3) and https://github.com/erdc/air-water-vv/2d/dambreak_Ubbink/README.rst - lidDrivenCavity, see http://foam.sourceforge.net/docs/Guides-a4/UserGuide.pdf (Section 2.1) -- dambreakGomez, see https://github.com/erdc-cm/air-water-vv/blob/master/3d/dambreak_Gomez/README.rst +- dambreakGomez, see https://github.com/erdc/air-water-vv/blob/master/3d/dambreak_Gomez/README.rst How to set up and run OpenFOAM cases ------------------------ diff --git a/README.rst b/README.rst index 33fef22c..dc1b0021 100644 --- a/README.rst +++ b/README.rst @@ -2,7 +2,7 @@ Verification and validation for air/water flow models ===================================================== -https://github.com/erdc-cm/air-water-vv +https://github.com/erdc/air-water-vv Organization of the test set ---------------------------- @@ -16,7 +16,7 @@ test consists of a single directory including - All data files describing the geometry and physical parameters of the problem - A set of input files for a code - (e.g. https://github.com/erdc-cm/proteus) + (e.g. https://github.com/erdc/proteus) Adding new test problems ------------------------ diff --git a/doc/source/chl.png b/doc/source/chl.png index 3f4fad58..5d51c44e 100644 Binary files a/doc/source/chl.png and b/doc/source/chl.png differ diff --git a/doc/source/hrw.png b/doc/source/hrw.png index c7a23dae..02332a52 100644 Binary files a/doc/source/hrw.png and b/doc/source/hrw.png differ diff --git a/inputTemplates/petsc.options.superlu_dist b/inputTemplates/petsc.options.superlu_dist index a14e35cf..1d7081cd 100644 --- a/inputTemplates/petsc.options.superlu_dist +++ b/inputTemplates/petsc.options.superlu_dist @@ -5,4 +5,5 @@ -mcorr_ksp_type preonly -mcorr_pc_type lu -mcorr_pc_factor_mat_solver_package superlu_dist -kappa_ksp_type preonly -kappa_pc_type lu -kappa_pc_factor_mat_solver_package superlu_dist -dissipation_ksp_type preonly -dissipation_pc_type lu -dissipation_pc_factor_mat_solver_package superlu_dist +-mesh_ksp_type preonly -mesh_pc_type lu -mesh_pc_factor_mat_solver_package superlu_dist diff --git a/private/proteus-mprans.yaml b/private/proteus-mprans.yaml index 8e1a4d58..f7cea748 100644 --- a/private/proteus-mprans.yaml +++ b/private/proteus-mprans.yaml @@ -4,7 +4,7 @@ dependencies: run: [proteus, numpy] sources: - - url: https://github.com/erdc-cm/proteus-mprans + - url: https://github.com/erdc/proteus-mprans key: git:fc5a90c2a50bb9716460876ad1ee91da598ff8f3 profile_links: