{"action":"created","issue":{"url":"https://api.github.com/repos/pism/pism/issues/424","repository_url":"https://api.github.com/repos/pism/pism","labels_url":"https://api.github.com/repos/pism/pism/issues/424/labels{/name}","comments_url":"https://api.github.com/repos/pism/pism/issues/424/comments","events_url":"https://api.github.com/repos/pism/pism/issues/424/events","html_url":"https://github.com/pism/pism/issues/424","id":381649940,"node_id":"MDU6SXNzdWUzODE2NDk5NDA=","number":424,"title":"LC elastic model","user":{"login":"sebhinck","id":25295341,"node_id":"MDQ6VXNlcjI1Mjk1MzQx","avatar_url":"https://avatars0.githubusercontent.com/u/25295341?v=4","gravatar_id":"","url":"https://api.github.com/users/sebhinck","html_url":"https://github.com/sebhinck","followers_url":"https://api.github.com/users/sebhinck/followers","following_url":"https://api.github.com/users/sebhinck/following{/other_user}","gists_url":"https://api.github.com/users/sebhinck/gists{/gist_id}","starred_url":"https://api.github.com/users/sebhinck/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/sebhinck/subscriptions","organizations_url":"https://api.github.com/users/sebhinck/orgs","repos_url":"https://api.github.com/users/sebhinck/repos","events_url":"https://api.github.com/users/sebhinck/events{/privacy}","received_events_url":"https://api.github.com/users/sebhinck/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":10,"created_at":"2018-11-16T15:38:29Z","updated_at":"2018-12-19T09:55:21Z","closed_at":null,"author_association":"NONE","body":"When enabling the LC elastic model I found strange artifacts in the uplift rate (see Figure) that seem to originate from some kind of unsymmetrical problem.\r\n![uplift_lc_elastic](https://user-images.githubusercontent.com/25295341/48630786-5f31f380-e9bd-11e8-8aba-102986942ccc.png)\r\n\r\nThis is why I had a closer look on how the elastic deformation is calculated. According to Bueler et al 2007 it is calculated as the convolution of the load with the Green's function, which is calculated once at the initialization of the LingleClarkSerial model (`m_load_response_matrix` line 218ff `src/earth/LingleClarkSerial.cc`). This matrix (see figure below) is only the upper right quadrant of the radially symmetrical Green's function (with the point load located in the origin). This matrix contains all information needed, but then the convolution function needs to take care of this.\r\n\r\n![load_response_matrix](https://user-images.githubusercontent.com/25295341/48630826-7b359500-e9bd-11e8-9a34-ee0f376effe7.png)\r\n\r\n\r\n`conv2_same` in `src/earth/matlablike.cc` line 26ff does only integrate the product of Green's function and load over the south-western (lower left) quadrant and thus neglects contributions of masses in the other quadrants. One simple but very computational demanding solution would be:\r\n```\r\nvoid conv2_same(Vec vA, int mA, int nA,  Vec vB, int mB, int nB,\r\n                Vec vresult) {\r\n\r\n  pism::petsc::VecArray2D\r\n    A(vA, mA, nA),\r\n    B(vB, mB, nB),\r\n    result(vresult, mA, nA);\r\n\r\n  for (int j=0; j < nA; j++) {\r\n    for (int i=0; i < mA; i++) {\r\n      double sum = 0.0;\r\n\r\n      for (int r = 0; r < mA; r++) {\r\n        for (int s = 0; s < nA; s++) {\r\n          sum += A(r, s) * B(std::abs(i - r), std::abs(j - s));\r\n        }\r\n      }\r\n      result(i, j) = sum;\r\n    }\r\n  }\r\n}\r\n```\r\nCompare the elastic deformation fields with and without this modification. Note the deformation being zero south and west of the ice sheet, which it shouldn't.\r\n\r\n![elastic_compare](https://user-images.githubusercontent.com/25295341/48630853-87b9ed80-e9bd-11e8-9c01-4bc39c6b87fe.png)\r\n\r\n\r\nI think there exist more elaborate functions to calculate the convolution based on FFT which are more efficient than this one here, and which is also mentioned in Bueler et al 2007.\r\n\r\nPlease correct me if I misunderstood something and I am wrong.\r\n\r\nI also have a question regarding #183 (make the on/off state of the elastic part of the \"-bed_def lc\" earth model configurable #183). What is the reason for disabling this part of the model by default though it is an Elastic Lithosphere Relaxing Asthenosphere (ELRA) model. This should probably be mentioned more explicitly in the user manual.\r\n\r\n"},"comment":{"url":"https://api.github.com/repos/pism/pism/issues/comments/448536177","html_url":"https://github.com/pism/pism/issues/424#issuecomment-448536177","issue_url":"https://api.github.com/repos/pism/pism/issues/424","id":448536177,"node_id":"MDEyOklzc3VlQ29tbWVudDQ0ODUzNjE3Nw==","user":{"login":"flapo099","id":26421314,"node_id":"MDQ6VXNlcjI2NDIxMzE0","avatar_url":"https://avatars2.githubusercontent.com/u/26421314?v=4","gravatar_id":"","url":"https://api.github.com/users/flapo099","html_url":"https://github.com/flapo099","followers_url":"https://api.github.com/users/flapo099/followers","following_url":"https://api.github.com/users/flapo099/following{/other_user}","gists_url":"https://api.github.com/users/flapo099/gists{/gist_id}","starred_url":"https://api.github.com/users/flapo099/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/flapo099/subscriptions","organizations_url":"https://api.github.com/users/flapo099/orgs","repos_url":"https://api.github.com/users/flapo099/repos","events_url":"https://api.github.com/users/flapo099/events{/privacy}","received_events_url":"https://api.github.com/users/flapo099/received_events","type":"User","site_admin":false},"created_at":"2018-12-19T09:55:21Z","updated_at":"2018-12-19T09:55:21Z","author_association":"NONE","body":"Thank you very much for this issue being reported. We are using heavily the -bed_def lc as we are looking into the self-graviational effect of ice and land on sea level. This fix for the v1.1 will be enough for us as we are using the new 2D sea level fonction."}}
{"actor":{"display_login":"flapo099"}}