I'm working on a synthetic setup for an experiment. It's based on the circular CalvingMIP geometry (https://github.com/JRowanJordan/CalvingMIP/wiki/Experiment-1), but I'm just using this as a starting point, my problem is not related to calving, moving front is disabled.
For my experiment, I need the thermal core enabled. But I'm getting unexpected results for the thermal steady state. The geometry of the ice sheet is rotational symmetric, as you can see in the plots.
I set the fields and BCs required by the thermal solver and then run the steady state thermal solution. I would expect the resulting temperature distribution to be symmetric as well. But as the plot shows, the result of basal temperature is severely asymmetric.
This is the script for the thermal solution, starting from loading the steady state calvingmip setup (I can supply the full calvingmip setup if needed):
md = load('saves/cmip-ex1-thermal-G20000-relax').md;
md = transientrestart(md);
T = 255;
md.basalforcings.geothermalflux = 0.05 * ones(size(md.mesh.x));
md.initialization.temperature = T * ones(size(md.mesh.x));
md.thermal.spctemperature = nan * ones(size(md.mesh.x));
md.thermal.spctemperature(logical(md.mesh.vertexonsurface)) = T;
md.thermal.isenthalpy = 1;
md.initialization.waterfraction = 0.0 * ones(size(md.mesh.x));
md.initialization.watercolumn = 0.0 * ones(size(md.mesh.x));
md.materials.rheology_n = 3 * ones(md.mesh.numberofelements, 1);
md.materials.rheology_B = paterson(T) * ones(size(md.mesh.x));
md.timestepping.start_time = 0;
md.timestepping.final_time = 0;
md.timestepping.time_step = 0;
md = solve(md, 'Thermal');
Am I missing something in the setup of the thermal solver? I tried to find any asymmetry in the setup, but haven't found anything, all fields are either constant in the whole domain or generated by symmetric functions.
I realize the resolution isn't great, and the asymmetry does improve with resolution. But it's still there, and in the same direction as before. In this plot you can also clearly see lines along the x and y axis, even though the mesh is irregular:
As a first idea I thought there might be a bias in the mesh or setup somehow. So I tried to rotate the geometry before running the thermal solver by adding the following lines to the script after transientrestart:
[md.mesh.x, md.mesh.y] = deal(-1 * md.mesh.y, md.mesh.x);
[md.mesh.x2d, md.mesh.y2d] = deal(-1 * md.mesh.y2d, md.mesh.x2d);
[md.initialization.vx, md.initialization.vy] = deal(-1 * md.initialization.vy, md.initialization.vx);
This should rotate the whole setup by 90 degrees counter clockwise. If the asymmetry is in the setup, the temperature asymmetry should now be rotated as well. But as this plot shows, the asymmetry remains the same. You can see that the geometry is rotated by the shape of the calving front:
I also tried to use a different mesh (squaremesh instead of triangle). I also tried the same with the Thule geometry of CalvingMIP, same result. The temperature is always colder towards the upper right corner (positive x and y direction).
This is why I think there might be a bias in the solver, but I realize that would be strange and I have no idea how that would happen or how to confirm it.
