diff --git a/.travis.yml b/.travis.yml index 7a17e53b..6fd60187 100644 --- a/.travis.yml +++ b/.travis.yml @@ -29,7 +29,7 @@ before_install: install: - git lfs fetch - git lfs checkout -- git clone https://github.com/erdc-cm/proteus +- git clone https://github.com/erdc/proteus - cd proteus - make hashdist - make stack diff --git a/2d/floatingStructures/floatingCaissonAddedMass/Heave.ipynb b/2d/floatingStructures/floatingCaissonAddedMass/Heave.ipynb index 2b886eed..d210c1ad 100644 --- a/2d/floatingStructures/floatingCaissonAddedMass/Heave.ipynb +++ b/2d/floatingStructures/floatingCaissonAddedMass/Heave.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 72, + "execution_count": 43, "metadata": {}, "outputs": [], "source": [ @@ -13,17 +13,26 @@ }, { "cell_type": "code", - "execution_count": 73, + "execution_count": 44, "metadata": {}, "outputs": [], "source": [ - "caisson=np.loadtxt(\"record_bridge.csv\",delimiter=\",\",skiprows=1)\n" + "caisson_am=np.loadtxt(\"record_caisson_am.csv\",delimiter=\",\",skiprows=1)\n", + "caisson=np.loadtxt(\"record_caisson.csv\",delimiter=\",\",skiprows=1)\n", + "caisson_am_nc=np.loadtxt(\"record_caisson_am_nc.csv\",delimiter=\",\",skiprows=1)\n", + "caisson_nc=np.loadtxt(\"record_caisson_nc.csv\",delimiter=\",\",skiprows=1)\n", + "caisson_chrono=np.loadtxt(\"record_caisson_chrono.csv\",delimiter=\",\",skiprows=1)\n", + "caisson_chrono_am=np.loadtxt(\"record_caisson_chrono_am.csv\",delimiter=\",\",skiprows=1)\n", + "caisson_chrono_am_nc=np.loadtxt(\"record_caisson_chrono_am_nc.csv\",delimiter=\",\",skiprows=1)\n", + "caisson_chrono_nc=np.loadtxt(\"record_caisson_chrono_nc.csv\",delimiter=\",\",skiprows=1)" ] }, { "cell_type": "code", - "execution_count": 74, - "metadata": {}, + "execution_count": 45, + "metadata": { + "scrolled": true + }, "outputs": [ { "data": { @@ -805,7 +814,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -817,33 +826,1694 @@ { "data": { "text/plain": [ - "[]" + "" ] }, - "execution_count": 74, + "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "plt.plot(caisson[:,0],caisson[:,2])" + "import math\n", + "fig = plt.figure()\n", + "#plt.plot(caisson[:,0],caisson[:,2],)\n", + "#plt.plot(caisson_am[:,0],caisson_am[:,2])\n", + "plt.plot(caisson_nc[:,0],caisson_nc[:,2])\n", + "#caisson_nc[:,0]= np.linspace(0.0,math.pi,caisson_nc.shape[0])\n", + "plt.plot(caisson_nc[:,0],-0.125*np.cos(caisson_nc[:,0]*math.pi)+caisson_nc[0,2]+0.125)\n", + "#plt.plot(caisson_am_nc[:,0],caisson_am_nc[:,2])\n", + "#plt.plot(caisson_chrono[:,0],caisson_chrono[:,4])\n", + "#plt.plot(caisson_chrono_am[:,0],caisson_chrono_am[:,4])\n", + "#plt.plot(caisson_chrono_am_nc[:,0],caisson_chrono_am_nc[:,4])\n", + "#plt.plot(caisson_chrono_am_nc[:,0],np.sin(caisson_chrono_am_nc[:,0]*4-3.14/2.0))\n", + "#plt.plot(caisson_chrono_nc[:,0],caisson_chrono_nc[:,4])\n", + "#plt.legend(['c','c+am','nc','nc+am','CH+c','CH+c+am','CH+nc','CH+nc+am'])\n", + "plt.legend(['custom','custom+am','chrono+am','chrono'])" ] }, { "cell_type": "code", - "execution_count": 75, + "execution_count": 46, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "0.0020000000000000018" + ] + }, + "execution_count": 46, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "plt.savefig(\"heave.png\")" + "plt.savefig(\"heave.png\")\n", + "1.19575-1.19375" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 47, "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " // select the cell after this one\n", + " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", + " IPython.notebook.select(index + 1);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fig = plt.figure()\n", + "#plt.plot(caisson[:,0],caisson[:,2],)\n", + "#plt.plot(caisson_am[:,0],caisson_am[:,2])\n", + "plt.plot(caisson_nc[:,0],caisson_nc[:,6])\n", + "plt.plot(caisson_nc[:,0],-(math.pi/4.0)*np.cos(caisson_nc[:,0]*math.pi))\n", + "#plt.plot(caisson_am_nc[:,0],caisson_am_nc[:,6])\n", + "#plt.legend(['c','c+am','nc','nc+am'])\n", + "#plt.plot(caisson_chrono[:,0],caisson_chrono[:,4])\n", + "#plt.plot(caisson_chrono_am[:,0],caisson_chrono_am[:,4])\n", + "#plt.plot(caisson_chrono_am_nc[:,0],caisson_chrono_am_nc[:,9])\n", + "#plt.legend(['am','Ch+am'])\n", + "#plt.plot(caisson_chrono_nc[:,0],caisson_chrono_nc[:,4])\n", + "#plt.legend(['c','c+am','nc','nc+am','CH+c','CH+c+am','CH+nc','CH+nc+am'])\n" + ] } ], "metadata": { diff --git a/2d/floatingStructures/floatingCaissonAddedMass/MeshRefinement.py b/2d/floatingStructures/floatingCaissonAddedMass/MeshRefinement.py new file mode 100644 index 00000000..5c0d3e05 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/MeshRefinement.py @@ -0,0 +1,75 @@ +def geometry_to_gmsh(domain): + import py2gmsh + from py2gmsh.Mesh import * + from py2gmsh.Entity import * + from py2gmsh.Fields import * + self = domain + lines_dict = {} + + mesh = Mesh() + + if self.boundaryTags: + for tag, flag in self.boundaryTags.items(): + phys = PhysicalGroup(nb=flag, name=tag) + mesh.addGroup(phys) + + for i, v in enumerate(self.vertices): + if domain.nd == 2: + p = Point([v[0], v[1], 0.]) + else: + p = Point(v) + mesh.addEntity(p) + g = mesh.groups.get(self.vertexFlags[i]) + if g: + g.addEntity(p) + nb_points = i+1 + for i in range(nb_points): + lines_dict[i] = {} + + for i, s in enumerate(self.segments): + lines_dict[s[0]][s[1]] = i + l = Line([mesh.points[s[0]+1], mesh.points[s[1]+1]]) + mesh.addEntity(l) + g = mesh.groups.get(self.segmentFlags[i]) + if g: + g.addEntity(l) + + for i, f in enumerate(self.facets): + if self.nd == 3 or (self.nd == 2 and i not in self.holes_ind): + lineloops = [] + for j, subf in enumerate(f): + lineloop = [] + # vertices in facet + for k, ver in enumerate(subf): + if ver in lines_dict[subf[k-1]].keys(): + lineloop += [lines_dict[subf[k-1]][ver]+1] + elif subf[k-1] in lines_dict[ver].keys(): + # reversed + lineloop += [(lines_dict[ver][subf[k-1]]+1)] + else: + l = Line([mesh.points[subf[k-1]+1], mesh.points[ver+1]]) + mesh.addEntity(l) + lineloop += [l.nb] + ll = LineLoop(mesh.getLinesFromIndex(lineloop)) + mesh.addEntity(ll) + lineloops += [ll.nb] + s = PlaneSurface([mesh.lineloops[loop] for loop in lineloops]) + mesh.addEntity(s) + g = mesh.groups.get(self.facetFlags[i]) + if g: + g.addEntity(s) + + for i, V in enumerate(self.volumes): + surface_loops = [] + hole_loops = [] + for j, sV in enumerate(V): + sl = SurfaceLoop(mesh.getSurfacesFromIndex((np.array(sV)+1).tolist())) + mesh.addEntity(sl) + surface_loops += [sl] + vol = Volume(surface_loops) + mesh.addEntity(vol) + g = mesh.groups.get(self.regionFlags[i]) + if g: + g.addEntity(vol) + + return mesh diff --git a/2d/floatingStructures/floatingCaissonAddedMass/added_massIB_n.py b/2d/floatingStructures/floatingCaissonAddedMass/added_massIB_n.py new file mode 100644 index 00000000..ffd0409e --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/added_massIB_n.py @@ -0,0 +1,56 @@ +from proteus.default_n import * +import added_massIB_p as physics +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools) +from proteus.mprans import AddedMass +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +triangleOptions = mesh.triangleOptions + +femSpaces = {0:ct.basis} + +stepController=StepControl.FixedStep + +numericalFluxType = AddedMass.NumericalFlux + +matrix = LinearAlgebraTools.SparseMatrix + + +multilevelLinearSolver = LinearSolvers.KSP_petsc4py +levelLinearSolver = LinearSolvers.KSP_petsc4py +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +nonlinearSmoother = NonlinearSolvers.AddedMassNewton +linearSmoother = LinearSolvers.NavierStokesPressureCorrection + +linear_solver_options_prefix = 'am_' + +multilevelNonlinearSolver = NonlinearSolvers.AddedMassNewton +levelNonlinearSolver = NonlinearSolvers.AddedMassNewton + +#linear solve rtolerance + +linTolFac = 0.0 +l_atol_res = 1.0e-10#0.001*ct.am_nl_atol_res +tolFac = 0.0 +nl_atol_res = ct.am_nl_atol_res + +nonlinearSolverConvergenceTest = 'rits' +levelNonlinearSolverConvergenceTest = 'rits' +linearSolverConvergenceTest = 'r-true' +maxNonlinearIts =1 +maxLineSearches=0 +periodicDirichletConditions=None +conservativeFlux=None +auxiliaryVariables=[ct.system.ProtChAddedMass] diff --git a/2d/floatingStructures/floatingCaissonAddedMass/added_massIB_p.py b/2d/floatingStructures/floatingCaissonAddedMass/added_massIB_p.py new file mode 100644 index 00000000..784dfda7 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/added_massIB_p.py @@ -0,0 +1,44 @@ +from proteus.default_p import * +from proteus.mprans import AddedMass +from proteus import Context +name = "added_mass" +ct = Context.get() +domain = ct.domain +nd = domain.nd +mesh = domain.MeshOptions +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T + +LevelModelType = AddedMass.LevelModel + +coefficients = AddedMass.Coefficients(nd=nd, + V_model=4, + barycenters=domain.barycenters, + flags_rigidbody=ct.flags_rigidbody) + +def getDBC_am(x,flag): +# if flag == ct.domain.boundaryTags['tank2D1_y+']: +# return lambda x,t: 0.0 +# else: + return None + +dirichletConditions = {0:getDBC_am} + +advectiveFluxBoundaryConditions = {} + +def getFlux_am(x, flag): + #the unit rigid motions will be applied internally + #leave this set to zero +# if flag == ct.domain.boundaryTags['tank2D1_y+']: +# return None +# else: + return lambda x,t: 0.0 + +diffusiveFluxBoundaryConditions = {0: {0:getFlux_am}} + +class dp_IC: + def uOfXT(self, x, t): + return 0.0 + +initialConditions = {0: dp_IC()} diff --git a/2d/floatingStructures/floatingCaissonAddedMass/added_mass_n.py b/2d/floatingStructures/floatingCaissonAddedMass/added_mass_n.py index 9fab73b8..a27ae72b 100644 --- a/2d/floatingStructures/floatingCaissonAddedMass/added_mass_n.py +++ b/2d/floatingStructures/floatingCaissonAddedMass/added_mass_n.py @@ -42,7 +42,7 @@ #linear solve rtolerance linTolFac = 0.0 -l_atol_res = ct.am_nl_atol_res +l_atol_res = 1.0e-10#0.001*ct.am_nl_atol_res tolFac = 0.0 nl_atol_res = ct.am_nl_atol_res diff --git a/2d/floatingStructures/floatingCaissonAddedMass/added_mass_p.py b/2d/floatingStructures/floatingCaissonAddedMass/added_mass_p.py index 5ec10895..f7be3cfc 100644 --- a/2d/floatingStructures/floatingCaissonAddedMass/added_mass_p.py +++ b/2d/floatingStructures/floatingCaissonAddedMass/added_mass_p.py @@ -18,6 +18,9 @@ flags_rigidbody=ct.flags_rigidbody) def getDBC_am(x,flag): +# if flag == ct.domain.boundaryTags['tank2D1_y+']: +# return lambda x,t: 0.0 +# else: return None dirichletConditions = {0:getDBC_am} @@ -25,8 +28,11 @@ def getDBC_am(x,flag): advectiveFluxBoundaryConditions = {} def getFlux_am(x, flag): - #the unit rigid motions will applied internally + #the unit rigid motions will be applied internally #leave this set to zero +# if flag == ct.domain.boundaryTags['tank2D1_y+']: +# return None +# else: return lambda x,t: 0.0 diffusiveFluxBoundaryConditions = {0: {0:getFlux_am}} diff --git a/2d/floatingStructures/floatingCaissonAddedMass/caisson2D_oscillation.py b/2d/floatingStructures/floatingCaissonAddedMass/caisson2D_oscillation.py index 614034cf..5ef3aced 100755 --- a/2d/floatingStructures/floatingCaissonAddedMass/caisson2D_oscillation.py +++ b/2d/floatingStructures/floatingCaissonAddedMass/caisson2D_oscillation.py @@ -3,13 +3,14 @@ from proteus import WaveTools as wt from proteus.Profiling import logEvent #import ChRigidBody as crb -from proteus.mbd import ChRigidBody as crb +from proteus.mbd import CouplingFSI as crb from math import * import numpy as np - opts=Context.Options([ - ("use_chrono", True, "use chrono (True) or custom solver"), + ("cf",1.0,"Use corrected form of added mass"), + ("nc_form",True,"Use non-conservative NSE form"), + ("use_chrono", False, "use chrono (True) or custom solver"), # predefined test cases ("water_level", 1.5, "Height of free surface above bottom"), # tank @@ -19,14 +20,14 @@ ("addedMass", True, "added mass"), ("caisson", True, "caisson"), ("caisson_dim", (0.5, 0.5), "Dimensions of the caisson"), - ("caisson_coords", (1.5, 1.5-0.05), "Dimensions of the caisson"), + ("caisson_coords", (1.5, 1.5+0.5), "Dimensions of the caisson"), ("free_x", (0.0, 1.0, 0.0), "Translational DOFs"), ("free_r", (0.0, 0.0, 0.0), "Rotational DOFs"), ("VCG", 0.05, "vertical position of the barycenter of the caisson"), ("caisson_mass", 125., "Mass of the caisson"), ("caisson_inertia", 4.05, "Inertia of the caisson"), ("rotation_angle", 0., "Initial rotation angle (in degrees)"), - ("chrono_dt", 0.00001, "time step of chrono"), + ("chrono_dt", 0.001, "time step of chrono"), # mesh refinement ("refinement", True, "Gradual refinement"), ("he", 0.03, "Set characteristic element size"), @@ -44,8 +45,7 @@ ("useRANS", 0, "RANS model"), ("sc", 0.25, "shockCapturing factor"), ("weak_factor", 10., "weak bc penalty factor"), - ("strong_dir", False, "strong dirichlet (True/False)"), - ]) + ("strong_dir", False, "strong dirichlet (True/False)")]) @@ -71,7 +71,13 @@ inertia = opts.caisson_inertia caisson_dim = opts.caisson_dim + caisson_name='caisson_chrono' + if opts.addedMass: + caisson_name+='_am' + if opts.nc_form: + caisson_name+='_nc' caisson = st.Rectangle(domain, dim=opts.caisson_dim, coords=[0.,0.], barycenter=np.zeros(3)) + caisson.name=caisson_name ang = rotation_angle caisson.setHoles([[0., 0.]]) caisson.holes_ind = np.array([0]) @@ -131,15 +137,26 @@ def calculate(self): body.It = np.array([[1., 0., 0.], [0., 1., 0.], [0., 0., inertia]]) - body.setRecordValues(filename='record_bridge', all_values=True) + filename='record_caisson' + if opts.addedMass: + filename+='_am' + if opts.nc_form: + filename+='_nc' + if opts.cf: + filename+='_cf' + body.setRecordValues(filename=filename, all_values=True) body.coords_system = caisson.coords_system # hack body.last_Aij = np.zeros((6,6),'d') + body.last_a = np.zeros((6,),'d') + body.a = np.zeros((6,),'d') body.last_Omega = np.zeros((3,3),'d') body.last_velocity = np.zeros((3,),'d') body.last_Q = np.eye(3) body.last_h = np.zeros((3,),'d') body.last_mom = np.zeros((6,),'d') body.last_u = np.zeros((18,),'d') + body.last_t = 0.0 + body.t = 0.0 body.h = np.zeros((3,),'d') body.velocity = np.zeros((3,),'d') body.mom = np.zeros((6,),'d') @@ -164,8 +181,11 @@ def calculate(self): # body.ProtChAddedMass = crb.ProtChAddedMass(system) def step(dt): + from math import ceil logEvent("Barycenter "+str(body.barycenter)) - def F(u): + n = max(1.0,ceil(dt/opts.chrono_dt)) + DT=dt/n + def F(u,theta): """The residual and Jacobian for discrete 6DOF motion with added mass""" v = u[:3] omega = u[3:6] @@ -175,90 +195,117 @@ def F(u): [ omega[2], 0.0, -omega[0]], [-omega[1], omega[0], 0.0]]) I = np.matmul(np.matmul(Q, body.It), Q.transpose()) - body.Aij = body.bodyAddedMass.model.levelModelList[-1].Aij[1].copy() + body.Aij = np.zeros((6,6),'d') + if opts.addedMass: + for i in range(1,5):#number of rigid body facets + body.Aij += body.bodyAddedMass.model.levelModelList[-1].Aij[i] avg_Aij=False if avg_Aij: - M = 0.5*(body.Aij + body.last_Aij) + M = body.Aij*theta + body.last_Aij*(1-theta) else: M = body.Aij.copy() for i in range(6): for j in range(6): M[i,j]*=body.free_dof[j]#only allow j accelerations to contribute to i force balance if j is free M[j,i]*=body.free_dof[j]#only allow j added mass forces if j is free + body.Aij[i,j]*=body.free_dof[j]#only allow j accelerations to contribute to i force balance if j is free + body.Aij[j,i]*=body.free_dof[j]#only allow j added mass forces if j is free body.FT[:3] = body.F body.FT[3:] = body.M + #cek debug + #body.FT[:] = 0.0 + #body.FT[1] = body.mass * 0.125*math.pi**2 * math.cos(body.t*math.pi) for i in range(3): M[i, i] += body.mass for j in range(3): M[3+i, 3+j] += I[i, j] r = np.zeros((18,),'d') - body_cons=True - if body_cons: - r[:6] = np.matmul(M, u[:6] - body.last_u[:6]) - dt*0.5*(body.FT+body.last_FT) + BE=True + CF=opts.cf + if BE: + r[:6] = np.matmul(M, u[:6]) - np.matmul(body.Aij, body.last_u[:6]) - body.last_mom - DT*body.FT - CF*np.matmul(body.Aij,body.last_a) + r[6:9] = h - body.last_h - DT*v + rQ = Q - body.last_Q - DT*np.matmul(Omega,Q) else: - r[:6] = np.matmul(M, u[:6]) - body.last_mom - dt*0.5*(body.FT+body.last_FT) - r[6:9] = h - body.last_h - dt*0.5*(v + body.last_velocity) - rQ = Q - body.last_Q - dt*0.25*np.matmul((Omega + body.last_Omega),(Q+body.last_Q)) + r[:6] = np.matmul(M, u[:6]) - np.matmul(body.Aij, body.last_u[:6]) - body.last_mom - DT*(body.FT*theta+body.last_FT*(1.0-theta)) - CF*np.matmul(body.Aij,body.last_a) + r[6:9] = h - body.last_h - DT*0.5*(v + body.last_velocity) + rQ = Q - body.last_Q - DT*0.25*np.matmul((Omega + body.last_Omega),(Q+body.last_Q)) r[9:18] = rQ.flatten() J = np.zeros((18,18),'d') J[:6,:6] = M #neglecting 0:6 dependence on Q for i in range(3): - J[6+i,i] = -dt*0.5 + if BE: + J[6+i,i] = -DT + else: + J[6+i,i] = -DT*0.5 J[6+i,6+i] = 1.0 for i in range(9): J[9+i,9+i] = 1.0 for i in range(3): for j in range(3): - J[9+i*3+j, 9+i+j*3] -= dt*0.25*(Omega+body.last_Omega)[i,j] - #cek, this indexing may be correct too, not sure - #for j in range(9,18): - # row = (j-9)/3 - # col = j%3 - # print("jjj", row, col) - # J[i,j] = 1.0 - dt*0.25*(Omega[row, col]+body.last_Omega[row, col]) - + if BE: + J[9+i*3+j, 9+i+j*3] -= DT*Omega[i,j] + else: + J[9+i*3+j, 9+i+j*3] -= DT*0.25*(Omega+body.last_Omega)[i,j] #neglecting 9:18 dependence on omega body.Omega[:] = Omega body.velocity[:] = v body.Q[:] = Q body.h[:] = h body.u[:] = u - body.mom[:] = np.matmul(M,u[:6]) + body.mom[:3] = body.mass*u[:3] + body.mom[3:6] = np.matmul(I,u[3:6]) + body.a[:] = u[:6] - body.last_u[:6] return r, J - # if body.ProtChAddedMass.model is not None: - # body.Aij[:] = body.ProtChAddedMass.model.Aij - # print("body Aij", body.Aij) nd = body.nd - u = np.zeros((18,),'d') - u[:] = body.last_u - r = np.zeros((18,),'d') - r,J = F(u) - its=0 - maxits=100 - while ((its==0 or np.linalg.norm(r) > 1.0e-8) and its < maxits): - u -= np.linalg.solve(J,r) - r,J = F(u) - its+=1 - logEvent("6DOF its = " + `its` + " residual = "+`r`) - logEvent("displacement, h = "+`body.h`) - logEvent("rotation, Q = "+`body.Q`) - logEvent("velocity, v = "+`body.velocity`) - logEvent("angular acceleration matrix, Omega = "+`body.Omega`) - body.last_Aij[:]=body.Aij - body.last_FT[:] = body.FT - body.last_Omega[:] = body.Omega - body.last_velocity[:] = body.velocity - body.last_Q[:] = body.Q - body.last_h[:] = body.h - body.last_u[:] = body.u - body.last_mom[:] = body.mom + Q_start=body.Q.copy() + h_start=body.last_h.copy() + for i in range(int(n)): + theta = (i+1)*DT/dt + body.t = body.last_t + DT + logEvent("6DOF theta {0}".format(theta)) + u = np.zeros((18,),'d') + u[:] = body.last_u + r = np.zeros((18,),'d') + r,J = F(u,theta) + its=0 + maxits=100 + while ((its==0 or np.linalg.norm(r) > 1.0e-8) and its < maxits): + u -= np.linalg.solve(J,r) + r,J = F(u,theta) + its+=1 + logEvent("6DOF its {0}".format(its)) + logEvent("6DOF res {0}".format(np.linalg.norm(r))) + body.last_Aij[:]=body.Aij + body.last_FT[:] = body.FT + body.last_Omega[:] = body.Omega + body.last_velocity[:] = body.velocity + body.last_Q[:] = body.Q + body.last_h[:] = body.h + body.last_u[:] = body.u + body.last_mom[:] = body.mom + body.last_t = body.t + body.last_a[:] = body.a # translate and rotate + body.h -= h_start + body.last_h[:] = 0.0 body.last_position[:] = body.position - body.rotation_matrix[:] = np.linalg.solve(body.last_Q,body.Q) + #body.rotation_matrix[:] = np.linalg.solve(Q_start,body.Q) + body.rotation_matrix[:] = np.matmul(np.linalg.inv(body.Q),Q_start) + body.rotation_euler[2] -= math.asin(body.rotation_matrix[1,0])#cek hack body.Shape.translate(body.h[:nd]) body.barycenter[:] = body.Shape.barycenter body.position[:] = body.Shape.barycenter + #logEvent("6DOF its = " + `its` + " residual = "+`r`) + logEvent("6DOF time {0}".format(body.t)) + logEvent("6DOF DT {0}".format(DT)) + logEvent("6DOF n {0}".format(n)) + logEvent("6DOF force {0}".format(body.FT[1])) + logEvent("displacement, h = {0}".format(body.h)) + logEvent("rotation, Q = {0}".format(body.Q)) + logEvent("velocity, v = {0}".format(body.velocity)) + logEvent("angular acceleration matrix, Omega = {0}".format(body.Omega)) body.step = step #body.scheme = 'Forward_Euler' @@ -406,7 +453,7 @@ def mesh_grading(start, he, grading): # Discretization -- input options useOldPETSc=False -useSuperlu = not True +useSuperlu = True spaceOrder = 1 useHex = False useRBLES = 0.0 @@ -419,15 +466,15 @@ def mesh_grading(start, he, grading): # 3 -- K-Omega, 1988 # Input checks if spaceOrder not in [1,2]: - print "INVALID: spaceOrder" + spaceOrder + print("INVALID: spaceOrder" + spaceOrder) sys.exit() if useRBLES not in [0.0, 1.0]: - print "INVALID: useRBLES" + useRBLES + print("INVALID: useRBLES" + useRBLES) sys.exit() if useMetrics not in [0.0, 1.0]: - print "INVALID: useMetrics" + print("INVALID: useMetrics") sys.exit() # Discretization @@ -435,22 +482,22 @@ def mesh_grading(start, he, grading): if spaceOrder == 1: hFactor=1.0 if useHex: - basis=C0_AffineLinearOnCubeWithNodalBasis - elementQuadrature = CubeGaussQuadrature(nd,3) - elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,3) + basis=C0_AffineLinearOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd,3) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,3) else: - basis=C0_AffineLinearOnSimplexWithNodalBasis - elementQuadrature = SimplexGaussQuadrature(nd,3) - elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,3) - #elementBoundaryQuadrature = SimplexLobattoQuadrature(nd-1,1) + basis=C0_AffineLinearOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,3) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,3) + #elementBoundaryQuadrature = SimplexLobattoQuadrature(nd-1,1) elif spaceOrder == 2: hFactor=0.5 if useHex: - basis=C0_AffineLagrangeOnCubeWithNodalBasis + basis=C0_AffineLagrangeOnCubeWithNodalBasis elementQuadrature = CubeGaussQuadrature(nd,4) elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,4) else: - basis=C0_AffineQuadraticOnSimplexWithNodalBasis + basis=C0_AffineQuadraticOnSimplexWithNodalBasis elementQuadrature = SimplexGaussQuadrature(nd,4) elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,4) @@ -458,8 +505,8 @@ def mesh_grading(start, he, grading): # Numerical parameters if opts.sc == 0.25: sc = 0.25 # default: 0.5. Test: 0.25 - sc_beta = 1. # default: 1.5. Test: 1. - epsFact_consrv_diffusion = 0.1 # default: 1.0. Test: 0.1. Safe: 10. + sc_beta = 1.5 # default: 1.5. Test: 1. + epsFact_consrv_diffusion = 10.0 # default: 1.0. Test: 0.1. Safe: 10. elif opts.sc == 0.5: sc = 0.5 sc_beta = 1.5 @@ -481,9 +528,9 @@ def mesh_grading(start, he, grading): vof_lag_shockCapturing = True vof_sc_uref = 1.0 vof_sc_beta = sc_beta - rd_shockCapturingFactor =sc + rd_shockCapturingFactor = 0.9 rd_lag_shockCapturing = False - epsFact_density = 3. + 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 = epsFact_consrv_diffusion @@ -529,12 +576,12 @@ def mesh_grading(start, he, grading): ns_nl_atol_res = max(1.0e-8,tolfac*he**2) vof_nl_atol_res = max(1.0e-8,tolfac*he**2) ls_nl_atol_res = max(1.0e-8,tolfac*he**2) -mcorr_nl_atol_res = max(1.0e-8,0.1*tolfac*he**2) -rd_nl_atol_res = max(1.0e-8,tolfac*he) +mcorr_nl_atol_res = max(1.0e-8,tolfac*he**2) +rd_nl_atol_res = max(1.0e-8,0.01*he) kappa_nl_atol_res = max(1.0e-8,tolfac*he**2) dissipation_nl_atol_res = max(1.0e-8,tolfac*he**2) mesh_nl_atol_res = max(1.0e-8,mesh_tol*he**2) -am_nl_atol_res = max(1.0e-8,mesh_tol*he**2) +am_nl_atol_res = 0.001#max(1.0e-8,mesh_tol*he**2) #turbulence ns_closure=0 #1-classic smagorinsky, 2-dynamic smagorinsky, 3 -- k-epsilon, 4 -- k-omega diff --git a/2d/floatingStructures/floatingCaissonAddedMass/caisson2D_oscillation_so.py b/2d/floatingStructures/floatingCaissonAddedMass/caisson2D_oscillation_so.py index 19dfb719..2b2b3695 100644 --- a/2d/floatingStructures/floatingCaissonAddedMass/caisson2D_oscillation_so.py +++ b/2d/floatingStructures/floatingCaissonAddedMass/caisson2D_oscillation_so.py @@ -13,7 +13,7 @@ elif '_so.pyc' in name_so[-7:]: name = name_so[:-7] else: - raise NameError, 'Split operator module must end with "_so.py"' + raise(NameError, 'Split operator module must end with "_so.py"') case = __import__(name) Context.setFromModule(case) @@ -47,18 +47,18 @@ #systemStepControllerType = ISO_fixed_MinAdaptiveModelStep if ct.dt_fixed: -# systemStepControllerType = Sequential_FixedStep - systemStepControllerType = Sequential_MinAdaptiveModelStep + systemStepControllerType = Sequential_FixedStep dt_system_fixed = ct.dt_fixed - stepExactSystem=False + systemStepExact=False else: # use CFL systemStepControllerType = Sequential_MinAdaptiveModelStep - stepExactSystem=False + systemStepExact=False needEBQ_GLOBAL = False needEBQ = False modelSpinUpList = [0] # for initial conditions of movemesh +archiveFlag = ArchiveFlags.EVERY_SEQUENCE_STEP if ct.opts.nsave == 0: if ct.dt_fixed > 0: diff --git a/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillation.py b/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillation.py new file mode 100755 index 00000000..a1cda50c --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillation.py @@ -0,0 +1,601 @@ +import math +from proteus import Domain, Context, Comm +from proteus.mprans import SpatialTools as st +from proteus import WaveTools as wt +from proteus.Profiling import logEvent +#import ChRigidBody as crb +from proteus.mbd import CouplingFSI as crb +from math import * +import numpy as np + +opts=Context.Options([ + ("cf",1.0,"Use corrected form of added mass"), + ("nc_form",True,"Use non-conservative NSE form"), + ("use_chrono", False, "use chrono (True) or custom solver"), + # predefined test cases + ("water_level", 1.5, "Height of free surface above bottom"), + # tank + ("tank_dim", (3., 3.,), "Dimensions of the tank"), + ("tank_sponge", (0., 0.), "Length of absorption zones (front/back, left/right)"), + # caisson + ("addedMass", True, "added mass"), + ("caisson", True, "caisson"), + ("caisson_dim", (0.5, 0.5), "Dimensions of the caisson"), + ("caisson_coords", (1.5, 1.5+0.5), "Dimensions of the caisson"), + ("free_x", (0.0, 1.0, 0.0), "Translational DOFs"), + ("free_r", (0.0, 0.0, 0.0), "Rotational DOFs"), + ("VCG", 0.05, "vertical position of the barycenter of the caisson"), + ("caisson_mass", 125., "Mass of the caisson"), + ("caisson_inertia", 4.05, "Inertia of the caisson"), + ("rotation_angle", 0., "Initial rotation angle (in degrees)"), + ("chrono_dt", 0.001, "time step of chrono"), + # mesh refinement + ("refinement", True, "Gradual refinement"), + ("he", 0.03, "Set characteristic element size"), + ("refinement_grading", np.sqrt(1.1*4./np.sqrt(3.))/np.sqrt(1.*4./np.sqrt(3)), "Grading of refinement/coarsening (default: 10% volume)"), + # numerical options + ("genMesh", True, "True: generate new mesh every time. False: do not generate mesh if file exists"), + ("use_gmsh", False, "True: use Gmsh. False: use Triangle/Tetgen"), + ("movingDomain", True, "True/False"), + ("T", 10.0, "Simulation time"), + ("dt_init", 0.001, "Initial time step"), + ("dt_fixed", None, "Fixed (maximum) time step"), + ("timeIntegration", "backwardEuler", "Time integration scheme (backwardEuler/VBDF)"), + ("cfl", 0.4 , "Target cfl"), + ("nsave", 5, "Number of time steps to save per second"), + ("useRANS", 0, "RANS model"), + ("sc", 0.25, "shockCapturing factor"), + ("weak_factor", 10., "weak bc penalty factor"), + ("strong_dir", False, "strong dirichlet (True/False)")]) + + + +# ----- CONTEXT ------ # + +# general options +waterLevel = water_level = opts.water_level +rotation_angle = np.radians(opts.rotation_angle) + + +# tank options +tank_dim = opts.tank_dim +tank_sponge = opts.tank_sponge + +# ----- DOMAIN ----- # + +domain = Domain.PlanarStraightLineGraphDomain() +# caisson options +if opts.caisson is True: + free_x = opts.free_x + free_r = opts.free_r + rotation = np.radians(opts.rotation_angle) + inertia = opts.caisson_inertia + + caisson_dim = opts.caisson_dim + caisson_name='caisson_chrono' + if opts.addedMass: + caisson_name+='_am' + if opts.nc_form: + caisson_name+='_nc' +# caisson = st.Rectangle(domain, dim=opts.caisson_dim, coords=[0.,0.], barycenter=np.zeros(3)) + caisson = st.Circle(domain, radius=opts.caisson_dim[0]/2., coords=[0.,0.], barycenter=np.zeros(3), nPoints=int(math.ceil(2.*math.pi*opts.caisson_dim[0]/opts.he))) + caisson.name=caisson_name + ang = rotation_angle + caisson.setHoles([[0., 0.]]) + caisson.holes_ind = np.array([0]) + print(caisson.regions+np.ones(2)) + print(caisson.regionFlags) + + trans = np.array([opts.caisson_coords[0], opts.caisson_coords[1]]) + print(trans+np.ones(2)) + caisson.translate(trans) + # system = crb.System(np.array([0., -9.81, 0.])) + # rotation = np.array([1, 0., 0., 0.]) + rotation_init = np.array([np.cos(ang/2.), 0., 0., np.sin(ang/2.)*1.]) + caisson.rotate(ang, pivot=caisson.barycenter) + for bc in caisson.BC_list: + bc.setNoSlip() + if opts.use_chrono: + system = crb.ProtChSystem()#np.array([0., -9.81, 0.])) + system.setTimeStep(opts.chrono_dt) + system.step_start = 10 + body = crb.ProtChBody(shape=caisson, + system=system) + chbod = body.ChBody + from proteus.mbd import pyChronoCore as pych + x, y, z = caisson.barycenter + pos = pych.ChVector(x, y, z) + e0, e1, e2, e3 = rotation_init + rot = pych.ChQuaternion(e0, e1, e2, e3) + inertia = pych.ChVector(1., 1., inertia) + chbod.SetPos(pos) + chbod.SetRot(rot) + chbod.SetMass(opts.caisson_mass) + chbod.SetInertiaXX(inertia) + body.setConstraints(free_x=np.array(opts.free_x), free_r=np.array(opts.free_r)) + # body.setInitialRot(rotation_init) + # body.rotation_init=np.array([np.cos(ang/2.), 0., 0., np.sin(ang/2.)*1.]) + body.setRecordValues(all_values=True) + else: + from proteus import AuxiliaryVariables + class BodyAddedMass(AuxiliaryVariables.AV_base): + def __init__(self): + AuxiliaryVariables.AV_base.__init__(self) + + def attachModel(self, model, ar): + self.model=model + return self + + def calculate(self): + pass + from proteus.mprans import BodyDynamics as bd + body = bd.RigidBody(shape=caisson) + system = body # so it is added in twp_navier_stokes_p.py + bodyAddedMass = BodyAddedMass() + body.bodyAddedMass = body.ProtChAddedMass = bodyAddedMass + body.setMass(opts.caisson_mass) + body.setConstraints(free_x=free_x, + free_r=free_r) + body.It = np.array([[1., 0., 0.], + [0., 1., 0.], + [0., 0., inertia]]) + filename='record_caisson' + if opts.addedMass: + filename+='_am' + if opts.nc_form: + filename+='_nc' + if opts.cf: + filename+='_cf' + body.setRecordValues(filename=filename, all_values=True) + body.coords_system = caisson.coords_system # hack + body.last_Aij = np.zeros((6,6),'d') + body.last_a = np.zeros((6,),'d') + body.a = np.zeros((6,),'d') + body.last_Omega = np.zeros((3,3),'d') + body.last_velocity = np.zeros((3,),'d') + body.last_Q = np.eye(3) + body.last_h = np.zeros((3,),'d') + body.last_mom = np.zeros((6,),'d') + body.last_u = np.zeros((18,),'d') + body.last_t = 0.0 + body.t = 0.0 + body.h = np.zeros((3,),'d') + body.velocity = np.zeros((3,),'d') + body.mom = np.zeros((6,),'d') + body.Omega = np.zeros((3,3),'d') + body.Q = np.eye(3,3) + body.u = np.zeros((18,),'d') + body.u[9:18]= body.Q.flatten() + body.last_u[:] = body.u + body.last_Q[:] = body.Q + body.FT = np.zeros((6,),'d') + body.last_FT = np.zeros((6,),'d') + body.free_dof = np.zeros((6,),'d') + body.free_dof[:3] = free_x + body.free_dof[3:] = free_r + # cek -- I did this with a different hack, see dummy AddedMassBody above + # body.Aij = np.zeros((6,6)) + + # # we create a ProtChAddedMass auxiliary variable to add to added_mass_n.py + # # the argument here should be a ProtChSystem, so we give it an empty system (hack) + # # this is just to get access to the AddedMass model from the body functions + # system = crb.ProtChSystem(g) + # body.ProtChAddedMass = crb.ProtChAddedMass(system) + + def step(dt): + from math import ceil + logEvent("Barycenter "+str(body.barycenter)) + n = max(1.0,ceil(dt/opts.chrono_dt)) + DT=dt/n + def F(u,theta): + """The residual and Jacobian for discrete 6DOF motion with added mass""" + v = u[:3] + omega = u[3:6] + h = u[6:9] + Q = u[9:18].reshape(3,3) + Omega = np.array([[ 0.0, -omega[2], omega[1]], + [ omega[2], 0.0, -omega[0]], + [-omega[1], omega[0], 0.0]]) + I = np.matmul(np.matmul(Q, body.It), Q.transpose()) + body.Aij = np.zeros((6,6),'d') + if opts.addedMass: + for i in range(1,5):#number of rigid body facets + body.Aij += body.bodyAddedMass.model.levelModelList[-1].Aij[i] + avg_Aij=False + if avg_Aij: + M = body.Aij*theta + body.last_Aij*(1-theta) + else: + M = body.Aij.copy() + for i in range(6): + for j in range(6): + M[i,j]*=body.free_dof[j]#only allow j accelerations to contribute to i force balance if j is free + M[j,i]*=body.free_dof[j]#only allow j added mass forces if j is free + body.Aij[i,j]*=body.free_dof[j]#only allow j accelerations to contribute to i force balance if j is free + body.Aij[j,i]*=body.free_dof[j]#only allow j added mass forces if j is free + body.FT[:3] = body.F + body.FT[3:] = body.M + #cek debug + #body.FT[:] = 0.0 + #body.FT[1] = body.mass * 0.125*math.pi**2 * math.cos(body.t*math.pi) + for i in range(3): + M[i, i] += body.mass + for j in range(3): + M[3+i, 3+j] += I[i, j] + r = np.zeros((18,),'d') + BE=True + CF=opts.cf + if BE: + r[:6] = np.matmul(M, u[:6]) - np.matmul(body.Aij, body.last_u[:6]) - body.last_mom - DT*body.FT - CF*np.matmul(body.Aij,body.last_a) + r[6:9] = h - body.last_h - DT*v + rQ = Q - body.last_Q - DT*np.matmul(Omega,Q) + else: + r[:6] = np.matmul(M, u[:6]) - np.matmul(body.Aij, body.last_u[:6]) - body.last_mom - DT*(body.FT*theta+body.last_FT*(1.0-theta)) - CF*np.matmul(body.Aij,body.last_a) + r[6:9] = h - body.last_h - DT*0.5*(v + body.last_velocity) + rQ = Q - body.last_Q - DT*0.25*np.matmul((Omega + body.last_Omega),(Q+body.last_Q)) + r[9:18] = rQ.flatten() + J = np.zeros((18,18),'d') + J[:6,:6] = M + #neglecting 0:6 dependence on Q + for i in range(3): + if BE: + J[6+i,i] = -DT + else: + J[6+i,i] = -DT*0.5 + J[6+i,6+i] = 1.0 + for i in range(9): + J[9+i,9+i] = 1.0 + for i in range(3): + for j in range(3): + if BE: + J[9+i*3+j, 9+i+j*3] -= DT*Omega[i,j] + else: + J[9+i*3+j, 9+i+j*3] -= DT*0.25*(Omega+body.last_Omega)[i,j] + #neglecting 9:18 dependence on omega + body.Omega[:] = Omega + body.velocity[:] = v + body.Q[:] = Q + body.h[:] = h + body.u[:] = u + body.mom[:3] = body.mass*u[:3] + body.mom[3:6] = np.matmul(I,u[3:6]) + body.a[:] = u[:6] - body.last_u[:6] + return r, J + nd = body.nd + Q_start=body.Q.copy() + h_start=body.last_h.copy() + for i in range(int(n)): + theta = (i+1)*DT/dt + body.t = body.last_t + DT + logEvent("6DOF theta {0}".format(theta)) + u = np.zeros((18,),'d') + u[:] = body.last_u + r = np.zeros((18,),'d') + r,J = F(u,theta) + its=0 + maxits=100 + while ((its==0 or np.linalg.norm(r) > 1.0e-8) and its < maxits): + u -= np.linalg.solve(J,r) + r,J = F(u,theta) + its+=1 + logEvent("6DOF its {0}".format(its)) + logEvent("6DOF res {0}".format(np.linalg.norm(r))) + body.last_Aij[:]=body.Aij + body.last_FT[:] = body.FT + body.last_Omega[:] = body.Omega + body.last_velocity[:] = body.velocity + body.last_Q[:] = body.Q + body.last_h[:] = body.h + body.last_u[:] = body.u + body.last_mom[:] = body.mom + body.last_t = body.t + body.last_a[:] = body.a + # translate and rotate + body.h -= h_start + body.last_h[:] = 0.0 + body.last_position[:] = body.position + #body.rotation_matrix[:] = np.linalg.solve(Q_start,body.Q) + body.rotation_matrix[:] = np.matmul(np.linalg.inv(body.Q),Q_start) + body.rotation_euler[2] -= math.asin(body.rotation_matrix[1,0])#cek hack + body.Shape.translate(body.h[:nd]) + body.barycenter[:] = body.Shape.barycenter + body.position[:] = body.Shape.barycenter + #logEvent("6DOF its = " + `its` + " residual = "+`r`) + logEvent("6DOF time {0}".format(body.t)) + logEvent("6DOF DT {0}".format(DT)) + logEvent("6DOF n {0}".format(n)) + logEvent("6DOF force {0}".format(body.FT[1])) + logEvent("displacement, h = {0}".format(body.h)) + logEvent("rotation, Q = {0}".format(body.Q)) + logEvent("velocity, v = {0}".format(body.velocity)) + logEvent("angular acceleration matrix, Omega = {0}".format(body.Omega)) + + body.step = step + #body.scheme = 'Forward_Euler' + +# ----- SHAPES ----- # +tank = st.Tank2D(domain, tank_dim) +tank.setSponge(x_n=tank_sponge[0], x_p=tank_sponge[1]) +if opts.caisson: + # let gmsh know that the caisson is IN the tank + tank.setChildShape(caisson, 0) + + +# ----- BOUNDARY CONDITIONS ----- # + +tank.BC['y+'].setAtmosphere() +tank.BC['y-'].setNoSlip() +tank.BC['x+'].setNoSlip() +tank.BC['x-'].setNoSlip() +tank.BC['sponge'].setNonMaterial() + +tank.BC['x-'].setFixedNodes() +tank.BC['x+'].setFixedNodes() +tank.BC['sponge'].setFixedNodes() +tank.BC['y+'].setTank() # sliding mesh nodes +tank.BC['y-'].setTank() #sliding mesh nodes + + + + +domain.MeshOptions.use_gmsh = opts.use_gmsh +domain.MeshOptions.genMesh = opts.genMesh +he = opts.he +domain.MeshOptions.he = he +st.assembleDomain(domain) +domain.use_gmsh = opts.use_gmsh +geofile='mesh'+str(opts.he) +domain.geofile=geofile + + +# MESH REFINEMENT + +if opts.use_gmsh: + import py2gmsh + from MeshRefinement import geometry_to_gmsh + mesh = geometry_to_gmsh(domain) + grading = opts.refinement_grading + he = opts.he + he_max = 10. + ecH = 3. + field_list = [] + + def mesh_grading(start, he, grading): + return '{he}*{grading}^(1+log((-1/{grading}*(abs({start})-{he})+abs({start}))/{he})/log({grading}))'.format(he=he, start=start, grading=grading) + + me01 = py2gmsh.Fields.MathEval(mesh=mesh) + dist = '(abs({zcoord}-y))'.format(zcoord=water_level) + me01.F = mesh_grading(he=he, start=dist, grading=grading) + field_list += [me01] + + me3 = py2gmsh.Fields.MathEval(mesh=mesh) + dist_z = '(abs(abs({z_p}-y)+abs(y-{z_n})-({z_p}-{z_n}))/2.)'.format(z_p=max(caisson.vertices[:,1]), z_n=min(caisson.vertices[:,1])) + dist_x = '(abs(abs({z_p}-x)+abs(x-{z_n})-({z_p}-{z_n}))/2.)'.format(z_p=max(caisson.vertices[:,0]), z_n=min(caisson.vertices[:,0])) + me3.F = '{he}*{grading}^(Sqrt({dist_x}^2+{dist_z}^2)/{he})'.format(he=he, grading=grading, dist_x=dist_x, dist_z=dist_z) + me3.F = mesh_grading(he=he, start=dist, grading=grading) + field_list += [me3] + + # background field + fmin = py2gmsh.Fields.Min(mesh=mesh) + fmin.FieldsList = field_list + mesh.setBackgroundField(fmin) + + # max element size + mesh.Options.Mesh.CharacteristicLengthMax = he_max + + mesh.writeGeo(geofile+'.geo') + + + +# passed in added_mass_p.py coefficients +max_flag = 0 +max_flag = max(domain.vertexFlags) +max_flag = max(domain.segmentFlags+[max_flag]) +max_flag = max(domain.facetFlags+[max_flag]) +flags_rigidbody = np.zeros(max_flag+1, dtype='int32') +if opts.use_chrono: + for s in system.subcomponents: + if type(s) is crb.ProtChBody: + for i in range(s.i_start, s.i_end): + flags_rigidbody[i] = 1 +else: + flags_rigidbody[1:5] = 1 + + +########################################## +# Numerical Options and other parameters # +########################################## + + +rho_0=998.2 +nu_0 =1.004e-6 +rho_1=1.205 +nu_1 =1.500e-5 +sigma_01=0.0 +g = [0., -9.81] + + + + +from math import * +from proteus import MeshTools, AuxiliaryVariables +import numpy +import proteus.MeshTools +from proteus import Domain +from proteus.Profiling import logEvent +from proteus.default_n import * +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.ctransportCoefficients import smoothedHeaviside_integral + + +#---------------------------------------------------- +# Boundary conditions and other flags +#---------------------------------------------------- +movingDomain=opts.movingDomain +checkMass=False +applyCorrection=True +applyRedistancing=True +freezeLevelSet=True + +#---------------------------------------------------- +# Time stepping and velocity +#---------------------------------------------------- +weak_bc_penalty_constant = opts.weak_factor/nu_0#Re +dt_init = opts.dt_init +T = opts.T +nDTout = int(opts.T*opts.nsave) +timeIntegration = opts.timeIntegration +if nDTout > 0: + dt_out= (T-dt_init)/nDTout +else: + dt_out = 0 +runCFL = opts.cfl +dt_fixed = opts.dt_fixed + +#---------------------------------------------------- + +# Discretization -- input options +useOldPETSc=False +useSuperlu = True +spaceOrder = 1 +useHex = False +useRBLES = 0.0 +useMetrics = 1.0 +useVF = 1.0 +useOnlyVF = False +useRANS = opts.useRANS # 0 -- None + # 1 -- K-Epsilon + # 2 -- K-Omega, 1998 + # 3 -- K-Omega, 1988 +# 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,3) + elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,3) + else: + basis=C0_AffineLinearOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd,3) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,3) + #elementBoundaryQuadrature = SimplexLobattoQuadrature(nd-1,1) +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) + + +# Numerical parameters +if opts.sc == 0.25: + sc = 0.25 # default: 0.5. Test: 0.25 + sc_beta = 1.5 # default: 1.5. Test: 1. + epsFact_consrv_diffusion = 10.0 # default: 1.0. Test: 0.1. Safe: 10. +elif opts.sc == 0.5: + sc = 0.5 + sc_beta = 1.5 + epsFact_consrv_diffusion = 10.0 # default: 1.0. Test: 0.1. Safe: 10. +else: + import sys + sys.quit() +ns_forceStrongDirichlet = opts.strong_dir +backgroundDiffusionFactor=0.01 +if useMetrics: + ns_shockCapturingFactor = sc + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = sc + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = sc_beta + vof_shockCapturingFactor = sc + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = sc_beta + 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 = epsFact_consrv_diffusion + redist_Newton = True#False + kappa_shockCapturingFactor = sc + kappa_lag_shockCapturing = False#True + kappa_sc_uref = 1.0 + kappa_sc_beta = sc_beta + dissipation_shockCapturingFactor = sc + dissipation_lag_shockCapturing = False#True + dissipation_sc_uref = 1.0 + dissipation_sc_beta = sc_beta +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#True + 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 + +tolfac = 0.001 +mesh_tol = 0.001 +ns_nl_atol_res = max(1.0e-8,tolfac*he**2) +vof_nl_atol_res = max(1.0e-8,tolfac*he**2) +ls_nl_atol_res = max(1.0e-8,tolfac*he**2) +mcorr_nl_atol_res = max(1.0e-8,tolfac*he**2) +rd_nl_atol_res = max(1.0e-8,0.01*he) +kappa_nl_atol_res = max(1.0e-8,tolfac*he**2) +dissipation_nl_atol_res = max(1.0e-8,tolfac*he**2) +mesh_nl_atol_res = max(1.0e-8,mesh_tol*he**2) +am_nl_atol_res = 0.001#max(1.0e-8,mesh_tol*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 + +def twpflowPressure_init(x, t): + p_L = 0.0 + phi_L = tank_dim[nd-1] - waterLevel + phi = x[nd-1] - waterLevel + return p_L -g[nd-1]*(rho_0*(phi_L - phi)+(rho_1 -rho_0)*(smoothedHeaviside_integral(epsFact_consrv_heaviside*opts.he,phi_L) + -smoothedHeaviside_integral(epsFact_consrv_heaviside*opts.he,phi))) diff --git a/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillationIB.py b/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillationIB.py new file mode 100755 index 00000000..a836c6e8 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillationIB.py @@ -0,0 +1,611 @@ +import math +from proteus import Domain, Context, Comm +from proteus.mprans import SpatialTools as st +from proteus import WaveTools as wt +from proteus.Profiling import logEvent +#import ChRigidBody as crb +from proteus.mbd import CouplingFSI as crb +from math import * +import numpy as np + +opts=Context.Options([ + ("cf",1.0,"Use corrected form of added mass"), + ("nc_form",True,"Use non-conservative NSE form"), + ("use_chrono", False, "use chrono (True) or custom solver"), + # predefined test cases + ("water_level", 1.5, "Height of free surface above bottom"), + # tank + ("tank_dim", (3., 3.,), "Dimensions of the tank"), + ("tank_sponge", (0., 0.), "Length of absorption zones (front/back, left/right)"), + # caisson + ("addedMass", True, "added mass"), + ("caisson", True, "caisson"), + ("caisson_dim", (0.5, 0.5), "Dimensions of the caisson"), + ("caisson_coords", (1.5, 1.5+0.5), "Dimensions of the caisson"), + ("free_x", (0.0, 1.0, 0.0), "Translational DOFs"), + ("free_r", (0.0, 0.0, 0.0), "Rotational DOFs"), + ("VCG", 0.05, "vertical position of the barycenter of the caisson"), + ("caisson_mass", 125., "Mass of the caisson"), + ("caisson_inertia", 4.05, "Inertia of the caisson"), + ("rotation_angle", 0., "Initial rotation angle (in degrees)"), + ("chrono_dt", 0.0001, "time step of chrono"), + # mesh refinement + ("refinement", True, "Gradual refinement"), + ("he", 0.03, "Set characteristic element size"), + ("refinement_grading", np.sqrt(1.1*4./np.sqrt(3.))/np.sqrt(1.*4./np.sqrt(3)), "Grading of refinement/coarsening (default: 10% volume)"), + # numerical options + ("genMesh", True, "True: generate new mesh every time. False: do not generate mesh if file exists"), + ("use_gmsh", False, "True: use Gmsh. False: use Triangle/Tetgen"), + ("movingDomain", False, "True/False"), + ("T", 10.0, "Simulation time"), + ("dt_init", 0.001, "Initial time step"), + ("dt_fixed",0.001, "Fixed (maximum) time step"), + ("timeIntegration", "backwardEuler", "Time integration scheme (backwardEuler/VBDF)"), + ("cfl", 0.4 , "Target cfl"), + ("nsave", 10, "Number of time steps to save per second"), + ("useRANS", 0, "RANS model"), + ("sc", 0.5, "shockCapturing factor"), + ("weak_factor", 10., "weak bc penalty factor"), + ("strong_dir", False, "strong dirichlet (True/False)")]) + +eps=1.0e-8 +def onBoundary(x): + if(x[0] < 0.0 + eps or + x[0] > opts.tank_dim[0] - eps or + x[1] < 0.0 + eps or + x[1] > opts.tank_dim[1] - eps): + return True + else: + return False + +# ----- CONTEXT ------ # + +# general options +waterLevel = water_level = opts.water_level +rotation_angle = np.radians(opts.rotation_angle) + + +# tank options +tank_dim = opts.tank_dim +tank_sponge = opts.tank_sponge + +# ----- DOMAIN ----- # + +domain = Domain.PlanarStraightLineGraphDomain() +domainDummy = Domain.PlanarStraightLineGraphDomain() +# caisson options +use_ball_as_particle = True +ball_center=np.zeros((1,3),'d') +ball_radius=np.array([opts.caisson_dim[0]/2.],'d') +ball_velocity=np.zeros((1,3),'d') +ball_angular_velocity=np.zeros((1,3),'d') +if opts.caisson is True: + free_x = opts.free_x + free_r = opts.free_r + rotation = np.radians(opts.rotation_angle) + inertia = opts.caisson_inertia + + caisson_dim = opts.caisson_dim + caisson_name='caisson_chrono' + if opts.addedMass: + caisson_name+='_am' + if opts.nc_form: + caisson_name+='_nc' +# caisson = st.Rectangle(domain, dim=opts.caisson_dim, coords=[0.,0.], barycenter=np.zeros(3)) + caisson = st.Circle(domainDummy, radius=opts.caisson_dim[0]/2., coords=[0.,0.], barycenter=np.zeros(3), nPoints=int(math.ceil(2.*math.pi*opts.caisson_dim[0]/opts.he))) + assert(opts.caisson_mass/math.pi*(opts.caisson_dim[0]/2.)**2 < 998.2) + caisson.name=caisson_name + ang = rotation_angle + #caisson.setHoles([[0., 0.]]) + #caisson.holes_ind = np.array([0]) + #print(caisson.regions+np.ones(2)) + #print(caisson.regionFlags) + + trans = np.array([opts.caisson_coords[0], opts.caisson_coords[1]]) + print(trans+np.ones(2)) + caisson.translate(trans) + # system = crb.System(np.array([0., -9.81, 0.])) + # rotation = np.array([1, 0., 0., 0.]) + rotation_init = np.array([np.cos(ang/2.), 0., 0., np.sin(ang/2.)*1.]) + caisson.rotate(ang, pivot=caisson.barycenter) +# for bc in caisson.BC_list: +# bc.setNoSlip() + if opts.use_chrono: + system = crb.ProtChSystem(np.array([0., -9.81, 0.])) + system.setTimeStep(opts.chrono_dt) + system.step_start = 10 + body = crb.ProtChBody(shape=caisson, + system=system) + chbod = body.ChBody + from proteus.mbd import pyChronoCore as pych + x, y, z = caisson.barycenter + pos = pych.ChVector(x, y, z) + e0, e1, e2, e3 = rotation_init + rot = pych.ChQuaternion(e0, e1, e2, e3) + inertia = pych.ChVector(1., 1., inertia) + chbod.SetPos(pos) + chbod.SetRot(rot) + chbod.SetMass(opts.caisson_mass) + chbod.SetInertiaXX(inertia) + body.setConstraints(free_x=np.array(opts.free_x), free_r=np.array(opts.free_r)) + # body.setInitialRot(rotation_init) + # body.rotation_init=np.array([np.cos(ang/2.), 0., 0., np.sin(ang/2.)*1.]) + body.setRecordValues(all_values=True) + else: + from proteus import AuxiliaryVariables + class BodyAddedMass(AuxiliaryVariables.AV_base): + def __init__(self): + AuxiliaryVariables.AV_base.__init__(self) + + def attachModel(self, model, ar): + self.model=model + return self + + def calculate(self): + pass + from proteus.mprans import BodyDynamics as bd + body = bd.RigidBody(shape=caisson) + ball_center[0,:]=body.Shape.barycenter + system = body # so it is added in twp_navier_stokes_p.py + bodyAddedMass = BodyAddedMass() + body.bodyAddedMass = body.ProtChAddedMass = bodyAddedMass + body.setMass(opts.caisson_mass) + body.setConstraints(free_x=free_x, + free_r=free_r) + body.It = np.array([[1., 0., 0.], + [0., 1., 0.], + [0., 0., inertia]]) + filename='record_caisson' + if opts.addedMass: + filename+='_am' + if opts.nc_form: + filename+='_nc' + if opts.cf: + filename+='_cf' + body.setRecordValues(filename=filename, all_values=True) + body.coords_system = caisson.coords_system # hack + body.last_Aij = np.zeros((6,6),'d') + body.last_a = np.zeros((6,),'d') + body.a = np.zeros((6,),'d') + body.last_Omega = np.zeros((3,3),'d') + body.last_velocity = np.zeros((3,),'d') + body.last_Q = np.eye(3) + body.last_h = np.zeros((3,),'d') + body.last_mom = np.zeros((6,),'d') + body.last_u = np.zeros((18,),'d') + body.last_t = 0.0 + body.t = 0.0 + body.h = np.zeros((3,),'d') + body.velocity = np.zeros((3,),'d') + body.mom = np.zeros((6,),'d') + body.Omega = np.zeros((3,3),'d') + body.Q = np.eye(3,3) + body.u = np.zeros((18,),'d') + body.u[9:18]= body.Q.flatten() + body.last_u[:] = body.u + body.last_Q[:] = body.Q + body.FT = np.zeros((6,),'d') + body.last_FT = np.zeros((6,),'d') + body.free_dof = np.zeros((6,),'d') + body.free_dof[:3] = free_x + body.free_dof[3:] = free_r + # cek -- I did this with a different hack, see dummy AddedMassBody above + # body.Aij = np.zeros((6,6)) + + # # we create a ProtChAddedMass auxiliary variable to add to added_mass_n.py + # # the argument here should be a ProtChSystem, so we give it an empty system (hack) + # # this is just to get access to the AddedMass model from the body functions + # system = crb.ProtChSystem(g) + # body.ProtChAddedMass = crb.ProtChAddedMass(system) + + def step(dt): + from math import ceil + logEvent("Barycenter "+str(body.barycenter)) + n = max(1.0,ceil(dt/opts.chrono_dt)) + DT=dt/n + def F(u,theta): + """The residual and Jacobian for discrete 6DOF motion with added mass""" + v = u[:3] + omega = u[3:6] + h = u[6:9] + Q = u[9:18].reshape(3,3) + Omega = np.array([[ 0.0, -omega[2], omega[1]], + [ omega[2], 0.0, -omega[0]], + [-omega[1], omega[0], 0.0]]) + I = np.matmul(np.matmul(Q, body.It), Q.transpose()) + body.Aij = np.zeros((6,6),'d') + if opts.addedMass: + body.Aij += body.bodyAddedMass.model.levelModelList[-1].coefficients.particle_Aij[0] + avg_Aij=False + if avg_Aij: + M = body.Aij*theta + body.last_Aij*(1-theta) + else: + M = body.Aij.copy() + for i in range(6): + for j in range(6): + M[i,j]*=body.free_dof[j]#only allow j accelerations to contribute to i force balance if j is free + M[j,i]*=body.free_dof[j]#only allow j added mass forces if j is free + body.Aij[i,j]*=body.free_dof[j]#only allow j accelerations to contribute to i force balance if j is free + body.Aij[j,i]*=body.free_dof[j]#only allow j added mass forces if j is free + body.FT[:3] = opts.caisson_mass*np.array([0.,-9.8,0.]) + body.bodyAddedMass.model.levelModelList[-1].coefficients.flowModel.coefficients.particle_netForces[0] + body.FT[3:] = body.bodyAddedMass.model.levelModelList[-1].coefficients.flowModel.coefficients.particle_netMoments[0] + #cek debug + #body.FT[:] = 0.0 + #body.FT[1] = body.mass * 0.125*math.pi**2 * math.cos(body.t*math.pi) + for i in range(3): + M[i, i] += body.mass + for j in range(3): + M[3+i, 3+j] += I[i, j] + r = np.zeros((18,),'d') + BE=True + CF=opts.cf + if BE: + r[:6] = np.matmul(M, u[:6]) - np.matmul(body.Aij, body.last_u[:6]) - body.last_mom - DT*body.FT - CF*np.matmul(body.Aij,body.last_a) + r[6:9] = h - body.last_h - DT*v + rQ = Q - body.last_Q - DT*np.matmul(Omega,Q) + else: + r[:6] = np.matmul(M, u[:6]) - np.matmul(body.Aij, body.last_u[:6]) - body.last_mom - DT*(body.FT*theta+body.last_FT*(1.0-theta)) - CF*np.matmul(body.Aij,body.last_a) + r[6:9] = h - body.last_h - DT*0.5*(v + body.last_velocity) + rQ = Q - body.last_Q - DT*0.25*np.matmul((Omega + body.last_Omega),(Q+body.last_Q)) + r[9:18] = rQ.flatten() + J = np.zeros((18,18),'d') + J[:6,:6] = M + #neglecting 0:6 dependence on Q + for i in range(3): + if BE: + J[6+i,i] = -DT + else: + J[6+i,i] = -DT*0.5 + J[6+i,6+i] = 1.0 + for i in range(9): + J[9+i,9+i] = 1.0 + for i in range(3): + for j in range(3): + if BE: + J[9+i*3+j, 9+i+j*3] -= DT*Omega[i,j] + else: + J[9+i*3+j, 9+i+j*3] -= DT*0.25*(Omega+body.last_Omega)[i,j] + #neglecting 9:18 dependence on omega + body.Omega[:] = Omega + body.velocity[:] = v + body.Q[:] = Q + body.h[:] = h + body.u[:] = u + body.mom[:3] = body.mass*u[:3] + body.mom[3:6] = np.matmul(I,u[3:6]) + body.a[:] = u[:6] - body.last_u[:6] + return r, J + nd = body.nd + Q_start=body.Q.copy() + h_start=body.last_h.copy() + for i in range(int(n)): + theta = (i+1)*DT/dt + body.t = body.last_t + DT + logEvent("6DOF theta {0}".format(theta)) + u = np.zeros((18,),'d') + u[:] = body.last_u + r = np.zeros((18,),'d') + r,J = F(u,theta) + its=0 + maxits=100 + while ((its==0 or np.linalg.norm(r) > 1.0e-8) and its < maxits): + u -= np.linalg.solve(J,r) + r,J = F(u,theta) + its+=1 + logEvent("6DOF its {0}".format(its)) + logEvent("6DOF res {0}".format(np.linalg.norm(r))) + body.last_Aij[:]=body.Aij + body.last_FT[:] = body.FT + body.last_Omega[:] = body.Omega + body.last_velocity[:] = body.velocity + body.last_Q[:] = body.Q + body.last_h[:] = body.h + body.last_u[:] = body.u + body.last_mom[:] = body.mom + body.last_t = body.t + body.last_a[:] = body.a + # translate and rotate + body.h -= h_start + body.last_h[:] = 0.0 + body.last_position[:] = body.position + #body.rotation_matrix[:] = np.linalg.solve(Q_start,body.Q) + body.rotation_matrix[:] = np.matmul(np.linalg.inv(body.Q),Q_start) + body.rotation_euler[2] -= math.asin(body.rotation_matrix[1,0])#cek hack + body.Shape.translate(body.h[:nd]) + body.barycenter[:] = body.Shape.barycenter + body.position[:] = body.Shape.barycenter + ball_center[0,:]=body.Shape.barycenter + body.bodyAddedMass.model.levelModelList[-1].coefficients.flowModel.coefficients.ball_center[0,:]=ball_center + #ball_radius=np.zeros([opts.caisson_dim[0]/2.],'d') + ball_velocity[0,:]=body.velocity + body.bodyAddedMass.model.levelModelList[-1].coefficients.flowModel.coefficients.ball_radius[0]=opts.caisson_dim[0]/2. + body.bodyAddedMass.model.levelModelList[-1].coefficients.flowModel.coefficients.ball_velocity[0,:]=body.velocity + body.bodyAddedMass.model.levelModelList[-1].coefficients.flowModel.coefficients.ball_angular_velocity[0,:]=0.0 + #ball_angular_velocity[0,:]=#hack! + #logEvent("6DOF its = " + `its` + " residual = "+`r`) + logEvent("6DOF time {0}".format(body.t)) + logEvent("6DOF DT {0}".format(DT)) + logEvent("6DOF n {0}".format(n)) + logEvent("6DOF force {0}".format(body.FT[1])) + logEvent("displacement, h = {0}".format(body.h)) + logEvent("rotation, Q = {0}".format(body.Q)) + logEvent("velocity, v = {0}".format(body.velocity)) + logEvent("angular acceleration matrix, Omega = {0}".format(body.Omega)) + + body.step = step + #body.scheme = 'Forward_Euler' + +# ----- SHAPES ----- # +tank = st.Tank2D(domain, tank_dim) +#tank.setSponge(x_n=tank_sponge[0], x_p=tank_sponge[1]) +#if opts.caisson: +# # let gmsh know that the caisson is IN the tank +# tank.setChildShape(caisson, 0) + + +# ----- BOUNDARY CONDITIONS ----- # + +tank.BC['y+'].setAtmosphere() +tank.BC['y-'].setNoSlip() +tank.BC['x+'].setNoSlip() +tank.BC['x-'].setNoSlip() +#tank.BC['sponge'].setNonMaterial() + +tank.BC['x-'].setFixedNodes() +tank.BC['x+'].setFixedNodes() +#tank.BC['sponge'].setFixedNodes() +tank.BC['y+'].setTank() # sliding mesh nodes +tank.BC['y-'].setTank() #sliding mesh nodes +boundaryTags={'left':4,'right':2,'top':3,'bottom':1,'front':1,'back':3} + + + +domain.MeshOptions.use_gmsh = opts.use_gmsh +domain.MeshOptions.genMesh = opts.genMesh +he = opts.he +domain.MeshOptions.he = he +st.assembleDomain(domain) +domain.use_gmsh = opts.use_gmsh +geofile='mesh'+str(opts.he) +domain.geofile=geofile + + +# MESH REFINEMENT + +if opts.use_gmsh: + import py2gmsh + from MeshRefinement import geometry_to_gmsh + mesh = geometry_to_gmsh(domain) + grading = opts.refinement_grading + he = opts.he + he_max = 10. + ecH = 3. + field_list = [] + + def mesh_grading(start, he, grading): + return '{he}*{grading}^(1+log((-1/{grading}*(abs({start})-{he})+abs({start}))/{he})/log({grading}))'.format(he=he, start=start, grading=grading) + + me01 = py2gmsh.Fields.MathEval(mesh=mesh) + dist = '(abs({zcoord}-y))'.format(zcoord=water_level) + me01.F = mesh_grading(he=he, start=dist, grading=grading) + field_list += [me01] + + me3 = py2gmsh.Fields.MathEval(mesh=mesh) + dist_z = '(abs(abs({z_p}-y)+abs(y-{z_n})-({z_p}-{z_n}))/2.)'.format(z_p=max(caisson.vertices[:,1]), z_n=min(caisson.vertices[:,1])) + dist_x = '(abs(abs({z_p}-x)+abs(x-{z_n})-({z_p}-{z_n}))/2.)'.format(z_p=max(caisson.vertices[:,0]), z_n=min(caisson.vertices[:,0])) + me3.F = '{he}*{grading}^(Sqrt({dist_x}^2+{dist_z}^2)/{he})'.format(he=he, grading=grading, dist_x=dist_x, dist_z=dist_z) + me3.F = mesh_grading(he=he, start=dist, grading=grading) + field_list += [me3] + + # background field + fmin = py2gmsh.Fields.Min(mesh=mesh) + fmin.FieldsList = field_list + mesh.setBackgroundField(fmin) + + # max element size + mesh.Options.Mesh.CharacteristicLengthMax = he_max + + mesh.writeGeo(geofile+'.geo') + + + +# passed in added_mass_p.py coefficients +max_flag = 0 +max_flag = max(domain.vertexFlags) +max_flag = max(domain.segmentFlags+[max_flag]) +max_flag = max(domain.facetFlags+[max_flag]) +flags_rigidbody = np.zeros(max_flag+1, dtype='int32') +#if opts.use_chrono: +# for s in system.subcomponents: +# if type(s) is crb.ProtChBody: +# for i in range(s.i_start, s.i_end): +# flags_rigidbody[i] = 1 +#else: +# flags_rigidbody[1:5] = 1 + + +########################################## +# Numerical Options and other parameters # +########################################## + + +rho_0=998.2 +nu_0 =1.004e-6 +rho_1=1.205 +nu_1 =1.500e-5 +#cek debug +#rho_1=rho_0 +#nu_1 =nu_0 +sigma_01=0.0 +g = [0., -9.81] + + + + +from math import * +from proteus import MeshTools, AuxiliaryVariables +import numpy +import proteus.MeshTools +from proteus import Domain +from proteus.Profiling import logEvent +from proteus.default_n import * +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.ctransportCoefficients import smoothedHeaviside_integral + + +#---------------------------------------------------- +# Boundary conditions and other flags +#---------------------------------------------------- +movingDomain=opts.movingDomain +checkMass=False +applyCorrection=True +applyRedistancing=True +freezeLevelSet=True + +#---------------------------------------------------- +# Time stepping and velocity +#---------------------------------------------------- +weak_bc_penalty_constant = opts.weak_factor/nu_0#Re +dt_init = opts.dt_init +T = opts.T +nDTout = int(opts.T*opts.nsave) +timeIntegration = opts.timeIntegration +if nDTout > 0: + dt_out= (T-dt_init)/nDTout +else: + dt_out = 0 +runCFL = opts.cfl +dt_fixed = opts.dt_fixed + +#---------------------------------------------------- + +# Discretization -- input options +useOldPETSc=False +useSuperlu = True +spaceOrder = 1 +useHex = False +useRBLES = 0.0 +useMetrics = 1.0 +useVF = 1.0 +useOnlyVF = False +useRANS = opts.useRANS # 0 -- None + # 1 -- K-Epsilon + # 2 -- K-Omega, 1998 + # 3 -- K-Omega, 1988 +# 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 +hFactor=1.0 +pbasis=basis=C0_AffineLinearOnSimplexWithNodalBasis +vbasis=C0_AffineQuadraticOnSimplexWithNodalBasis +elementQuadrature = SimplexGaussQuadrature(nd,5) +elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,5) + + +# Numerical parameters +if opts.sc == 0.25: + sc = 0.25 # default: 0.5. Test: 0.25 + sc_beta = 1.5 # default: 1.5. Test: 1. + epsFact_consrv_diffusion = 10.0 # default: 1.0. Test: 0.1. Safe: 10. +elif opts.sc == 0.5: + sc = 0.5 + sc_beta = 1.5 + epsFact_consrv_diffusion = 10.0 # default: 1.0. Test: 0.1. Safe: 10. +else: + import sys + sys.quit() +ns_forceStrongDirichlet = opts.strong_dir +backgroundDiffusionFactor=0.01 +if useMetrics: + ns_shockCapturingFactor = sc + ns_lag_shockCapturing = True + ns_lag_subgridError = True + ls_shockCapturingFactor = sc + ls_lag_shockCapturing = True + ls_sc_uref = 1.0 + ls_sc_beta = sc_beta + vof_shockCapturingFactor = sc + vof_lag_shockCapturing = True + vof_sc_uref = 1.0 + vof_sc_beta = sc_beta + 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 = epsFact_consrv_diffusion + redist_Newton = True#False + kappa_shockCapturingFactor = sc + kappa_lag_shockCapturing = False#True + kappa_sc_uref = 1.0 + kappa_sc_beta = sc_beta + dissipation_shockCapturingFactor = sc + dissipation_lag_shockCapturing = False#True + dissipation_sc_uref = 1.0 + dissipation_sc_beta = sc_beta +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#True + 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 + +tolfac = 0.001 +mesh_tol = 0.001 +ns_nl_atol_res = max(1.0e-8,tolfac*he**2) +vof_nl_atol_res = max(1.0e-8,tolfac*he**2) +ls_nl_atol_res = max(1.0e-8,tolfac*he**2) +mcorr_nl_atol_res = max(1.0e-8,tolfac*he**2) +rd_nl_atol_res = max(1.0e-8,0.1*he) +kappa_nl_atol_res = max(1.0e-8,tolfac*he**2) +dissipation_nl_atol_res = max(1.0e-8,tolfac*he**2) +mesh_nl_atol_res = max(1.0e-8,mesh_tol*he**2) +am_nl_atol_res = 0.001#max(1.0e-8,mesh_tol*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 + +def twpflowPressure_init(x, t): + p_L = 0.0 + phi_L = tank_dim[nd-1] - waterLevel + phi = x[nd-1] - waterLevel + return p_L -g[nd-1]*(rho_0*(phi_L - phi)+(rho_1 -rho_0)*(smoothedHeaviside_integral(epsFact_consrv_heaviside*opts.he,phi_L) + -smoothedHeaviside_integral(epsFact_consrv_heaviside*opts.he,phi))) diff --git a/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillationIB_so.py b/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillationIB_so.py new file mode 100644 index 00000000..47539cfe --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillationIB_so.py @@ -0,0 +1,71 @@ +""" +Split operator module for two-phase flow +""" + +import os +from proteus.default_so import * +from proteus import Context + +# Create context from main module +name_so = os.path.basename(__file__) +if '_so.py' in name_so[-6:]: + name = name_so[:-6] +elif '_so.pyc' in name_so[-7:]: + name = name_so[:-7] +else: + raise(NameError, 'Split operator module must end with "_so.py"') + +case = __import__(name) +Context.setFromModule(case) +ct = Context.get() + +# List of p/n files +pnList = [("vof3p_p", "vof3p_n"),#0 + ("ls3p_p", "ls3p_n"),#1 + ("redist3p_p", "redist3p_n"),#2 + ("ls_consrv3p_p", "ls_consrv3p_n"),#3 + ("rans3p_p", "rans3p_n"),#4 + ("pressureincrement_p", "pressureincrement_n"),#5 + ("pressure_p", "pressure_n"),#6 + ("added_massIB_p","added_massIB_n"),#7 + ("pressureInitial_p", "pressureInitial_n")]#8 +modelSpinUpList = [8] +#systemStepControllerType = ISO_fixed_MinAdaptiveModelStep +#if ct.dt_fixed: +# systemStepControllerType = Sequential_FixedStep +# dt_system_fixed = ct.dt_fixed +# systemStepExact=False +#else: # use CFL +# systemStepControllerType = Sequential_MinAdaptiveModelStep +# systemStepExact=False +#class Sequential_PS(Sequential_FixedStep): +# def __init__(self,modelList,system=defaultSystem,stepExact=True): +# Sequential_FixedStep.__init__(self,modelList,system,stepExact) +# self.modelList = modelList[:len(pnList)] + +class Sequential_FixedStepPS(Sequential_FixedStep): + def __init__(self,modelList,system=defaultSystem,stepExact=True): + Sequential_FixedStep.__init__(self,modelList,system,stepExact) + self.modelList = modelList[:len(pnList)-1] +systemStepControllerType = Sequential_FixedStepPS +systemStepExact=False +dt_system_fixed = ct.opts.dt_fixed + +needEBQ_GLOBAL = False +needEBQ = False + +#modelSpinUpList = [0] # for initial conditions of movemesh +#archiveFlag = ArchiveFlags.EVERY_SEQUENCE_STEP + +if ct.opts.nsave == 0: + if ct.dt_fixed > 0: + archiveFlag = ArchiveFlags.EVERY_USER_STEP + if ct.dt_init < ct.dt_fixed: + tnList = [0., ct.dt_init, ct.dt_fixed, ct.T] + else: + tnList = [0., ct.dt_fixed, ct.T] + else: + tnList = [0., ct.dt_init, ct.T] +else: + tnList=[0.0,ct.dt_init]+[ct.dt_init+ i*ct.dt_out for i in range(1,ct.nDTout+1)] + diff --git a/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillation_so.py b/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillation_so.py new file mode 100644 index 00000000..2b2b3695 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/circle_caisson2D_oscillation_so.py @@ -0,0 +1,74 @@ +""" +Split operator module for two-phase flow +""" + +import os +from proteus.default_so import * +from proteus import Context + +# Create context from main module +name_so = os.path.basename(__file__) +if '_so.py' in name_so[-6:]: + name = name_so[:-6] +elif '_so.pyc' in name_so[-7:]: + name = name_so[:-7] +else: + raise(NameError, 'Split operator module must end with "_so.py"') + +case = __import__(name) +Context.setFromModule(case) +ct = Context.get() + +# List of p/n files +pnList = [] + +# moving mesh +if ct.movingDomain: + pnList += [("moveMesh_p", "moveMesh_n")] + +# added mass +if ct.opts.addedMass: + pnList += [("added_mass_p","added_mass_n")] + +# Navier-Stokes and VOF +pnList += [("twp_navier_stokes_p", "twp_navier_stokes_n"), + ("vof_p", "vof_n")] + +# Level set +if not ct.useOnlyVF: + pnList += [("ls_p", "ls_n"), + ("redist_p", "redist_n"), + ("ls_consrv_p", "ls_consrv_n")] + +# Turbulence +if ct.useRANS > 0: + pnList += [("kappa_p", "kappa_n"), + ("dissipation_p", "dissipation_n")] + +#systemStepControllerType = ISO_fixed_MinAdaptiveModelStep +if ct.dt_fixed: + systemStepControllerType = Sequential_FixedStep + dt_system_fixed = ct.dt_fixed + systemStepExact=False +else: # use CFL + systemStepControllerType = Sequential_MinAdaptiveModelStep + systemStepExact=False + +needEBQ_GLOBAL = False +needEBQ = False + +modelSpinUpList = [0] # for initial conditions of movemesh +archiveFlag = ArchiveFlags.EVERY_SEQUENCE_STEP + +if ct.opts.nsave == 0: + if ct.dt_fixed > 0: + archiveFlag = ArchiveFlags.EVERY_USER_STEP + if ct.dt_init < ct.dt_fixed: + tnList = [0., ct.dt_init, ct.dt_fixed, ct.T] + else: + tnList = [0., ct.dt_fixed, ct.T] + else: + tnList = [0., ct.dt_init, ct.T] +else: + tnList=[0.0,ct.dt_init]+[ct.dt_init+ i*ct.dt_out for i in range(1,ct.nDTout+1)] + diff --git a/2d/floatingStructures/floatingCaissonAddedMass/cylinder2D.py b/2d/floatingStructures/floatingCaissonAddedMass/cylinder2D.py new file mode 100644 index 00000000..f5379e57 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/cylinder2D.py @@ -0,0 +1,535 @@ + +import numpy as np +from math import cos, ceil, pi +from proteus import (Domain, + Context, + WaveTools as wt) +from proteus.mprans import SpatialTools as st +from proteus.Profiling import logEvent +import proteus.TwoPhaseFlow.TwoPhaseFlowProblem as TpFlow +#from proteus.mbd import ChRigidBody as crb #chrono update +from proteus.mbd import CouplingFSI as cfsi +#from proteus.mbd import pyChronoCore as pych +import pychrono as pych +from proteus.Profiling import logEvent + +from proteus import Comm +comm=Comm.init() + +# ____ _ _ ___ _ _ +# / ___|___ _ __ | |_ _____ _| |_ / _ \ _ __ | |_(_) ___ _ __ ___ +# | | / _ \| '_ \| __/ _ \ \/ / __| | | | | '_ \| __| |/ _ \| '_ \/ __| +# | |__| (_) | | | | || __/> <| |_ | |_| | |_) | |_| | (_) | | | \__ \ +# \____\___/|_| |_|\__\___/_/\_\\__| \___/| .__/ \__|_|\___/|_| |_|___/ +# |_| +# Context options +# used in command line directly with option -C +# e.g.: parun [...] -C "g=(0.,-9.81,0.) rho_0=998.2 genMesh=False" +# +# only change/add context options in the "other options section" +# other sections have variables used in _p and _n files + +context_options = [] +# physical constants +context_options += [ +# ("rho_0", 998.2, "Water density"), + ("rho_0", 1037.3, "Water (bottom fluid) density"), + ("nu_0", 1.004e-6, "Water kinematic viscosity m/sec^2"), +# ("rho_1", 1.205, "Air Densiy"), + ("rho_1", 1014.1, "Air (top fluid) Densiy"), + ("nu_1", 1.5e-5, "Air kinematic viscosity m/sec^2"), + ("sigma_01", 0., "Surface tension"), + ("g", (0, -9.81, 0), "Gravitational acceleration vector"), + ] +# run time options +context_options += [ + ("T", 0.5 ,"Simulation time in s"), + ("dt_init", 0.001 ,"Value of initial time step"), + ("dt_fixed", None, "Value of maximum time step"), + ("archiveAllSteps", False, "archive every steps"), + ("dt_output", 0.05, "number of saves per second"), + ("runCFL", 0.5 ,"Target CFL value"), + ("cfl", 0.5 ,"Target CFL value"), + ] + +context_options += [ + ("water_level", 3.5, "Height of free surface above bottom"), + ("ecH", 1.5, "how many elements around phi=0"), + # tank + ('tank_wavelength_scale', True, 'if True: tank_x=value*wavelength, tank_y=value*wavelength'), + ('tank_x', 5., 'Length of tank'), + ('tank_z', 6., 'Height of tank'), + ('tank_sponge_lambda_scale', True, 'True: sponge length relative to wavelength'), + ('tank_sponge_xn', 0., 'length of sponge layer x-'), + ('tank_sponge_xp', 0., 'length of sponge layer x+'), + ('tank_sponge_gen', 'x-', 'sponge layers to be gen zones (if waves)'), + ('tank_sponge_abs', 'x+', 'sponge layers to be abs zones'), + ('IC', 'AtRest', 'Initial Conditions: Perturbed or AtRest'), + # chrono options + ("sampleRate", 0., "sampling rate for chrono. 0 for every timestep"), + # sphere + ("fixStructure", False, "fix structure in place"), + ("free_x", (1., 1., 0.), "Translational DOFs"), + ("free_r", (0., 0., 1.), "Rotational DOFs"), + # waves + ("waves", False, "Generate waves (True/False)"), + ("wave_period", 1.2, "Period of the waves"), + ("wave_height",0.10, "Height of the waves"), + ("wave_dir", (1., 0., 0.), "Direction of the waves (from left boundary)"), + # mesh refinement +# ("he", 0.05, "Set characteristic element size"), + ("he", 0.5, "Set characteristic element size"), + # numerical options + ("genMesh", True, "True: generate new mesh every time. False: do not generate mesh if file exists"), + ("use_gmsh", False, "use_gmsh"), + ("refinement", True, "ref"), + ("refinement_freesurface", 0.05, "ref"), + ("refinement_grading", 1.2, "ref"), + ("movingDomain", True, "True/False"), + ("addedMass", False, "True/False"), + ("chrono_dt", 1e-4, "time step in chrono"), + ("timeIntegration", "backwardEuler", "Time integration scheme (backwardEuler/VBDF)"), + ("useRANS", 0, "RANS model"), + ("useSuperlu", True, "RANS model"), + ] +# instantiate context options +opts=Context.Options(context_options) + + +rho_0 = opts.rho_0 +nu_0 = opts.nu_0 +rho_1 = opts.rho_1 +nu_1 = opts.nu_1 +sigma_01= opts.sigma_01 +#g = np.array(opts.g) +pyg = np.array(opts.g) +he = opts.he + +# ----- CONTEXT ------ # + +# general options +water_level = opts.water_level + +# tank options +tank_dim = np.array([opts.tank_x, opts.tank_z]) +free_x = np.array(opts.free_x) +free_r = np.array(opts.free_r) +chrono_dt = opts.chrono_dt + +wavelength=1. +# general options + +wave_to_initial_conditions = True + +if opts.waves is True: + mwl = depth = opts.water_level + direction = np.array(opts.wave_dir) + height = opts.wave_height + period = opts.wave_period + BCoeffs = np.zeros(3) + YCoeffs = np.zeros(3) + wave = wt.MonochromaticWaves(period=period, + waveHeight=height, + mwl=mwl, + depth=depth, + g=g, + waveDir=direction, + wavelength=wavelength, + waveType='Linear', + Ycoeff=YCoeffs, + Bcoeff=BCoeffs, + Nf=len(BCoeffs), + fast=False) + wavelength = wave.wavelength + + +# ____ _ +# | _ \ ___ _ __ ___ __ _(_)_ __ +# | | | |/ _ \| '_ ` _ \ / _` | | '_ \ +# | |_| | (_) | | | | | | (_| | | | | | +# |____/ \___/|_| |_| |_|\__,_|_|_| |_| +# Domain +# All geometrical options go here (but not mesh options) + +#domain = Domain.PiecewiseLinearComplexDomain() +#domain2 = Domain.PiecewiseLinearComplexDomain() +domain = Domain.PlanarStraightLineGraphDomain() +domain2 = Domain.PlanarStraightLineGraphDomain() + +nd=2 + +# ----- SHAPES ----- # + +# TANK +if opts.tank_wavelength_scale and opts.waves: + tank_dim[0:1] *= wavelength +tank = st.Tank2D(domain, tank_dim) +sponges = {'x-': opts.tank_sponge_xn, + 'x+': opts.tank_sponge_xp} +if opts.tank_sponge_lambda_scale and opts.waves: + for key in sponges: + sponges[key] *= wave.wavelength +tank.setSponge(x_n=sponges['x-'], x_p=sponges['x+']) +# to change individual BC functions, example: +# tank.BC['x-'].u_dirichlet.uOfXT = lambda x, t: 3.*t + + +# SPHERE +eps=2.0; +coords = np.array([tank_dim[0]/2., tank_dim[1]/2.+eps]) +barycenter = np.array([tank_dim[0]/2., tank_dim[1]/2.+eps]) +sphere = st.Circle(domain, + radius=0.5, + coords=coords, + barycenter=barycenter, + nPoints=int(ceil(2.*pi*tank_dim[0]/opts.he))) +sphere.setHoles([sphere.coords]) +# for gmsh: +sphere.holes_ind = np.array([0]) +tank.setChildShape(sphere, 0) + +domain.MeshOptions.he = 0.5 +st.assembleDomain(domain) # must be called after defining shapes +# ____ _ ____ _ _ _ _ +# | __ ) ___ _ _ _ __ __| | __ _ _ __ _ _ / ___|___ _ __ __| (_) |_(_) ___ _ __ ___ +# | _ \ / _ \| | | | '_ \ / _` |/ _` | '__| | | | | / _ \| '_ \ / _` | | __| |/ _ \| '_ \/ __| +# | |_) | (_) | |_| | | | | (_| | (_| | | | |_| | |__| (_) | | | | (_| | | |_| | (_) | | | \__ \ +# |____/ \___/ \__,_|_| |_|\__,_|\__,_|_| \__, |\____\___/|_| |_|\__,_|_|\__|_|\___/|_| |_|___/ +# |___/ +# Boundary Conditions + +for key, bc in sphere.BC.items(): + bc.setNoSlip() + +for key, bc in tank.BC.items(): + # fix the nodes on the wall of tank + # in case of moving domain + bc.setFixedNodes() +tank.BC['y+'].setAtmosphere() +tank.BC['y-'].setFreeSlip() +if opts.waves: + tank.BC['x-'].setUnsteadyTwoPhaseVelocityInlet(wave = wave, vert_axis = 2, smoothing= 3.0*opts.he) +else: + tank.BC['x-'].setFreeSlip() +#tank.BC['x+'].setUnsteadyTwoPhaseVelocityInlet(wave = wave, vert_axis = 2, smoothing= 3.0*opts.he) +tank.BC['x+'].setFreeSlip() +tank.BC['sponge'].setNonMaterial() +dragAlpha = 0.5/nu_0 +gen = opts.tank_sponge_gen +abso = opts.tank_sponge_abs +if opts.waves is True: + dragAlpha = 5*(2*np.pi/period)/nu_0 + smoothing = opts.he*3 + tank.setGenerationZones(x_n=('x-' in gen and sponges['x-'] != 0), + x_p=('x+' in gen and sponges['x+'] != 0), + waves=wave, + smoothing=smoothing, + dragAlpha=dragAlpha) +tank.setAbsorptionZones(x_n=('x-' in abso and sponges['x-'] != 0), + x_p=('x+' in abso and sponges['x+'] != 0), + dragAlpha=dragAlpha) + +# ___ _ _ _ _ ____ _ _ _ _ +# |_ _|_ __ (_) |_(_) __ _| | / ___|___ _ __ __| (_) |_(_) ___ _ __ ___ +# | || '_ \| | __| |/ _` | | | | / _ \| '_ \ / _` | | __| |/ _ \| '_ \/ __| +# | || | | | | |_| | (_| | | | |__| (_) | | | | (_| | | |_| | (_) | | | \__ \ +# |___|_| |_|_|\__|_|\__,_|_| \____\___/|_| |_|\__,_|_|\__|_|\___/|_| |_|___/ +# Initial Conditions + +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.ctransportCoefficients import smoothedHeaviside_integral + +smoothing = opts.ecH * opts.he +if opts.IC == 'Perturbed': + k = 0.18 + nd = domain.nd + def signedDistance(x): + water_level = (wave.eta(x,0)+wave.mwl) + phi_z = x[2]-water_level + return water_level,phi_z + + def eta_IC(x,t): + return wave.eta(x,0) + + def vel_u(x,t): + return wave.u(x,0) + + class P_IC: + #def __init__(self,waterLevel): + # self.waterLevel = opts.waterLine_z + def uOfXT(self,x,t): + if signedDistance(x)[1] < 0: + H = 0 + dyn = -(eta_IC(x,t)*((cosh(k*(water_level-(water_level-x[2]))))/(cosh(water_level*k)))*rho_0*g[2]) + hyd = -(tank_dim[2] - signedDistance(x)[0])*rho_1*g[2]-(water_level - x[2])*rho_0*g[2] + pressure = dyn + hyd + p_air = 0.0 + elif 0 < signedDistance(x)[1] <= smoothing: + H = smoothedHeaviside(smoothing/2. , signedDistance(x)[1] - smoothing / 2.) + x_max = x[:] + x_max[0]=x[0] + x_max[1]=x[1] + x_max[2]=x[2] - signedDistance(x)[1] + dyn = -(eta_IC(x,t)*((cosh(k*(water_level-(-eta_IC(x,t)))))/(cosh(water_level*k)))*rho_0*g[2]) + hyd = -(tank_dim[2] - signedDistance(x_max)[0])*rho_1*g[2]-(water_level - x_max[2])*rho_0*g[2] + pressure = dyn + hyd + p_air = -(tank_dim[2] - signedDistance(x)[0])*rho_1*g[2] + else: + H = 1 + p_air = -(tank_dim[2] - signedDistance(x)[0])*rho_1*g[2] + pressure = 0.0 + p = H * p_air + (1-H) * pressure + return p + class U_IC: + def uOfXT(self,x,t): + return tank.BC['x-'].u_dirichlet.uOfXT(x, t) + class V_IC: + def uOfXT(self,x,t): + return tank.BC['x-'].v_dirichlet.uOfXT(x, t) + class W_IC: + def uOfXT(self,x,t): + return tank.BC['x-'].w_dirichlet.uOfXT(x, t) + class VF_IC: + def uOfXT(self, x, t): + return tank.BC['x-'].vof_dirichlet.uOfXT(x, t) + class PHI_IC: + def uOfXT(self, x, t): + return x[nd-1] - signedDistance(x)[0] + +elif opts.IC == 'AtRest': + class P_IC: + def uOfXT(self, x, t): + p_L = 0.0 + phi_L = tank_dim[nd-1] - water_level + phi = x[nd-1] - water_level + return p_L -pyg[nd-1]*(rho_0*(phi_L - phi) + +(rho_1 -rho_0)*(smoothedHeaviside_integral(smoothing,phi_L) + -smoothedHeaviside_integral(smoothing,phi))) + class U_IC: + def uOfXT(self, x, t): + return 0.0 + class V_IC: + def uOfXT(self, x, t): + return 0.0 + class W_IC: + def uOfXT(self, x, t): + return 0.0 + class VF_IC: + def uOfXT(self, x, t): + return smoothedHeaviside(smoothing,x[nd-1]-water_level) + class PHI_IC: + def uOfXT(self, x, t): + return x[nd-1] - water_level + +# instanciating the classes for *_p.py files +initialConditions = {'pressure': P_IC(), + 'vel_u': U_IC(), + 'vel_v': V_IC(), + 'vel_w': W_IC(), + 'vof': VF_IC(), + 'ncls': PHI_IC(), + 'rdls': PHI_IC()} + +# ____ _ +# / ___| |__ _ __ ___ _ __ ___ +# | | | '_ \| '__/ _ \| '_ \ / _ \ +# | |___| | | | | | (_) | | | | (_) | +# \____|_| |_|_| \___/|_| |_|\___/ +# Chrono + +g = pych.ChVectorD(0.,-9.81,0.0) +system = cfsi.ProtChSystem(sampleRate=opts.sampleRate) +system.ChSystem.Set_G_acc(g) +system.update_substeps = True # update drag and added mass forces on cable during substeps +#system.chrono_processor = 0 +#system.parallel_mode = False +system.setTimeStep(chrono_dt) +system.build_kdtree = True +timestepper = "Euler" +if timestepper == "HHT": + system.setTimestepperType("HHT") + + +body = cfsi.ProtChBody(system) +body.attachShape(sphere) +body.setConstraints(free_x=free_x, free_r=free_r) +#mass = 4*np.pi/3*sphere.radius**3*(opts.rho_0 + 1) +mass = 4*np.pi/3*sphere.radius*1040.0 +Ixx = Iyy = Izz = 2./5.*mass*sphere.radius**2 +body.ChBody.SetMass(mass) +inert = pych.ChVectorD(Ixx, Iyy, Izz) +body.ChBody.SetInertiaXX(inert) +body.setRecordValues(all_values=True) +if opts.fixStructure: + body.ChBody.SetBodyFixed(True) + +# _ _ _ +# | \ | |_ _ _ __ ___ ___ _ __(_) ___ ___ +# | \| | | | | '_ ` _ \ / _ \ '__| |/ __/ __| +# | |\ | |_| | | | | | | __/ | | | (__\__ \ +# |_| \_|\__,_|_| |_| |_|\___|_| |_|\___|___/ +# Numerics + +outputStepping = TpFlow.OutputStepping( + final_time=opts.T, + dt_init=opts.dt_init, + # cfl=opts.cfl, + dt_output=opts.dt_output, + nDTout=None, + dt_fixed=opts.dt_fixed, +) + +myTpFlowProblem = TpFlow.TwoPhaseFlowProblem( + ns_model=None, + ls_model=None, + nd=domain.nd, + cfl=opts.cfl, + outputStepping=outputStepping, + structured=False, + he=he, + nnx=None, + nny=None, + nnz=None, + domain=domain, + initialConditions=initialConditions, + boundaryConditions=None, # set with SpatialTools, + useSuperlu=opts.useSuperlu, +) + +myTpFlowProblem.movingDomain = opts.movingDomain + +params = myTpFlowProblem.Parameters + +# line below needed for relaxation zones +# (!) hack +myTpFlowProblem.Parameters.Models.rans2p.auxiliaryVariables += domain.auxiliaryVariables['twp'] +myTpFlowProblem.archiveAllSteps = True + +# MESH PARAMETERS +params.mesh.genMesh = opts.genMesh +params.mesh.he = opts.he + +# PHYSICAL PARAMETERS +params.physical.densityA = opts.rho_0 # water +params.physical.densityB = opts.rho_1 # air +params.physical.kinematicViscosityA = opts.nu_0 # water +params.physical.kinematicViscosityB = opts.nu_1 # air +params.physical.gravity = np.array(opts.g) +params.physical.surf_tension_coeff = opts.sigma_01 + +# MODEL PARAMETERS +ind = -1 +if opts.movingDomain: + params.Models.moveMeshElastic.index = ind+1 + ind += 1 +params.Models.rans2p.index = ind+1 +ind += 1 +params.Models.vof.index = ind+1 +ind += 1 +params.Models.ncls.index = ind+1 +ind += 1 +params.Models.rdls.index = ind+1 +ind += 1 +params.Models.mcorr.index = ind+1 +ind += 1 +if opts.addedMass is True: + params.Models.addedMass.index = ind+1 + ind += 1 + +# auxiliary variables +params.Models.rans2p.auxiliaryVariables += [system] +#params.Models.rans2p.weak_bc_penalty_constant = 10./nu_0#Re + +if opts.addedMass is True: + # passed in added_mass_p.py coefficients + params.Models.addedMass.auxiliaryVariables += [system.ProtChAddedMass] + max_flag = 0 + max_flag = max(domain.vertexFlags) + max_flag = max(domain.segmentFlags+[max_flag]) + max_flag = max(domain.facetFlags+[max_flag]) + flags_rigidbody = np.zeros(max_flag+1, dtype='int32') + for s in system.subcomponents: + if type(s) is cfsi.ProtChBody: + for i in range(s.i_start, s.i_end): + flags_rigidbody[i] = 1 + params.Models.addedMass.flags_rigidbody = flags_rigidbody + +# __ __ _ ___ _ _ +# | \/ | ___ ___| |__ / _ \ _ __ | |_(_) ___ _ __ ___ +# | |\/| |/ _ \/ __| '_ \ | | | | '_ \| __| |/ _ \| '_ \/ __| +# | | | | __/\__ \ | | | | |_| | |_) | |_| | (_) | | | \__ \ +# |_| |_|\___||___/_| |_| \___/| .__/ \__|_|\___/|_| |_|___/ +# |_| + +he = opts.he + +mesh_fileprefix = 'mesh'+str(int(he*1000)) +domain.MeshOptions.he = he +domain.MeshOptions.setTriangleOptions() +domain.use_gmsh = opts.use_gmsh +domain.MeshOptions.genMesh = opts.genMesh +domain.MeshOptions.use_gmsh = opts.use_gmsh +domain.MeshOptions.setOutputFiles(name=mesh_fileprefix) + +st.assembleDomain(domain) # must be called after defining shapes + +# prescribed_init = False + +# if prescribed_init: +# logEvent('Calculating chrono prescribed motion before starting simulation with dt='+str(time_init_dt)+' for '+str(time_init)+' seconds (this migth take some time)') +# system.calculate_init() +# system.setTimeStep(time_init_dt) +# system.calculate(time_init) # run chrono before fluid sim for intended time to executed prescribed motion +# for i in range(int(time_init/1e-3)): +# system.calculate(1e-3) # run chrono before fluid sim for intended time to executed prescribed motion +# logEvent('finished prescribed motion with body at position '+str(body.ChBody.GetPos())) +# system.setTimeStep(opts.chrono_dt) + +if opts.use_gmsh and opts.refinement is True: + grading = np.cbrt(opts.refinement_grading*12/np.sqrt(2))/np.cbrt(1.*12/np.sqrt(2)) # convert change of volume to change of element size + import py2gmsh + from py2gmsh import geometry2mesh + mesh = geometry2mesh(domain) + grading = np.cbrt(opts.refinement_grading*12/np.sqrt(2))/np.cbrt(1.*12/np.sqrt(2)) # convert change of volume to change of element size + he = opts.he + he_max = 100. + he_max_water = 100. + field_list = [] + def mesh_grading(start, he, grading): + return '(abs({start})*({grading}-1)+{he})/{grading}'.format(he=he, start=start, grading=grading) + def dist_plane(xn, xp, plane='x'): + x_range = abs(xp-xn) + dist = '0.5*(abs({plane}-({xn}))+abs({plane}-({xp}))-{x_range})'.format(xn=xn, xp=xp, x_range=x_range, plane=plane) + return dist + + box = opts.wave_height/2. + me1 = py2gmsh.Field.MathEval(mesh=mesh) + dist_z = dist_plane(xn=water_level-box, xp=water_level+box, plane='z') + #dist = 'sqrt(({dist_x})^2+({dist_y})^2+({dist_z})^2)'.format(dist_x=dist_x, dist_y=dist_y, dist_z=dist_z) + dist = dist_z + me1.F = mesh_grading(start=dist, he=he, grading=grading) + field_list += [me1] + + me3 = py2gmsh.Field.MathEval(mesh=mesh) + dist = 'Sqrt(({x_center}-x)^2+({y_center}-y)^2+({z_center}-z)^2)-radius'.format(x_center=sphere.coords[0], + y_center=sphere.coords[1], + z_center=sphere.coords[2]) + me3.F = mesh_grading(he=he, start=dist, grading=grading) + field_list += [me3] + + # background field + fmin = py2gmsh.Field.Min(mesh=mesh) + fmin.FieldsList = field_list + mesh.setBackgroundField(fmin) + + # max element size + mesh.Options.Mesh.CharacteristicLengthMax = he_max + mesh.Coherence = True + + mesh.writeGeo(mesh_fileprefix+'.geo') + +# system.calculate_init() +# for i in range(100000): +# print(i, system.GetChTime(), body.ChBody.GetPos()) +# system.calculate(0.01) +# print(m1.getTensionBack(), m2.getTensionBack(), m3.getTensionBack()) diff --git a/2d/floatingStructures/floatingCaissonAddedMass/ls3p_n.py b/2d/floatingStructures/floatingCaissonAddedMass/ls3p_n.py new file mode 100644 index 00000000..2cce7b75 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/ls3p_n.py @@ -0,0 +1,82 @@ +from proteus.default_n import * +import ls3p_p as physics +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools) +from proteus.mprans import NCLS +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +# time stepping +runCFL = ct.runCFL +if ct.timeIntegration == "VBDF": + timeIntegration = TimeIntegration.VBDF + timeOrder = 2 +else: + timeIntegration = TimeIntegration.BackwardEuler_cfl +stepController = StepControl.Min_dt_controller + +# mesh options +nLevels = ct.nLevels +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +restrictFineSolutionToAllMeshes = mesh.restrictFineSolutionToAllMeshes +triangleOptions = mesh.triangleOptions + + + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +femSpaces = {0: ct.basis} + +massLumping = False +conservativeFlux = None +numericalFluxType = NCLS.NumericalFlux + +subgridError = NCLS.SubgridError(coefficients=physics.coefficients, + nd=ct.domain.nd) +shockCapturing = NCLS.ShockCapturing(physics.coefficients, + ct.domain.nd, + shockCapturingFactor=ct.ls_shockCapturingFactor, + lag=ct.ls_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = LinearAlgebraTools.SparseMatrix + +if ct.useOldPETSc: + multilevelLinearSolver = LinearSolvers.PETSc + levelLinearSolver = LinearSolvers.PETSc +else: + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py +if ct.useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + +linear_solver_options_prefix = 'ncls_' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.001 +l_atol_res = 0.001*ct.ls_nl_atol_res +nl_atol_res = ct.ls_nl_atol_res +useEisenstatWalker = False#True + +maxNonlinearIts = 50 +maxLineSearches = 0 + +auxiliaryVariables = ct.domain.auxiliaryVariables['ls'] diff --git a/2d/floatingStructures/floatingCaissonAddedMass/ls3p_p.py b/2d/floatingStructures/floatingCaissonAddedMass/ls3p_p.py new file mode 100644 index 00000000..03a29619 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/ls3p_p.py @@ -0,0 +1,37 @@ +from proteus.default_p import * +from proteus.mprans import NCLS +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = domain.nd +mesh = domain.MeshOptions + + +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T + +LevelModelType = NCLS.LevelModel + +coefficients = NCLS.Coefficients(V_model=4, + RD_model=2, + ME_model=1, + checkMass=False, + useMetrics=ct.useMetrics, + epsFact=ct.epsFact_consrv_heaviside, + sc_uref=ct.ls_sc_uref, + sc_beta=ct.ls_sc_beta, + movingDomain=ct.movingDomain) + +dirichletConditions = {0: lambda x, flag: None} + +advectiveFluxBoundaryConditions = {} + +diffusiveFluxBoundaryConditions = {0: {}} + +class PHI_IC: + def uOfXT(self, x, t): + return x[nd-1] - ct.waterLevel + +initialConditions = {0: PHI_IC()} diff --git a/2d/floatingStructures/floatingCaissonAddedMass/ls_consrv3p_n.py b/2d/floatingStructures/floatingCaissonAddedMass/ls_consrv3p_n.py new file mode 100644 index 00000000..13ac0ca8 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/ls_consrv3p_n.py @@ -0,0 +1,71 @@ +from proteus.default_n import * +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools, + NumericalFlux) +import ls_consrv3p_p as physics +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +#time stepping +runCFL = ct.runCFL +timeIntegrator = TimeIntegration.ForwardIntegrator +timeIntegration = TimeIntegration.NoIntegration + +#mesh options +nLevels = ct.nLevels +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +restrictFineSolutionToAllMeshes = mesh.restrictFineSolutionToAllMeshes +triangleOptions = mesh.triangleOptions + + + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +femSpaces = {0: ct.basis} + +subgridError = None +massLumping = False +numericalFluxType = NumericalFlux.DoNothing +conservativeFlux = None +shockCapturing = None + +fullNewtonFlag = True +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = LinearAlgebraTools.SparseMatrix + +if ct.useOldPETSc: + multilevelLinearSolver = LinearSolvers.PETSc + levelLinearSolver = LinearSolvers.PETSc +else: + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py + +if ct.useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + +linear_solver_options_prefix = 'mcorr_' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.001 +l_atol_res = 0.001*ct.mcorr_nl_atol_res +nl_atol_res = ct.mcorr_nl_atol_res +useEisenstatWalker = False#True + +maxNonlinearIts = 50 +maxLineSearches = 0 diff --git a/2d/floatingStructures/floatingCaissonAddedMass/ls_consrv3p_p.py b/2d/floatingStructures/floatingCaissonAddedMass/ls_consrv3p_p.py new file mode 100644 index 00000000..d2fd8ba9 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/ls_consrv3p_p.py @@ -0,0 +1,37 @@ +from proteus.default_p import * +from proteus.mprans import MCorr +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + + +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T + +LevelModelType = MCorr.LevelModel + +coefficients = MCorr.Coefficients(LSModel_index=1, + V_model=4, + me_model=3, + VOFModel_index=0, + applyCorrection=ct.applyCorrection, + nd=nd, + checkMass=True, + useMetrics=ct.useMetrics, + epsFactHeaviside=ct.epsFact_consrv_heaviside, + epsFactDirac=ct.epsFact_consrv_dirac, + epsFactDiffusion=ct.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/floatingStructures/floatingCaissonAddedMass/pressureInitial_n.py b/2d/floatingStructures/floatingCaissonAddedMass/pressureInitial_n.py new file mode 100644 index 00000000..d860f0ec --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/pressureInitial_n.py @@ -0,0 +1,50 @@ +from __future__ import absolute_import +from proteus import * +from proteus.default_n import * +try: + from .pressureInitial_p import * +except: + from pressureInitial_p import * + +timeIntegration=NoIntegration + +triangleOptions = triangleOptions + +femSpaces = {0:pbasis} + +stepController=FixedStep + +#numericalFluxType = NumericalFlux.ConstantAdvection_Diffusion_SIPG_exterior #weak boundary conditions (upwind ?) +matrix = LinearAlgebraTools.SparseMatrix + +if useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + parallelPartitioningType = parallelPartitioningType + nLayersOfOverlapForParallel = nLayersOfOverlapForParallel + nonlinearSmoother = None + linearSmoother = None + +numericalFluxType = NumericalFlux.ConstantAdvection_exterior + +linear_solver_options_prefix = 'pinit_' + +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +#linear solve rtolerance + +linTolFac = 0.0 +l_atol_res = 1.0e-10 +tolFac = 0.0 +nl_atol_res = 1.0e-10 +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' +maxLineSearches=0 +periodicDirichletConditions=None + +conservativeFlux=None diff --git a/2d/floatingStructures/floatingCaissonAddedMass/pressureInitial_p.py b/2d/floatingStructures/floatingCaissonAddedMass/pressureInitial_p.py new file mode 100644 index 00000000..3aad7078 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/pressureInitial_p.py @@ -0,0 +1,48 @@ +from __future__ import absolute_import +from builtins import object +from math import * +from proteus import * +from proteus.default_p import * +try: + from .circle_caisson2D_oscillationIB import * +except: + from circle_caisson2D_oscillationIB import * + +from proteus.mprans import PresInit + +#domain = ctx.domain +#nd = ctx.nd +name = "pressureInitial" + +coefficients=PresInit.Coefficients(nd=nd, + modelIndex=8, + fluidModelIndex=4, + pressureModelIndex=6) + +#LevelModelType = PresInit.LevelModel + +#pressure increment should be zero on any pressure dirichlet boundaries +def getDBC_pInit(x,flag): + if flag == boundaryTags['top']: + return lambda x,t: 0.0 + +#the advectiveFlux should be zero on any no-flow boundaries +def getAdvectiveFlux_pInit(x,flag): + if flag != boundaryTags['top']: + return lambda x,t: 0.0 + +def getDiffusiveFlux_pInit(x,flag): + if flag != boundaryTags['top']: + return lambda x,t: 0.0 + +class getIBC_pInit(object): + def __init__(self): + pass + def uOfXT(self,x,t): + return 0.0 + +initialConditions = {0:getIBC_pInit()} + +dirichletConditions = {0:getDBC_pInit } +advectiveFluxBoundaryConditions = {0:getAdvectiveFlux_pInit} +diffusiveFluxBoundaryConditions = {0:{0:getDiffusiveFlux_pInit}} diff --git a/2d/floatingStructures/floatingCaissonAddedMass/pressure_n.py b/2d/floatingStructures/floatingCaissonAddedMass/pressure_n.py new file mode 100644 index 00000000..e9c64d27 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/pressure_n.py @@ -0,0 +1,52 @@ +from proteus import * +from proteus.default_n import * +from pressure_p import * + +timeIntegration = TimeIntegration.NoIntegration + +triangleOptions = triangleOptions + + +femSpaces = {0:pbasis} + + +# stepController = StepControl.Min_dt_cfl_controller +# runCFL= 0.99 +# runCFL= 0.5 + +stepController=FixedStep + +#matrix type +numericalFluxType = NumericalFlux.ConstantAdvection_exterior +matrix = LinearAlgebraTools.SparseMatrix + +linear_solver_options_prefix = 'pressure_' + +if useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU +else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + parallelPartitioningType = parallelPartitioningType + nLayersOfOverlapForParallel = nLayersOfOverlapForParallel + nonlinearSmoother = None + linearSmoother = None + +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +#linear solve rtolerance + +linTolFac = 0.0 +l_atol_res = 1.0e-10 +tolFac = 0.0 +nl_atol_res = 1.0e-10 +maxLineSearches = 0 +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +periodicDirichletConditions=None + +conservativeFlux=None diff --git a/2d/floatingStructures/floatingCaissonAddedMass/pressure_p.py b/2d/floatingStructures/floatingCaissonAddedMass/pressure_p.py new file mode 100644 index 00000000..0775ac4a --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/pressure_p.py @@ -0,0 +1,29 @@ +from math import * +from proteus import * +from proteus.default_p import * +from circle_caisson2D_oscillationIB import * +from proteus.mprans import Pres + +name = "pressure" + +coefficients=Pres.Coefficients(modelIndex=6, + fluidModelIndex=4, + pressureIncrementModelIndex=5, + useRotationalForm=False) + +LevelModelType = Pres.LevelModel + +def getDBC_p(x,flag): + return None + +def getFlux(x,flag): + return lambda x,t: 0.0 + +class getIBC_p: + def uOfXT(self,x,t): + return -(L[1] - x[1])*rho_1*g[1] + +initialConditions = {0:getIBC_p()} + +dirichletConditions = {0:getDBC_p } # pressure bc are explicitly set +advectiveFluxBoundaryConditions = {0:getFlux} diff --git a/2d/floatingStructures/floatingCaissonAddedMass/pressureincrement_n.py b/2d/floatingStructures/floatingCaissonAddedMass/pressureincrement_n.py new file mode 100644 index 00000000..840c9730 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/pressureincrement_n.py @@ -0,0 +1,55 @@ +from proteus import * +from proteus.default_n import * +from pressureincrement_p import * + +timeIntegration = TimeIntegration.NoIntegration + +triangleOptions = triangleOptions + +femSpaces = {0:pbasis} + +stepController=FixedStep + +numericalFluxType = PresInc.NumericalFlux +matrix = LinearAlgebraTools.SparseMatrix + +openTop=True +if openTop: + if useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + else: + multilevelLinearSolver = KSP_petsc4py + levelLinearSolver = KSP_petsc4py + parallelPartitioningType = parallelPartitioningType + nLayersOfOverlapForParallel = nLayersOfOverlapForParallel + nonlinearSmoother = None + linearSmoother = None +else: + linearSmoother = LinearSolvers.NavierStokesPressureCorrection # pure neumann laplacian solver + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py + +linearSmoother = LinearSolvers.NavierStokesPressureCorrection # pure neumann laplacian solver +multilevelLinearSolver = LinearSolvers.KSP_petsc4py +levelLinearSolver = LinearSolvers.KSP_petsc4py + +linear_solver_options_prefix = 'phi_' + +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +#linear solve rtolerance + +linTolFac = 0.0 +l_atol_res = 1.0e-10 +tolFac = 0.0 +nl_atol_res = 1.0e-10 +nonlinearSolverConvergenceTest = 'r' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' +maxLineSearches=0 +periodicDirichletConditions=None + +#conservativeFlux = {0:'point-eval'} #'point-eval','pwl-bdm-opt' +conservativeFlux=None diff --git a/2d/floatingStructures/floatingCaissonAddedMass/pressureincrement_p.py b/2d/floatingStructures/floatingCaissonAddedMass/pressureincrement_p.py new file mode 100644 index 00000000..583c4426 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/pressureincrement_p.py @@ -0,0 +1,42 @@ +from math import * +from proteus import * +from proteus.default_p import * +from circle_caisson2D_oscillationIB import * + + +#domain = ctx.domain +#nd = ctx.nd +name = "pressureincrement" + +from proteus.mprans import PresInc +rho_s=rho_1 +coefficients=PresInc.Coefficients(rho_f_min = (1.0-1.0e-8)*rho_1, + rho_s_min = (1.0-1.0e-8)*rho_s, + nd = nd, + modelIndex=5, + fluidModelIndex=4) + +LevelModelType = PresInc.LevelModel + +#pressure increment should be zero on any pressure dirichlet boundaries +def getDBC_phi(x,flag): + return None + +#the advectiveFlux should be zero on any no-flow boundaries +def getAdvectiveFlux_qt(x,flag): + if onBoundary(x): + return lambda x,t: 0.0 + +def getDiffusiveFlux_phi(x,flag): + if onBoundary(x): + return lambda x,t: 0.0 + +class getIBC_phi: + def uOfXT(self,x,t): + return 0.0 + +initialConditions = {0:getIBC_phi()} + +dirichletConditions = {0:getDBC_phi} +advectiveFluxBoundaryConditions = {0:getAdvectiveFlux_qt} +diffusiveFluxBoundaryConditions = {0:{0:getDiffusiveFlux_phi}} diff --git a/2d/floatingStructures/floatingCaissonAddedMass/rans3p_n.py b/2d/floatingStructures/floatingCaissonAddedMass/rans3p_n.py new file mode 100644 index 00000000..360041be --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/rans3p_n.py @@ -0,0 +1,69 @@ +from proteus import * +from proteus.default_n import * +from rans3p_p import * +from circle_caisson2D_oscillationIB import * +timeDiscretization='be' +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 = {0:time_tol,1:time_tol} + rtol_u = {0:time_tol,1:time_tol} +else: + timeIntegration = BackwardEuler_cfl + stepController = Min_dt_cfl_controller + +femSpaces = {0:vbasis, + 1:vbasis} +#femSpaces = {0:C0_AffineQuadraticOnSimplexWithNodalBasis, +# 1:C0_AffineQuadraticOnSimplexWithNodalBasis} + +massLumping = False +numericalFluxType = None +conservativeFlux = None + +numericalFluxType = RANS3PF.NumericalFlux +subgridError = RANS3PF.SubgridError(coefficients,nd,lag=ns_lag_subgridError,hFactor=hFactor) +shockCapturing = RANS3PF.ShockCapturing(coefficients,nd,ns_shockCapturingFactor,lag=ns_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton + +nonlinearSmoother = None + +linearSmoother = None#SimpleNavierStokes2D + +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 = 'rans3p_' +nonlinearSolverConvergenceTest = 'rits' +levelNonlinearSolverConvergenceTest = 'rits' +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:'point-eval'} +#auxiliaryVariables=[pointGauges,lineGauges] +auxiliaryVariables=[body] diff --git a/2d/floatingStructures/floatingCaissonAddedMass/rans3p_p.py b/2d/floatingStructures/floatingCaissonAddedMass/rans3p_p.py new file mode 100644 index 00000000..d8a8137f --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/rans3p_p.py @@ -0,0 +1,105 @@ +from proteus import * +from proteus.default_p import * +from circle_caisson2D_oscillationIB import * +from proteus.mprans import RANS3PF + +LevelModelType = RANS3PF.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 + + +mylist_of_sdfs=[]; + +coefficients = RANS3PF.Coefficients(epsFact=epsFact_viscosity, + sigma=0.0, + USE_SUPG = True, + ARTIFICIAL_VISCOSITY=2, + rho_0 = rho_0, + nu_0 = nu_0, + rho_1 = rho_1, + nu_1 = nu_1, + g=g, + nd=nd, + ME_model=4, + PRESSURE_model=6, + SED_model=None, + VOF_model=0, + VOS_model=None, + LS_model=1, + Closure_0_model=None, + Closure_1_model=None, + 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, + dragAlpha=0.0, + PSTAB=0.0, + nParticles=1, + particle_epsFact=0.33, + particle_alpha=0, + particle_beta=0, + particle_penalty_constant=100.0, + particle_sdfList=[], + particle_velocityList=[], + use_sbm=0, + use_ball_as_particle = 1, + ball_center=ball_center, + ball_radius=ball_radius, + ball_velocity=ball_velocity, + ball_angular_velocity=ball_angular_velocity) +#=============================================================================== +# BC +#=============================================================================== + +def getDBC_u(x,flag): + if onBoundary(x): + return lambda x,t: 0.0 + +def getDBC_v(x,flag): + if onBoundary(x): + return lambda x,t: 0.0 + +dirichletConditions = {0:getDBC_u, + 1:getDBC_v} + +def getAFBC_u(x,flag): + return None + +def getAFBC_v(x,flag): + return None + +def getDFBC_u(x,flag): + return None + +def getDFBC_v(x,flag): + return None + +advectiveFluxBoundaryConditions = {0:getAFBC_u, + 1:getAFBC_v} + +diffusiveFluxBoundaryConditions = {0:{0:getDFBC_u}, + 1:{1:getDFBC_v}} + +class AtRest: + def uOfXT(self,x,t): + return 0.0 + +initialConditions = {0:AtRest(), + 1:AtRest()} diff --git a/2d/floatingStructures/floatingCaissonAddedMass/redist3p_n.py b/2d/floatingStructures/floatingCaissonAddedMass/redist3p_n.py new file mode 100644 index 00000000..d43a68b7 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/redist3p_n.py @@ -0,0 +1,96 @@ +from proteus.default_n import * +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools, + NumericalFlux) +from proteus.mprans import RDLS +import redist3p_p as physics +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +# time stepping +runCFL = ct.runCFL + +# mesh options +nLevels = ct.nLevels +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +restrictFineSolutionToAllMeshes = mesh.restrictFineSolutionToAllMeshes +triangleOptions = mesh.triangleOptions + + + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +femSpaces = {0: ct.basis} +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +massLumping = False +numericalFluxType = NumericalFlux.DoNothing +conservativeFlux = None +subgridError = RDLS.SubgridError(coefficients=physics.coefficients, + nd=ct.domain.nd) +shockCapturing = RDLS.ShockCapturing(coefficients=physics.coefficients, + nd=ct.domain.nd, + shockCapturingFactor=ct.rd_shockCapturingFactor, + lag=ct.rd_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +nonlinearSmoother = NonlinearSolvers.NLGaussSeidel +linearSmoother = None + +matrix = LinearAlgebraTools.SparseMatrix + +if ct.useOldPETSc: + multilevelLinearSolver = LinearSolvers.PETSc + levelLinearSolver = LinearSolvers.PETSc +else: + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py +if ct.useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + +if ct.redist_Newton: + timeIntegration = TimeIntegration.NoIntegration + stepController = StepControl.Newton_controller + maxNonlinearIts = 25 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'r' + levelNonlinearSolverConvergenceTest = 'r' + linearSolverConvergenceTest = 'r-true' +else: + timeIntegration = TimeIntegration.BackwardEuler_cfl + stepController = RDLS.PsiTC + runCFL = 0.5 + psitc['nStepsForce'] = 6 + psitc['nStepsMax'] = 25 + psitc['reduceRatio'] = 3.0 + psitc['startRatio'] = 1.0 + rtol_res[0] = 0.0 + atol_res[0] = ct.rd_nl_atol_res + useEisenstatWalker = False#True + maxNonlinearIts = 1 + maxLineSearches = 0 + nonlinearSolverConvergenceTest = 'rits' + levelNonlinearSolverConvergenceTest = 'rits' + linearSolverConvergenceTest = 'r-true' + +linear_solver_options_prefix = 'rdls_' + +nl_atol_res = ct.rd_nl_atol_res +tolFac = 0.0 +linTolFac = 0.001 +l_atol_res = 0.001*ct.rd_nl_atol_res +useEisenstatWalker = False#True diff --git a/2d/floatingStructures/floatingCaissonAddedMass/redist3p_p.py b/2d/floatingStructures/floatingCaissonAddedMass/redist3p_p.py new file mode 100644 index 00000000..5fe7ebb9 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/redist3p_p.py @@ -0,0 +1,41 @@ +from proteus.default_p import * +from proteus.mprans import RDLS +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = domain.nd +mesh = domain.MeshOptions + + +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T + +""" +The redistancing equation in the sloshbox test problem. +""" + +LevelModelType = RDLS.LevelModel + +coefficients = RDLS.Coefficients(applyRedistancing=ct.applyRedistancing, + epsFact=ct.epsFact_redistance, + nModelId=1, + rdModelId=2, + useMetrics=ct.useMetrics, + backgroundDiffusionFactor=ct.backgroundDiffusionFactor) + +def getDBC_rd(x, flag): + pass + +dirichletConditions = {0: getDBC_rd} +weakDirichletConditions = {0: RDLS.setZeroLSweakDirichletBCsSimple} + +advectiveFluxBoundaryConditions = {} +diffusiveFluxBoundaryConditions = {0: {}} + +class PHI_IC: + def uOfXT(self, x, t): + return x[nd-1] - ct.waterLevel + +initialConditions = {0: PHI_IC()} diff --git a/2d/floatingStructures/floatingCaissonAddedMass/redist_n.py b/2d/floatingStructures/floatingCaissonAddedMass/redist_n.py index b4bd6070..ed68bf6e 100644 --- a/2d/floatingStructures/floatingCaissonAddedMass/redist_n.py +++ b/2d/floatingStructures/floatingCaissonAddedMass/redist_n.py @@ -68,7 +68,7 @@ maxNonlinearIts = 25 maxLineSearches = 0 nonlinearSolverConvergenceTest = 'rits' - levelNonlinearSolverConvergenceTest = 'rits' + levelNonlinearSolverConvergenceTest = 'r' linearSolverConvergenceTest = 'r-true' else: timeIntegration = TimeIntegration.BackwardEuler_cfl diff --git a/2d/floatingStructures/floatingCaissonAddedMass/rt.sh b/2d/floatingStructures/floatingCaissonAddedMass/rt.sh new file mode 100755 index 00000000..7379a7bb --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/rt.sh @@ -0,0 +1,8 @@ +mpiexec -np 4 /home/cekees/proteus/scripts/parun caisson2D_oscillation_so.py -C "use_chrono=True cfl=0.33 strong_dir=False caisson_dim=[0.5,0.5] caisson_coords=[1.5,1.49] caisson_mass=125.0 he=0.06 T=5.0 addedMass=True sc=0.5 nc_form=True" -O ~/irb/bridge/petsc.options.superlu_dist -l 5 -v +mpiexec -np 4 /home/cekees/proteus/scripts/parun caisson2D_oscillation_so.py -C "use_chrono=True cfl=0.33 strong_dir=False caisson_dim=[0.5,0.5] caisson_coords=[1.5,1.49] caisson_mass=125.0 he=0.06 T=5.0 addedMass=True sc=0.5 nc_form=False" -O ~/irb/bridge/petsc.options.superlu_dist -l 5 -v +mpiexec -np 4 /home/cekees/proteus/scripts/parun caisson2D_oscillation_so.py -C "use_chrono=True cfl=0.33 strong_dir=False caisson_dim=[0.5,0.5] caisson_coords=[1.5,1.49] caisson_mass=125.0 he=0.06 T=5.0 addedMass=False sc=0.5 nc_form=True" -O ~/irb/bridge/petsc.options.superlu_dist -l 5 -v +mpiexec -np 4 /home/cekees/proteus/scripts/parun caisson2D_oscillation_so.py -C "use_chrono=True cfl=0.33 strong_dir=False caisson_dim=[0.5,0.5] caisson_coords=[1.5,1.49] caisson_mass=125.0 he=0.06 T=5.0 addedMass=False sc=0.5 nc_form=False" -O ~/irb/bridge/petsc.options.superlu_dist -l 5 -v +mpiexec -np 4 /home/cekees/proteus/scripts/parun caisson2D_oscillation_so.py -C "use_chrono=False cfl=0.33 strong_dir=False caisson_dim=[0.5,0.5] caisson_coords=[1.5,1.49] caisson_mass=125.0 he=0.06 T=5.0 addedMass=True sc=0.5 nc_form=True" -O ~/irb/bridge/petsc.options.superlu_dist -l 5 -v +mpiexec -np 4 /home/cekees/proteus/scripts/parun caisson2D_oscillation_so.py -C "use_chrono=False cfl=0.33 strong_dir=False caisson_dim=[0.5,0.5] caisson_coords=[1.5,1.49] caisson_mass=125.0 he=0.06 T=5.0 addedMass=True sc=0.5 nc_form=False" -O ~/irb/bridge/petsc.options.superlu_dist -l 5 -v +mpiexec -np 4 /home/cekees/proteus/scripts/parun caisson2D_oscillation_so.py -C "use_chrono=False cfl=0.33 strong_dir=False caisson_dim=[0.5,0.5] caisson_coords=[1.5,1.49] caisson_mass=125.0 he=0.06 T=5.0 addedMass=False sc=0.5 nc_form=True" -O ~/irb/bridge/petsc.options.superlu_dist -l 5 -v +mpiexec -np 4 /home/cekees/proteus/scripts/parun caisson2D_oscillation_so.py -C "use_chrono=False cfl=0.33 strong_dir=False caisson_dim=[0.5,0.5] caisson_coords=[1.5,1.49] caisson_mass=125.0 he=0.06 T=5.0 addedMass=False sc=0.5 nc_form=False" -O ~/irb/bridge/petsc.options.superlu_dist -l 5 -v diff --git a/2d/floatingStructures/floatingCaissonAddedMass/twp_navier_stokes_n.py b/2d/floatingStructures/floatingCaissonAddedMass/twp_navier_stokes_n.py index 554826fd..a6e36655 100644 --- a/2d/floatingStructures/floatingCaissonAddedMass/twp_navier_stokes_n.py +++ b/2d/floatingStructures/floatingCaissonAddedMass/twp_navier_stokes_n.py @@ -83,7 +83,7 @@ useEisenstatWalker = False#True maxNonlinearIts = 50 maxLineSearches = 0 -conservativeFlux = {0:'pwl-bdm-opt'} +#conservativeFlux = {0:'pwl-bdm-opt'} if ct.opts.caisson is True: auxs = [ct.system] else: diff --git a/2d/floatingStructures/floatingCaissonAddedMass/twp_navier_stokes_p.py b/2d/floatingStructures/floatingCaissonAddedMass/twp_navier_stokes_p.py index 146d1a3c..f2d41dea 100644 --- a/2d/floatingStructures/floatingCaissonAddedMass/twp_navier_stokes_p.py +++ b/2d/floatingStructures/floatingCaissonAddedMass/twp_navier_stokes_p.py @@ -69,7 +69,8 @@ dragAlphaTypes=dragAlphaTypes, dragBetaTypes=dragBetaTypes, epsFact_solid=epsFact_solid, - barycenters=ct.domain.barycenters) + barycenters=ct.domain.barycenters, + NONCONSERVATIVE_FORM=float(ct.opts.nc_form)) dirichletConditions = {0: lambda x, flag: domain.bc[flag].p_dirichlet.init_cython(), diff --git a/2d/floatingStructures/floatingCaissonAddedMass/vof3p_n.py b/2d/floatingStructures/floatingCaissonAddedMass/vof3p_n.py new file mode 100644 index 00000000..7d3f0a0c --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/vof3p_n.py @@ -0,0 +1,82 @@ +from proteus.default_n import * +from proteus import (StepControl, + TimeIntegration, + NonlinearSolvers, + LinearSolvers, + LinearAlgebraTools) +import vof3p_p as physics +from proteus.mprans import VOF +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = ct.domain.nd +mesh = domain.MeshOptions + +# time stepping +runCFL = ct.runCFL +if ct.timeIntegration == "VBDF": + timeIntegration = TimeIntegration.VBDF + timeOrder = 2 +else: + timeIntegration = TimeIntegration.BackwardEuler_cfl +stepController = StepControl.Min_dt_controller + +# mesh options +nLevels = ct.nLevels +parallelPartitioningType = mesh.parallelPartitioningType +nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel +restrictFineSolutionToAllMeshes = mesh.restrictFineSolutionToAllMeshes +triangleOptions = mesh.triangleOptions + + + +elementQuadrature = ct.elementQuadrature +elementBoundaryQuadrature = ct.elementBoundaryQuadrature + +femSpaces = {0: ct.basis} + +massLumping = False +numericalFluxType = VOF.NumericalFlux +conservativeFlux = None +subgridError = VOF.SubgridError(coefficients=physics.coefficients, + nd=ct.domain.nd) +shockCapturing = VOF.ShockCapturing(physics.coefficients, + ct.domain.nd, + shockCapturingFactor=ct.vof_shockCapturingFactor, + lag=ct.vof_lag_shockCapturing) + +fullNewtonFlag = True +multilevelNonlinearSolver = NonlinearSolvers.Newton +levelNonlinearSolver = NonlinearSolvers.Newton + +nonlinearSmoother = None +linearSmoother = None + +matrix = LinearAlgebraTools.SparseMatrix + +if ct.useOldPETSc: + multilevelLinearSolver = LinearSolvers.PETSc + levelLinearSolver = LinearSolvers.PETSc +else: + multilevelLinearSolver = LinearSolvers.KSP_petsc4py + levelLinearSolver = LinearSolvers.KSP_petsc4py + +if ct.useSuperlu: + multilevelLinearSolver = LinearSolvers.LU + levelLinearSolver = LinearSolvers.LU + +linear_solver_options_prefix = 'vof_' +levelNonlinearSolverConvergenceTest = 'r' +linearSolverConvergenceTest = 'r-true' + +tolFac = 0.0 +linTolFac = 0.001 +l_atol_res = 0.001*ct.vof_nl_atol_res +nl_atol_res = ct.vof_nl_atol_res +useEisenstatWalker = False#True + +maxNonlinearIts = 50 +maxLineSearches = 0 + +auxiliaryVariables = ct.domain.auxiliaryVariables['vof'] diff --git a/2d/floatingStructures/floatingCaissonAddedMass/vof3p_p.py b/2d/floatingStructures/floatingCaissonAddedMass/vof3p_p.py new file mode 100644 index 00000000..c9078014 --- /dev/null +++ b/2d/floatingStructures/floatingCaissonAddedMass/vof3p_p.py @@ -0,0 +1,45 @@ +from proteus.default_p import * +from proteus.ctransportCoefficients import smoothedHeaviside +from proteus.mprans import VOF3P +from proteus import Context + +ct = Context.get() +domain = ct.domain +nd = domain.nd +mesh = domain.MeshOptions + + +genMesh = mesh.genMesh +movingDomain = ct.movingDomain +T = ct.T + +LevelModelType = VOF3P.LevelModel +if ct.useOnlyVF: + RD_model = None + LS_model = None +else: + RD_model = 3 + LS_model = 2 + +coefficients = VOF3P.Coefficients(LS_model=1, + V_model=4, + RD_model=2, + ME_model=0, + checkMass=True, + useMetrics=ct.useMetrics, + epsFact=ct.epsFact_vof, + sc_uref=ct.vof_sc_uref, + sc_beta=ct.vof_sc_beta, + movingDomain=ct.movingDomain) + +dirichletConditions = {0: lambda x, flag: domain.bc[flag].vof_dirichlet.init_cython()} + +advectiveFluxBoundaryConditions = {0: lambda x, flag: domain.bc[flag].vof_advective.init_cython()} + +diffusiveFluxBoundaryConditions = {0: {}} + +class VF_IC: + def uOfXT(self, x, t): + return smoothedHeaviside(ct.epsFact_consrv_heaviside*ct.he,x[nd-1]-ct.waterLevel) + +initialConditions = {0: VF_IC()} 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/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: