update slides

Elmer/4b_PROGNOSTIC/1DVar_OUTPUT_BMB_.log.names

+File started at: 2020/11/20 12:06:45
+ 
+Elmer version: 8.4
+Elmer revision: 8860d158
+Elmer Compilation Date: 2020-11-10
+ 
+Variables in columns of matrix:1DVar_OUTPUT_BMB_.dat
+   1: Time
+   2: Volume
+   3: Volume Above Floatation
+   4: Volume rate of change
+   5: SMB Flux
+   6: BMB Flux
+   7: Residual Flux
+   8: Ice Discharge
+   9: Ice flux at Grounding Line
+  10: Grounded ice area
+  11: Floating ice area
+  12: Ice Free area

Slides/Slides.Rmd

@@ -542,7 +542,7 @@ Solver 3
     Nonlinear System Max Iterations = 10
     Nonlinear System Convergence Tolerance  = 1.0e-09
 !!
-    Stabilization Method = Stabilized
+    Stabilization Method = Stabilized  
     Transient Stabilisation = logical true
     Apply Dirichlet = Logical True
 
@@ -550,10 +550,30 @@ Solver 3
     Flow Solution Name = String "SSAVelocity"
 End
 ```
+
+---
+
+- stabilisation:
+  - The classical Galerkin method do not work to solve advection equation and a stabilisation method is required
+  - Stabilisation method is chosen with the keyword **Stabilization Method**
+    - **Stabilized** (default) is a Streamline Upwing Petrov-Galerkin method (SUPG)
+      - **Transient Stabilisation = logical true** use a different formulation for the stabilisation parameter
+    - **Bubbles** is a residual free bubbles method
+  - In general I set the bdf order to 2 for the timestepping
+   ```text
+    Simulation
+    ...
+      Timestepping Method = "bdf"
+      BDF Order = 2
+    ...
+    End
+    ```
+
 ---
 
 - Limiters: 
   - enforce min/max limits for the ice thickness
+  - make the system non-linear
   - Keywords related to limiters
   
 ```text
@@ -592,5 +612,126 @@ End
 
 - Boundary conditions
   - Dirichlet boundary condition are required at inflow boundaries
+  
+## Transient simulation
+
+The following solvers are solved in sequence during a time step
+
+  1. **Flotation**: Update geometry and mask for current $H$ 
+  
+  2. **SSA**: compute the advection field
+  
+  3. **Thickness Solver**: compute new thickness
+ 
+An iterative scheme can be used within time steps to couple the solvers; leading to a semi-implicit method, so this should improve the robustness 
+```text
+Simulation
+ ...
+  Steady State Min Iterations = 1
+  Steady State Max Iterations = 3
+ ...
+End
+```
+
+---
+
+- Tips:
+  - Save a **pvd** file that contains informations about physical time: 
+  ```text
+    Simulation
+      ...
+      Post File = "$name$.vtu"
+      vtu: Vtu Time Collection = Logical True
+      ...
+    End
+  ```
+
+## Remarks
+
+- For the **dSMB** and **BMB** experiments the mass balance evolves with time. To visualise the total mass balance we
+create an intermediate variable using the **UpdateExport** Solver
+```text
+Solver 1
+ Exec Solver = Before Timestep
+ Equation = "Update"
+ Procedure = "ElmerIceSolvers" "UpdateExport"
+ Variable = -nooutput "up"
+
+! This solver update exported variables from their definition in
+! the body forces
+ Exported Variable 1 = smba
+End
+```
+
+```text
+Body Force 1
+! update the total smb (see solver 1)
+  smba = Variable Time, smb, dsmb
+   real procedure "SMB_A" "SMB_ANOM"
+End
+```  
+
+---
+
+- For the **BMB** experiment we have written a user function to prescribe a quadratic dependence to zb
+```text
+Body Force 1
+  bmb = Variable GroundedMask,Zb
+   real procedure "SUB_MELT" "SUB_MELT"
+  SUB_MELT meltA = Real 0.0
+  SUB_MELT meltB = Real -1.25e-4
+End
+```  
+
+- The **bmb** variable is nodal and we prescribe melt only if the node is floting. If we directly apply this by default a partially floting element will have some melt resulting from the integration of the nodal variables. We can apply a no-melt parameterisation by checking if we are in a floting element or partially flotaing element:
+```text
+Body Force 1
+   Bottom Surface Accumulation = Variable bmb
+    real procedure "SUB_MELT" "DISC_MELT"
+End
+```  
+---
+
+```
+ FUNCTION DISC_MELT
+ ...
+  CurEl => Model % CurrentElement
+  NodeIndexes => CurEl % NodeIndexes(:)
+
+  GM => VariableGet(Model%Mesh%Variables,"groundedmask",UnFoundFatal=.TRUE.)
+  GMValues => GM % Values(:)
+  GMPerm => GM % Perm(:)
+
+  IF (ANY(GMValues(GMPerm(NodeIndexes(:))).GT.-0.5)) THEN
+    VarOut=0._dp
+  ELSE
+    VarOut = VarIn
+  END IF
+END FUNCTION DISC_MELT
+```   
+
+## Results
+
+- The solver **Scalar_OUTPUT.F90** in the **src** computes global variables (volume, fluxes, area): 
+
+<div class="centered">
+```{r, echo=FALSE, message=FALSE,warning=FALSE,results='hide', out.width="60%"}
+library(ggplot2)
+ctrl <- read.table("../Elmer/4b_PROGNOSTIC/1DVar_OUTPUT_CTRL_.log", header = FALSE)
+bmb <- read.table("../Elmer/4b_PROGNOSTIC/1DVar_OUTPUT_BMB_.log", header = FALSE)
+dsmb <- read.table("../Elmer/4b_PROGNOSTIC/1DVar_OUTPUT_dSMB_.log", header = FALSE)
+header <- c("Time","Volume", "Volume Above Floatation","Volume rate of change","SMB Flux","BMB Flux","Residual Flux","Ice Discharge","Ice flux at Grounding Line","Grounded ice area","Floating ice area","Ice Free area")
+names(ctrl) <- header
+names(bmb) <- header
+names(dsmb) <- header
+
+ggplot() + geom_line(data=ctrl,aes(x=Time,y=Volume*0.910/361.8e9,colour="ctrl")) +
+ geom_line(data=bmb,aes(x=Time,y=Volume*0.910/361.8e9,colour="bmb")) + geom_line(data=dsmb,aes(x=Time,y=Volume*0.910/361.8e9,colour="dsmb")) +
+  scale_color_manual(values = c('ctrl' = 'black','bmb'='blue','dsmb' = 'red')) +
+ theme_bw() + ylab("volume (mm sle)") + xlab("Time (a)") + labs(color = 'exp.')
+
+```
+</div>
+
 
 # inverse problems
\ No newline at end of file