from .context import phaseflow

import fenics
import scipy.optimize as opt


def verify_pci_position(w):

    true_pci_pos = 0.075551957640682 # From Kai's MATLAB script
    
    def theta(x):
    
        wval = w(fenics.Point(x))
        
        return wval[2]
    
    pci_pos = opt.newton(theta, 0.075)
    
    assert(abs(pci_pos - true_pci_pos) < 0.002)
        
        
def test_stefan_problem():

    theta_h = 1.
    
    theta_c = -1.
    
    w = phaseflow.run(
        write_output = False,
        Pr = 1.,
        Ste = 1.,
        g = [0.],
        mesh = fenics.UnitIntervalMesh(1000),
        initial_values_expression = (
            "0.",
            "0.",
            "("+str(theta_h)+" - "+str(theta_c)+")*near(x[0],  0.) "+str(theta_c)),
        boundary_conditions = [
            {'subspace': 0, 'value_expression': [0.], 'degree': 3, 'location_expression': "near(x[0],  0.) | near(x[0],  1.)", 'method': "topological"},
            {'subspace': 2, 'value_expression': theta_h, 'degree': 2, 'location_expression': "near(x[0],  0.)", 'method': "topological"},
            {'subspace': 2, 'value_expression': theta_c, 'degree': 2, 'location_expression': "near(x[0],  1.)", 'method': "topological"}],
        regularization = {'a_s': 2., 'theta_s': 0.01, 'R_s': 0.005},
        final_time = 0.01,
        time_step_bounds = 0.001,
        linearize = False)

    verify_pci_position(w)
    

if __name__=='__main__':
    
    test_stefan_problem()

@@ -100,16 +100,26 @@
 osition(
+Ste, R_s, 
 w):%0D%0A%0D%0A 
@@ -121,16 +121,47 @@
 %0D%0A%0D%0A    
+data = %5B%0D%0A        %7B'Ste': 1., '
 true_pci
@@ -164,18 +164,21 @@
 _pci_pos
+':  
  
-=
  0.07555
@@ -187,16 +187,80 @@
 57640682
+%7D, %0D%0A        %7B'Ste': 0.1, 'true_pci_pos':    0.037826726426565%7D%5D
  # From 
@@ -428,24 +428,107 @@
 075)%0D%0A    %0D%0A
+    for record in data:%0D%0A    %0D%0A        if Ste == record%5B'Ste'%5D:%0D%0A        %0D%0A        
     assert(a
@@ -540,16 +540,24 @@
 i_pos - 
+record%5B'
 true_pci
@@ -564,17 +564,17 @@
 _pos
+'%5D
 ) %3C 
-0.002
+R_s
 )%0D%0A 
@@ -667,24 +667,84 @@
  -1.%0D%0A    %0D%0A
+    R_s = 0.005%0D%0A    %0D%0A    for Ste in %5B1., 0.1%5D:%0D%0A    %0D%0A    
     w = phas
@@ -755,16 +755,20 @@
 w.run(%0D%0A
+    
         
@@ -798,16 +798,20 @@
         
+    
 Pr = 1.,
@@ -824,19 +824,28 @@
     
+    
 Ste = 
-1.,%0D%0A
+Ste,%0D%0A    
     
@@ -855,24 +855,28 @@
  g = %5B0.%5D,%0D%0A
+    
         mesh
@@ -918,16 +918,20 @@
         
+  
+  
 initial_
@@ -957,39 +957,47 @@
  (%0D%0A            
+    
 %220.%22,%0D%0A
+    
             %220.%22
@@ -1011,16 +1011,20 @@
         
+  
+  
 %22(%22+str(
@@ -1096,16 +1096,20 @@
         
+    
 boundary
@@ -1121,24 +1121,28 @@
 itions = %5B%0D%0A
+    
             
@@ -1287,32 +1287,36 @@
 %7D,%0D%0A            
+ 
+   
 %7B'subspace': 2, 
@@ -1430,32 +1430,36 @@
 %7D,%0D%0A            
+   
+ 
 %7B'subspace': 2, 
@@ -1578,16 +1578,20 @@
         
+    
 regulari
@@ -1639,17 +1639,19 @@
 s': 
-0.005%7D,%0D%0A
+R_s%7D,%0D%0A    
     
@@ -1682,16 +1682,20 @@
         
+    
 time_ste
@@ -1721,16 +1721,20 @@
         
+    
 lineariz
@@ -1751,16 +1751,20 @@
 %0D%0A%0D%0A    
+   
+ 
 verify_p
@@ -1775,16 +1775,26 @@
 osition(
+Ste, R_s, 
 w)%0D%0A    

Test multiple Stefan numbers