Add Fluidity cornerflow 2d example

examples/fluidity/README.md

@@ -16,3 +16,4 @@ This table lists the corresponding PyDRex tests:
 | Fluidity example | PyDRex test |
 | --- | --- |
 | `advection2d` | `test_simple_shear_2d.TestOlivineA.test_dudz_pathline` |
+| `corner2d` | `test_corner_flow_2d.TestOlivineA.test_steady4` |

examples/fluidity/advection2d/advection2d.flml

@@ -98,7 +98,7 @@
 
     # Construct 3D velocity gradient matrix and position vector.
     position = np.asarray(X) if len(X) == 3 else np.insert(X, 1, 0)
-    velocity_gradient_init = _velocity_grad(position)
+    velocity_gradient_init = _velocity_grad(t, position)
 
     # Create new mineral with random initial orientations and undeformed initial state.
     # The extra 22 values are stored after the actual CPO data and they are:
@@ -153,12 +153,8 @@
     velocity_gradient = fields["VelocityGradient"]
     if velocity_gradient.shape == (2, 2):
         velocity_gradient = utils.add_dim(velocity_gradient, 1)
-    _log.info("deformation gradient: %s", deformation_gradient_prev.flatten())
-    _log.info("position: %s", cb_interp_position(position_prev, position, t, dt)(t))
-    _log.info("velocity gradient: %s",
-        cb_interp_velocity_grad(velocity_gradient_prev, velocity_gradient, t, dt)(t).flatten()
-    )
 
+    _log.CONSOLE_LOGGER.setLevel("DEBUG")
     deformation_gradient = mineral.update_orientations(
         PARAMS,
         deformation_gradient_prev,
@@ -331,7 +327,7 @@
             <string_value type="code" language="python" lines="20">def val(X, t):
     from advection2d import velocity
 
-    return velocity(X, t)</string_value>
+    return velocity(t, X)</string_value>
           </python>
         </value>
         <output/>
@@ -393,7 +389,7 @@
               <string_value type="code" language="python" lines="20">def val(X, t):
     from advection2d import velocity_gradient
 
-    return velocity_gradient(X, t)</string_value>
+    return velocity_gradient(t, X)</string_value>
             </python>
           </anisotropic_asymmetric>
         </value>

examples/fluidity/advection2d/advection2d.py

@@ -33,7 +33,7 @@ def cb_interp_velocity_grad(velocity_grad_prev, velocity_grad, t, dt):
     """Return ∇v interpolator for CPO integration between `velocity_grad_prev` and `velocity_grad`."""
 
     @nb.njit(fastmath=True)
-    def _interp_vel_grad(int_time):
+    def _interp_vel_grad(int_time, int_position):
         return velocity_grad_prev + (int_time - t + dt) / dt * (
             velocity_grad - velocity_grad_prev
         )
@@ -58,9 +58,9 @@ def _interp_vel_grad(int_time):
 _velocity, _velocity_grad = simple_shear_2d("X", "Z", STRAIN_RATE)
 
 
-def velocity(X, t):
-    return remove_dim(_velocity(add_dim(X, 1)), 1)
+def velocity(t, X):
+    return remove_dim(_velocity(t, add_dim(X, 1)), 1)
 
 
-def velocity_gradient(X, t):
-    return remove_dim(_velocity_grad(add_dim(X, 1)), 1)
+def velocity_gradient(t, X):
+    return remove_dim(_velocity_grad(t, add_dim(X, 1)), 1)

examples/fluidity/corner2d/Makefile

@@ -0,0 +1,66 @@
+# Make options to keep what's left of our sanity.
+.SHELLFLAGS += -u
+.ONESHELL:
+MAKEFLAGS += --warn-undefined-variables
+MAKEFLAGS += --no-builtin-rules
+# Path options, SIM_NAME is the name of the .msh and .py files.
+SIM_NAME := corner2d
+SRC_DIR := corner2d
+OUT_DIR := _out
+# Initial geometry options (in units of metres), passed to `pydrex-mesh`.
+WIDTH := 1e6
+DEPTH := 2e5
+HALF_WIDTH := $(shell python3 -c 'print($(WIDTH)/2)')
+HALF_DEPTH := $(shell python3 -c 'print($(DEPTH)/2)')
+RESOLUTION_HI := 1e2
+RESOLUTION_LO := 1e4
+# Initial conditions, used to parametrise velocity fields and spawn particles.
+# INIT_HORIZ can be a space-delimited array of values (for multiple particles).
+# NOTE: PLATE_SPEED is in cm/yr here.
+PLATE_SPEED := 2.0
+INIT_HORIZ := 3.13e4 9.74e4 2.02e5 3.97e5
+INIT_DEPTH := 2e5
+# Recrystallisation parameters.
+STRESS_EXPONENT := 1.5
+DEFORMATION_EXPONENT := 3.5
+GBM_MOBILITY := 10
+GBS_THRESHOLD := 0.3
+NUCLEATION_EFFICIENCY := 5
+
+$(OUT_DIR)/fluidity.log-0: $(OUT_DIR)/$(SIM_NAME).flml \
+	$(OUT_DIR)/$(SIM_NAME).msh $(OUT_DIR)/$(SIM_NAME).py $(OUT_DIR)/$(SIM_NAME).ini
+	@echo "********** Running fluidity with verbose logging enabled..."
+	./envcheck.sh -f
+	cd $(OUT_DIR) && fluidity -v2 -l $(SIM_NAME).flml
+
+$(OUT_DIR)/$(SIM_NAME).ini: $(OUT_DIR)/$(SIM_NAME).py
+	@echo "********** Setting up initial conditions and recryst. parameters..."
+	echo "[initial conditions]" > $@
+	echo "PLATE_SPEED = $(PLATE_SPEED)" >> $@
+	echo "INIT_HORIZ = $(INIT_HORIZ)" >> $@
+	echo "INIT_DEPTH = $(INIT_DEPTH)" >> $@
+	echo "STRESS_EXPONENT = $(STRESS_EXPONENT)" >> $@
+	echo "DEFORMATION_EXPONENT = $(DEFORMATION_EXPONENT)" >> $@
+	echo "GBM_MOBILITY = $(GBM_MOBILITY)" >> $@
+	echo "GBS_THRESHOLD = $(GBS_THRESHOLD)" >> $@
+	echo "NUCLEATION_EFFICIENCY = $(NUCLEATION_EFFICIENCY)" >> $@
+
+$(OUT_DIR)/$(SIM_NAME).py: $(SRC_DIR)/$(SIM_NAME).py
+	@echo "********** Copying python velocity callables..."
+	mkdir -p $(@D)
+	cp -f $< $@
+
+$(OUT_DIR)/$(SIM_NAME).flml: $(SRC_DIR)/$(SIM_NAME).flml
+	@echo "********** Copying serial flml file..."
+	mkdir -p $(@D)
+	cp -f $< $@
+
+$(OUT_DIR)/$(SIM_NAME).msh:
+	@echo "********** Building the mesh file..."
+	./envcheck.sh -m
+	pydrex-mesh -k="rectangle" -a xy $(WIDTH),$(DEPTH) \
+		-c $(HALF_WIDTH),-$(HALF_DEPTH) -r NW:$(RESOLUTION_HI),SE:$(RESOLUTION_LO) $@
+
+.PHONY: clean
+clean:
+	rm -rf $(OUT_DIR)